Amandine, Daniel, Elisabeth Merci À Aleix

Giacomo

Lucie

Mounib, Marion Nancy Maria

Piernicola, Rituparna, Stéphane, Suvadip, Thomas Vannary Pascal

Keywords: Neuroscience, in vivo neural networks, neural ensembles, spikes, inverse problem, coupling spatial statistics, Bayesian . . . . . . . . . . . .

The neural code, i.e. how interconnected neurons can perform complex operations, allowing the quick adaptation of animals to their environment, remains an open question and an intensive field of research both in experimental and computational neurosciences. Advances in molecular biology and microscopy have recently made it possible to monitor the activity of individual neurons in living animals and, in the case of small animals containing only a few thousands of neurons, to measure the activity of the entire nervous system. However, the mathematical framework that would bridge the gap between single neuron activity and the emergent computational properties of neuronal ensembles is missing. In the following thesis manuscript, we introduce a sequential statistical processing pipeline that efficiently and robustly extracts neuronal ensembles from calcium imagery of neuronal activity. In particular, we develop a Bayesian inference framework based on a biologically interpretable model to extract neuronal ensembles characterized by noise, asynchrony and overlapping. The provided tool demonstrates that a Gibbs sampling routine can efficiently estimate statistical parameters and hidden variables to uncover neuronal ensembles based on synchronization patterns both on synthetic data and on various experimental datasets from mice and zebrafish visual cortex to Hydra vulgaris. The thesis equally develops a point process statistical framework to quantify how neuronal ensembles encode evoked stimuli or spontaneous behaviors in living animals. This versatile tool is also used for the inference of the functional connectivity of neuronal activity or the automatically calibration procedure of the spike inference algorithms applied to calcium recordings. For the providing algorithms to be largely spread in the neurobiologist community, results are supported by interpretable biological estimates, statistical evidence, rigorous mathematical proofs, and free-available software. Our contributive implementation, that goes from pixel intensity to estimated neuronal ensembles, equally identify from the synchronous firing patterns of neuronal ensembles, neurons with specific roles that can be used to predict, improve, or alter the behaviors of living animals. The provided framework unravels the emergence of collective properties from the recording of extremely varying individual signals that make the neural code still elusive.

Résumé : Le code neuronal, c'est-à-dire la manière dont les neurones interconnectés peuvent effectuer des opérations complexes, permettant l'adaptation rapide des animaux à leur environnement, reste une question ouverte et un champ de recherche intensif tant en neurosciences expérimentales qu'en neurosciences computationnelles. Les progrès de la biologie moléculaire et de la microscopie ont récemment permis de surveiller l'activité de neurones individuels chez un animal vivant et, dans le cas de petits animaux ne contenant que quelques milliers de neurones, de mesurer l'activité de l'ensemble du système nerveux. Cependant, le cadre mathématique qui permettrait de combler le fossé entre l'activité d'un seul neurone et les propriétés computationnelles émergentes des ensembles neuronaux fait défaut. Dans le manuscrit de thèse suivant, nous présentons un pipeline de traitement statistique séquentiel qui permet d'extraire efficacement et de manière robuste des ensembles neuronaux à partir de l'imagerie calcique de l'activité neuronale. En particulier, nous développons un cadre d'inférence bayésienne basé sur un modèle biologiquement interprétable pour extraire des ensembles neuronaux caractérisés par du bruit, de l'asynchronisme et du recouvrement. L'outil fourni démontre qu'une procédure d'échantillonnage de Gibbs peut estimer efficacement les paramètres statistiques et les variables latentes pour extraire les ensembles neuronaux basés sur un modèle de synchronisation à la fois sur des données synthétiques et sur des données expérimentales allant de stimulations du cortex visuel de la souris et du poisson zèbre à l'activité spontanée de Hydra vulgaris. La thèse développe également un cadre statistique de processus ponctuel pour quantifier la façon dont les ensembles neuronaux encodent les stimuli évoqués ou les comportements spontanés chez les animaux vivants. Cet outil polyvalent est également utilisé pour l'inférence de la connectivité fonctionnelle de l'activité neuronale ou la procédure de calibration automatique des algorithmes d'inférence de pics appliqués aux enregistrements calciques. Pour que les algorithmes fournis soient largement diffusés dans la communauté des neurobiologistes, les résultats doivent être étayés par des estimations biologiques interprétables, des preuves statistiques, des démonstrations mathématiques rigoureuses et des logiciels libres d'accès. Notre implémentation contributive, qui va de l'intensité des pixels aux ensembles neuronaux estimés, identifie également, à partir des schémas d'activation synchrone des ensembles neuronaux, les neurones ayant des rôles spécifiques qui peuvent être utilisés pour prédire, améliorer ou modifier les comportements d'animaux vivants. Le cadre fourni permet de démontrer l'émergence de propriétés collectives à partir de l'enregistrement de signaux individuels extrêmement variables, qui rendent le code neuronal encore insaisissable.

Mots-clés : Neurosciences, Réseaux de neurones in vivo, ensembles neuronaux, potentiels d'action, problème inverse, statistiques spatiales de couplage, inférence bayésienne, stimulation et comportements accueil régulier à Paris et la cuisine de bon plats ivoiriens, ont fait preuve d'une gentillesse et bienveillance familiale qui m'ont énormément aidé en cette fin de thèse.

Enfin un très grand merci à mon grand-père, Michel pour son affection et son soutien au cours de la thèse, et qui à la lecture de ses remerciements ne pourra s'empêcher de me céder la Sirocco.

Je remercie également ma belle-famille pour leur enthousiasme et leur accueil pendant ce travail. Merci à Pascale et Jean-Marc, Ambre, Roxane et Xavier.

Je remercie également Julien et sa famille, Amélie, Anaé et même, oui, même Gabriel pour leur accueil et leur gentillesse au cours de cette thèse.

Que serait notre famille de sang sans notre famille de coeur ? Je remercie infiniment mes amis d'avoir été un soutien indéfectible pendant cette thèse. En particulier, un grand merci à mes frères Christopher et Gabriel (ordre alphabétique). Je vous remercie d'avoir toujours été présents, parfois acerbes, parfois jugeants, humiliants ou moqueurs mais toujours là pour me remonter le moral, et me donner le sourire avec une rigueur et un dévouement à toute épreuve. Je vous remercie du fond du coeur pour cette longue aventure échelonnée par des sorties dans les mines galactiques, les voyages du bout du monde au fin fond des Cévennes ou des Vosges, ces heures de jeux de société, d'escape Game ou de randonnées. Cette thèse vous est également dédiée.

Un grand merci à mon amie de toujours Inès, qui par son accueil, son âme chaleureuse bienveillante et calme m'a également apporté du réconfort. Son oreille attentive m'était précieuse dans les moments de doutes que j'ai traversé par le biais de nos repas et de nos balades. Merci pour ton aide et ta douceur.

A ma conjointe Maëlle. . . . Je te remercie sincèrement du fond du coeur pour ta bienveillance, ta patience et ta générosité. Merci d'avoir toujours cru en moi, d'avoir été ma première supportrice, de m'avoir écouté pendant des heures, des mois et même des années durant. Je te remercie d'avoir toujours eu foi alors que rien ne marchait, d'avoir cru alors que moi-même je n'espérais plus. Tu as écouté des explications scientifiques pendant des centaines heures, discuté avec moi les analogies, accepté de contempler, parfois sous la contrainte, mes images ou mes résultats. Merci également pour tout le reste. De m'avoir permis de m'aérer, penser à d'autres choses, bref de m'avoir permis de continuer à vivre à côté.

II

Chapter I

Overview 1 Context i Summary

Neurons are the basic computing units of the brain. However, the neural code, i.e. how interconnected neurons can perform complex operations, allowing the quick adaptation of animals to their environment, remains an open question and an intensive field of research in experimental and computational neurosciences. A promising approach to "break" the neural code, i.e. to relate the sequential activity of neurons or ensemble of neurons to animal's behavior, consists in labeling individual neurons with fluorescent indicator (genetically-encoded calcium indicator typically), and following neuron's activity in small model organisms such as the zebrafish [1] or Hydra [2]. Such organisms, with their relatively simple nervous organization and reduced repertoire of behaviors, can be systematically analyzed and perhaps completely understood. While state-of-the art genetic editing and time-lapse microscopy allow the long-term imaging of individual cells in freely-behaving animals, the mathematical paradigm, that would bridge the gap between single cell activity and the emergent computational properties of neuronal ensembles, is missing. The main objective of this thesis is the development of an innovative statistical framework to extract, from individual neuronal fluorescent activity, emergent computational ensembles whose activity is correlated with external stimulations or spontaneous behaviors. This will provide important insights into the computational strategy of interacting neurons, unraveling the mathematical principles that underly the robustness and efficiency of brain computation.

ii Goals

The goal of this thesis is to develop statistical tools to characterize neuronal ensembles from individual calcium fluorescence traces of neurons in free-behaving animals. Decoding neuronal ensembles from the recording of neuron calcium individual traces corresponds to a multi-step processing pipeline (see Figure I.1). The interdependent steps are divided as follows :

1) The extraction of neuron individual fluorescence traces from the calcium image recordings, 2) The inference of the exact time of action potential emissions, 3) The estimation of the network topology, 4) The detection of the neuronal ensembles, 5) The inference of computation and interaction rules between ensembles, 6) The prediction of the animal behavioral states from ensemble activity. Each step requires the development of robust tools to account for the risk of propagating errors. For example, erroneous tracking or spike deconvolution will decrease the accuracy of the estimation of neuronal ensembles.

Contributions of the thesis i Contributions

Our work has been structured around the different steps proposed by the original processing pipeline (see Figure I.1). In this work, we demonstrate the capabilities of our tools on several model animals : the mice visual cortex in response to stimulations (data provided by Yuste's lab), the larval optical tectum of zebrafish (data in free-access from [3]) and the Hydra vulgaris neural activity in response to free behaviors (data from [2]). The contributions are structured as follows around the different steps of pipeline:

1)
Step 2 -Inferring spikes from fluorescence traces: The first contribution is the benchmarking of the most promising state-of-the-art spike inference techniques on a calcium fluorescent simulator to estimate the most efficient method on data whose properties match Hydra vulgaris' firing statistics and non-linearities. The proposed versatile simulator embeds non-linear dynamics, non-homogeneous Poisson firing statistics with several firing regimes. The simulator aims to assess the robustness and efficiency of the methods applied to the Hydra vulgaris case, and to point out potential limitations.

2) Step 2 & 3 -Inferring spikes from fluorescence traces & The estimation of a network topology:

The second contribution is the development of a statistical tool to infer neuron-to-neuron statistical coupling in order to extract the network topology while being robust to false positives. Its use as a functional connectivity method facilitates spike inference by filtering out neuron activation artifacts.

3) Step 2 -Inferring spikes from fluorescence traces: The third contribution of this work is the application of spatial statistics framework to filter out spike artifacts in the presence of a decorrelated noise. 4) Directly from step 2 to step 4 -Neuronal ensemble detection: The fourth contribution aims to extend existing tools for estimating functional neuronal ensembles by using a Bayesian Inference framework. It would allow analyzing the emergent synchronization properties at network scale from the activation of individual cells. This framework uses synchronicity as a biological interpretable neuronal ensemble definition, handling overlapping, and embeds the point process statistical coupling as an evaluation of assemblies relevance in response to stimulations. The spatial statistics framework is adapted to measure the correlation between stimulations and neuronal ensemble activity.

ii Publications

The publications associated to this work are:

• The benchmarking of the existing Spike Inference Techniques for calcium fluorescent traces so as to find the most relevant ones to apply in the Hydra vulgaris case where statistical firing properties, noise level, non-linear baseline dynamics largely pollute the data.

S. Kubler • The implementation of an overlapping neural ensemble detection algorithm by Bayesian Inference with an explicit definition of neural ensembles defined as a synchronization process. The statistical coupling framework can be adapted to infer the coupling between neuron ensemble activities and stimulations or behaviors.

S. Kubler, J.-C. Olivo-Marin, et T. Lagache, " Bayesian Inference of overlapping neuronal ensembles ", Nature Computational Science, submitted, under-review.

Experimental dataset of in vivo neural network analysis

In our work, we used different in vivo experimental setups detailed hereafter.

i Mice visual cortex dataset from Yuste's laboratory

The mice brain visual cortex V1 of a transgenic animal is imaged using a two-photon microscope. Static grating stimuli with different directions are applied in front of the mice to scan its receptive field. Group of neurons of the visual cortex are expected to fire synchronously in response to the stimuli. [4] [5] A grey screen is assumed to be the resting state to record the spontaneous neural activity. The neuron spiking activity matrix is extracted during the whole experiment. Its aim is to cluster the neurons based on their synchronous neural activity and to correlate the group existence with stimuli. ii Zebrafish optic tectum dataset from Mölter et al [3] The experiment is similar for the larval zebrafish tectum, Danio rerio, stimulated using visual stimuli of spots displayed at specific locations. These locations are 15°spaced leading to several fluorescence responses. The animal is genetically-engineered to express elavl3:H2B-GCaMP6s and neuronal ensemble responses have been demonstrated.

iii Hydra vulgaris' spontaneous behavior dataset from [2] Hydra vulgaris is a freshwater polyp animal that belongs to the Cnidarian class. This tubular species has the simplest functional nervous system of the entire animal kingdom that supports though behavioral advanced functions. Its nervous system is decentralized within the whole animal, with neurobiological features very similar to mammalian species [6] because of the existence of neuron, sensory and ganglion cells, neurites, chemical [7] and electrical synapses (gap junctions) [8], or neurotransmitters and neuromodulators (neuropeptides) [9]. b) The neural response in the optic tectum is recorded using a 2-photon microscope. A wave activation reponse is observed for each stimuli but the strength of the response and the overlapping largely vary according to the stimulation. c) A very clear spike response can be recorded for some specific stimulation.

All these elements allow Hydra to react to light variation and vibrations [10], temperature [11], chemical signaling [6] and interact with its environment to perform coordinated movements, catch its preys, or regulate its osmolarity [12].

The animal is small, transparent, can be easily genetically engineered by creating a transgenic Hydra line that expresses the fluorescent calcium indicator GCaMP6s making the microscopy of its fluorescent neurons straightforward.

The nervous system of Hydra is composed by few hundreds to few thousands neurons distributed inside two independent nerve nets, in the ectoderm and endoderm with a radial symmetry. Its neural circuitry is divided into 4 different functional networks that coordinate a limited mechanosensory behavior repertoire [13] [14]. Rhythmic Potential 1 (RP1), Rhythmic Potential 2 (RP2), Contraction Burst (CB) and Subtentacle Network (STN) encode respectively elongations, radial contractions, longitudinal contraction, and nodding [2]. It is expected that sequential activation patterns of such networks could encode behavioral tasks like feeding, response to light, temperature or starvation, somersaulting and egestion.

Hydra vulgaris is genetically engineered to express Green-Fluorescent Protein (GFP) by a micro-injection of plasmids in fertilized eggs. The adult transgenic free-behaving animal is placed in a restricted environment between two coverslips separated by 100 µm to image its neural activity during spontaneous activity at 10-33 Hz.

Even if RP1, RP2 and CB can have coactive neurons firing in the same 100ms window, Hydra's nervous system is divided into networks that are non-overlapping structurally and functionally [2]. It means that individual neurons only participate in a unique network at a time. Networks encode specific behaviors : RP1 and RP2 are respectively located in the Each part of the animal plays a fundamental role to interact with its environment [15]. c) To make the microscopy more straightforward the animal is constrained between two coverslips to place it entirely in the field of view of the microscope [2].

ectoderm and endoderm associated with elongation and radial contraction of the animal. Changes in their activation frequency were associated with the behavior triggering. Slow potential propagation that occurs simultaneously to asymmetric longitudinal contraction of the body column were reported during the activation of STN. STN can either be initiated from the body column to the tentacles or from the hypostome region traveling downward. Multiple conduction speeds have equally been observed in Hydra in vivo recordings in CB network [2]. Fast propagation seems related to longitudinal contraction while slow propagation could just increase the excitability of ganglion and sensory neurons. The neural network findings of Hydra highlight many similarities with other animals of this phylum. The existence of nerve nets, the spontaneous longitudinal contractions and characteristic sequence in feeding behavior are part of them.

The thesis manuscript is divided as follows. In a first part, an introduction describes the convergence of neuroscience research toward neuronal communication at the population cell level in the in vivo neural networks of entire free-behaving animals. In the second part, a literature review is carried out to go from individual cell calcium recordings to the detection of neuronal ensembles that encode behaviors. In a third part, the general issues of the processing pipeline to free-behaving animal case are highlighted and the contributions of the manuscript reminded. In the fourth part, the contributions of the manuscript are concatenated and discussed in part five. Finally, we conclude on the research topic and provide future prospects to go beyond in the interpretation of the neural code.

Introduction 1 Neuronal network i From Brain to neuronal network

The brain is an inter-individual varying complex structure that embeds a huge diversity of cognitive functions : memory, thinking, feeling, dreaming, learning, moving, decisionmaking are all performed by the brain, supported by the same local chemical and electrical transmission machinery. According to its development and environment, its structures can evolve, adapt, restructure dynamically in a structural and functional way. In neuroscience, the establishment of theoretical principles that can be easily transposed to different cortical areas, different individuals or animal models, or different time and space scales, is complicated and often requires the statistical analysis of simplified models, studied at lower scale, in experiments where the number of study parameters is deliberately limited. Indeed, brain is a complex organ that far exceeds our technological capacity to finely analyze the interactions of around 86 billion cells (∼ 10 11 ) connected by more than 100 trillion (∼ 1.10 20 ) connections [16] [17]. In neurosciences, the question of how information is encoded by interacting individual neurons remains unknown.

Understanding the brain, usually, means defining functional properties of the cognition processes in relation to behaviors, unveiling the intrinsic algorithms of brain computation and being able to explain how neurons and their connections support such algorithm execution [18]. Brain computation is supported by the interaction of up to billions of neurons [16] [17] in mammals making its understanding complex and elusive. Changing scale by going from entire brain to single neurons analysis (see Figure II.1) makes sense since it is assumed that common general scale-free working principles exist. It means that neural properties could be more easily extracted at a lower scale and transposable to a higher one. This changing scale analysis is called a bottom-up approach. Such approaches have unveiled similar structural and organizing principles between in vivo neuronal networks, corresponding to few (hundreds to thousands) interacting neurons, and brain. The communication and information exchanges are performed through hierarchical functional and anatomical modules related to specific tasks involving a huge diversity of molecular, cellular, and neuronal phenomena [19]. Hierarchy and modularity are the two cornerstones that support robust, compartmentalized, complex invariant-scale architecture able to communicate effectively. It shapes the network dynamics and demonstrates the close relationship between topology and dynamics, structure, and neural function. Deciphering the neural code would be a great progress biologically and technologically speaking. The first application would be the direct use of such knowledge in medicine to some neuropathological diseases. Indeed, corroborating studies highlight that topological network properties at brain level are correlated with cognitive abilities [21] like verbal fluency [22], IQ [23] or working memory [24] accuracy for instance and that such structure alterations could lead to diverse clinical states as schizophrenia [25] [26] [27] [28] or Alzheimer's disease [29] [30] [START_REF] Fu | Different Functional Network Connectivity Patterns in Epilepsy: A Rest-State fMRI Study on Mesial Temporal Lobe Epilepsy and Benign Epilepsy With Centrotemporal Spike[END_REF] [32] [START_REF] Wang | Alterations in Brain Network Topology and Structural-Functional Connectome Coupling Relate to Cognitive Impairment[END_REF]. Interestingly, at the neuronal network level, the synchronization of firing patterns in neural networks could equally be linked to epilepsy [START_REF] Fu | Different Functional Network Connectivity Patterns in Epilepsy: A Rest-State fMRI Study on Mesial Temporal Lobe Epilepsy and Benign Epilepsy With Centrotemporal Spike[END_REF] [START_REF] Klugah-Brown | Altered Dynamic Functional Network Connectivity in Frontal Lobe Epilepsy[END_REF]. Understanding the way in vivo neural networks finely work could open a way of new therapeutics.

ii Imaging neuronal network in vivo

Gold-standard techniques of imaging neuronal networks in vivo target brain regions like visual cortex in mice (see Figure II.2). Recent progress of fluorescent microscopy and genetically engineering enabled the imagery of a whole population of neurons at individual cell level in simpler animal models with limited repertoire of behaviors. Such imaging of neuron population is easier since it decreases the level of complexity and simplifies the way to record neural activity. Imaging the simultaneous activity of entire nervous systems in simple model organisms such that the worm [START_REF] Venkatachalam | Pan-neuronal imaging in roaming Caenorhabditis elegans[END_REF] [START_REF] Nguyen | Whole-brain calcium imaging with cellular resolution in freely behaving Caenorhabditis elegans[END_REF] [37] [START_REF] Schrödel | Brain-wide 3D imaging of neuronal activity in Caenorhabditis elegans with sculpted light[END_REF], fruit fly [START_REF] Chhetri | Whole-animal functional and developmental imaging with isotropic spatial resolution[END_REF] [40] and zebrafish [START_REF] Ahrens | Brain-wide neuronal dynamics during motor adaptation in zebrafish[END_REF] [42] [43] [44] [45] [START_REF] Tomer | SPED Light Sheet Microscopy: Fast Mapping of Biological System Structure and Function[END_REF] has benefited from breakthrough in microscopy, probe development and genetic engineering. In particular, the recent advances in fluorescent calcium microscopy techniques make it possible to image the activity of a whole model animal at the scale of individual neuron cells. In this thesis, we will focus on Hydra vulgaris, a small fresh water cnidarian, that has been recently genetically engineered to express calcium fluorescent indicators at single neuron level [2].

Two-Photon Microscope

Mice Visual Cortex (V1)

Hydra Vulgaris

Confocal Spinnning

Disk Microscope

Two-Photon Microscope

Larval Zebrafish (OT) ) or larval zebrafish (data from [3], image from [START_REF] Diana | Bayesian inference of neuronal assemblies[END_REF]) using GCamP calcium imaging.

It is worth noting that Hydra vulgaris is very promising for the analysis of in vivo neural networks in free-behaving animal. Indeed, this organism, with a decentralized nervous system possesses a limitless regenerative capacity and its transparency allows the complete imaging of entire nervous system [START_REF] Han | Comprehensive machine learning analysis of Hydra behavior reveals a stable basal behavioral repertoire[END_REF]. In addition, the neuronal substrates of its limited behavioral repertoire might be completely understood in a near future. If we highlight its ability to regenerate from an aggregate of single cells [START_REF] Lovas | Ensemble synchronization in the reassembly of Hydra's nervous system[END_REF], it also appears promising for understanding neuronal development from a functional and structural point of view.

Neuronal imaging in unrestrained animals has expanded the range of behaviors analyzable at micro-circuitry level in several model organisms [START_REF]Method of the Year 2018: Imaging in freely behaving animals[END_REF]. Traditionally, animals are constrained to analyze the neural response in accordance with restricted behaviors to deliberately limit the study parameters and give to experimenters exquisite control over the sensory environment of their studies e.g, visual cortex of mice which responds to pre-defined stimuli (patterns). However, such procedures restrict animal ability to exhibit more complex behaviors such as social interaction, motion, biophysical regulation. The recent interest of neural code analysis in free-behaving animal has equally been supported by new microscopy modalities like light-sheet microscopy. The combined neural activity recording with behavioral monitoring may uncover correlated effects able to predict and anticipate how neural circuitry encode a decision and a motion [START_REF] Wallace | Circuit interrogation in freely moving animals[END_REF] [52] [START_REF] Pereira | Fast animal pose estimation using deep neural networks[END_REF].

From single neurons to population coding i Neurons

The pioneer of modern neuroscience is S. Ramón y Cajal who was the first around 1900, to discover and sketch new cell type with specific body shape and long-wire extension called neurons [54] [55]. This cell type extremely dense in the nervous central system of vertebrates (10 4 cells for kilometers of wiring per cubic millimeter at human cortex) was claimed to be the functional unit of brain computation. Functionally, the entity is divided into distinct parts called dendrites, soma, and axon that respectively play a role of input signal collector, central processing unit and output transmission channel. ii

Conveying information between neurons

The neuron activity is encoded through short electrical pulses called action potentials or spikes that typically last between 1-2 ms with an amplitude about 100 mV. They correspond to electrical depolarizations of membrane potential that spread along somas, dendrites or axons and allow distant neurons to communicate with each other. The underlying mechanisms to transmit these pulses vary according to the nature of the neuron involved, and the biological local topology and structure. A recovery delay called refractory period follows a spike event to allow the biological transmitter reservoirs to recharge and prevent two spikes from being emitted too closely. Neuroscientists generally agree that action potential shape does not vary much and does not, thus, encode any information. Major theories have have been developed to describe the way the information could be embedded in such spiking sequences called spike trains (see Figure II.5) corresponding to the sequential activation pattern of a single neuron.

Information is transmitted from one neuron to another through a junction called synapse (see Figure II.6). The afferent neuron is called pre-synaptic cell and the receiving neuron is called post-synaptic cell. The gap between the axon of a presynaptic neuron and the dendrite or soma of a post-synaptic cell constitutes the synapse. The most common type of synapse in vertebrate brain is the chemical synapse meaning that information will be conveyed by the release of chemical neurotransmitters through the synaptic cleft. The presynaptic membrane depolarization linked to an action potential arrival triggers a complex chain of biochemical signaling in pre-synaptic vesicles. Neurotransmitters are released from the presynaptic vesicles in the synaptic cleft and bind to post-synaptic receptors, triggering the opening of ions b) also exist to transmit information directly between two close cellular bodies through ion channels. Both chemical and electrical gap junctions have been observed in Hydra vulgaris [7] [8]. It is worth mentioning that other molecules called neuromodulators regulate the neuronal communication such as large family of neuropeptides. These molecules excite or inhibit the neuronal activity by diffusing in entire areas without targeting a specific site. In particular, these peptides seem to be the main players in ancient animals like Hydra vulgaris [9] where they could allow a differentiation of neural microcircuitry and be the intrinsic mechanism of neuronal communication. Finally, action potential has a stereotyped shape and transmit at a constant speed involving a self-maintained communication process that requires the support of external agents. One of them are the glial cells that support and regulate neuronal transmission.

For a long while, the neuroscience community exhaustively studied neuronal transmission by estimating the time variation of the membrane potential and creating generative mathematical models able to account for different levels of complexity and generate dif- In parallel, easier mathematical frameworks were applied at mesoscopic scale to describe neuronal information transmission. The mathematical model of Leaky-Integrate-and-Fire (LIF) neuron provides, thus, a representation of neuronal transmission as an integration of inputs via a summation and a thresholding process that support the representation of the neural networks of artificial intelligence. By doing such analysis, it is possible to informatically simulate data of entire neuronal populations with excitatory and inhibitory effects. The direct applicability of this kind of analysis to the in vivo case is, however, limited by: 1) the sampling frequency that does not allow the recording at an action potential resolution. (sampling frequency are around 10 Hz for calcium fluorescent microscopy leading to a 100-millisecond recording period against a 10 millisecond action potential dynamics) 2) the lack of an exact anatomical description of an entire neural network due to poor spatial resolution. Since Adrian's early observations, several theories have been developed to model how neurons can encode information. Rate coding or frequency coding claims that information is encoded by action potential frequency [START_REF] Adrian | The impulses produced by sensory nerve endings: Part 3. Impulses set up by Touch and Pressure[END_REF]. This theory accounts for the heterogeneity and diversity of action potential responses of a neuron given a stimulus. When intensity of the stimulus increases, action potential rate generally increases non-linearly as well. Thus, an estimation of firing rate or instantaneous firing rate function through the spike density are usually performed using peristimulus time histogram or kernel density estimation techniques [START_REF] Shinomoto | Estimating the Firing Rate[END_REF] [62]. They extract correlation activity patterns through rate coding. This encoding model was limited since unable to explain experiments showing that changes in correlation was independent of spike rate [START_REF] Salinas | Correlated neuronal activity and the flow of neural information[END_REF], experiments where neurons timing and duration of nonfiring neuron was also carrying information [START_REF] Forrest | Intracellular calcium dynamics permit a Purkinje neuron model to perform toggle and gain computations upon its inputs[END_REF] or experiments where low frequency events were conveying information like for Purkinje neurons [START_REF] Abrams | Signals and circuits in the Purkinje neuron[END_REF] [66] [START_REF] Popa | Purkinje Cell Representations of Behavior: Diary of a Busy Neuron[END_REF]. More generally rate coding fails to explain fast organism response especially when the response is faster than 100 milliseconds or with only few spikes involved [START_REF] Dear | A possible neuronal basis for representation of acoustic scenes in auditory cortex of the big brown bat[END_REF] [69] [START_REF] Thorpe | Speed of processing in the human visual system[END_REF].

The limitation of the rate coding to account for neuroscience experiments with fast dynamics led to another coding theory called time coding or spike coding. In this theory, the temporal pattern of precise action potential emission is thought to convey information, especially in sensitive neurons where the response is fast. While rate coding, basically, counts and averages the number of spikes that fall into a time bin window, time coding asserts that the exact date of emission action potentials and high-order statistics, like inter-spike intervals, offer a higher dimensionality to encode information. Timescale of both theories and the statistical tools to model the processes are different. A spike train, thus, owns a better encoding ability with time coding and the diversity of responses to a same stimulus could represent fundamental information.

However, these classical coding theories are limited to single neurons while information is known to be processed by population of neurons.

iii From single neurons to neural networks

Synaptic plasticity

Neurons are densely connected with synapses, forming functional ensembles. Each synapse has an individual efficacy that describes the relative amplitude of the neuronal post-synaptic response (EPSP) related to the presynaptic afferent action potential. Electrophysiological experiments demonstrated the ability to alter the synaptic strength either positively through a process known as Long-Term Potentiation (LTP) or negatively via the Long Term Depression (LTD). This dynamic change of synaptic strength with respect to activity is called synaptic plasticity. This process would be the fundamental support of learning and memory at neuronal level. The changes can be imprinted at long-or short-time scales and demonstrate neuron's ability to modify the synaptic efficacy based on their functional activity. The main biophysical mechanisms that account for synaptic plasticity are the modulation of the number of post-synaptic neurotransmitter receptors and the restructuring of the shape of the receiving neuron's dendrite. LTP corresponds to the increased synaptic strength over hours when injecting high-frequency pulses in-presynaptic fibers to evoke a post-synaptic firing response. Another experimental protocol inducing a prolonged increase of synaptic strength is the Spike-Timing-Dependant-Plasticity (STDP). It corresponds to an increase of the synaptic efficacy in a neural network by activating groups of neurons with a correct time delay. It induces an artificial causality that increases synaptic coupling efficacy on long time-periods. Synaptic plasticity underlies the emergence of neuronal networks, i.e. groups of highly connected neurons that are thought to encode and convey information.

Neuronal network and Hebbian rule

Over time, synaptic coupling parameters evolve, are adjusted and optimized in a process called learning (see Figure II.7). Several theories to describe mathematically this learning rule have been formulated. The most famous one is called Hebbian learning rule [71]. The learning rule formulated by D. Hebb in 1949 as "When an axon of cell A is near enough to excite a cell B and repeatedly or persistently takes part in firing it, some growth process or metabolic change takes place in one or both cells such that A's efficiency, as one of the cells firing B, is increased". This theory provided one of the most famous neural axiom of neuroscience history : neurons that "fire together wire together" [START_REF] Shatz | The Developing Brain[END_REF]. Various explanation models emerged to account for this learning and adaptation phenomenon such as the covariance rule [START_REF] Sejnowski | Storing covariance with nonlinearly interacting neurons[END_REF], the Oja's rule [START_REF] Oja | Simplified neuron model as a principal component analyzer[END_REF] or the Bienenstock-Cooper-Munro rule [START_REF] Bienenstock | Theory for the development of neuron selectivity: orientation specificity and binocular interaction in visual cortex[END_REF] to explain how synaptic coupling or synaptic transmission efficacy is driven by correlations in the firing activity of pre-and post-synaptic neurons. All these frameworks were a breakthrough in neuroscience since they theorized the fundamental mechanisms underlying the emergence and reinforcement of neuronal networks or ensembles. Such a plasticity, at the mesoscopic scale, might be the origin of a stabilizing of neuronal activity patterns in the brain [71]. If these stabilizing patterns correspond to cognitive representation of behavior, it means that the non-supervised adaptation of synaptic weights could be the process to learn the behaviors and create a mental representation from sensory inputs for example. Hebb's postulate embeds two fundamental notions: locality and joint activity. It means that the change of synaptic efficacy can only depend on local synaptic variables. Joint activity means that preand post-synaptic neurons need to be simultaneously active for a synaptic weight change to occur. The recurrent simultaneous coactivity of several neurons of a neuronal population leads to the framework of population coding within a new neural doctrine: neural information is encoded through neuronal ensembles.

Hebbian rule explains how stable groups of interconnected neurons, known as neuronal ensembles, emerge during activity and neural network development.

The connectivity of neural networks

Neural networks are non-homogeneous but hierarchical, subdivided into specialized computation modules [START_REF] Bullmore | Complex brain networks: graph theoretical analysis of structural and functional systems[END_REF] [77] [78] [79] [80] [START_REF] Rubinov | Neurobiologically Realistic Determinants of Self-Organized Criticality in Networks of Spiking Neurons[END_REF]. The heterogeneity of micro structural properties implies that responses to stimuli vary largely depending on processing, behavior, or cognitive operations. Neural networks are complex [START_REF] Tononi | A measure for brain complexity: relating functional segregation and integration in the nervous system[END_REF], able to perform either functional segregation that enables the rapid extraction of information and the generation of coherent neural states, and functional integration that integrates information through interplay, interaction, and exchange between multiple networks that can encode several tasks or cognitive function simultaneously. Thus, no clear bijection is assumed between a (sub-)module and a cognitive, sensory, or motor process. Many physical and experimental models have been developed to explain how topology of neural networks emerge from single neuron connectivity (see Appendix A). Understanding cognitive functions, thus, implies to understand the exact role of the in vivo networks and the way they shape a mental representation. The network connectivity also known as connectome supports this neural code and must be deciphered [START_REF] Friston | Beyond Phrenology: What Can Neuroimaging Tell Us About Distributed Circuitry?[END_REF] [START_REF] Eickhoff | Functional Connectivity[END_REF].

The connectome is based on three different kinds of connectivity: anatomical (structural), functional and effective [START_REF] Friston | Functional and effective connectivity in neuroimaging: A synthesis[END_REF] [START_REF] Sporns | Classes of network connectivity and dynamics[END_REF]. The anatomical connectivity corresponds to hard wired connections between neural networks. It is based on physical links interconnecting neurons that convey information at a biological level through axons, dendrites, chemical or electrical synapses... The functional connectivity corresponds to all interactions between networks that embed coactivity dynamics not directly related to physical connections. If anatomical connections correspond to the hardware part of the computer, functional connectivity is the software that will process the information making neuron populations interdependent and implement cognitive functions. The effective connectivity summarizes a more complex influence of one node to another. Its main difference with functional connectivity is the integration of the direction flow in the form of excitatory/inhibitory effects, causality (who triggers who) and evolving time dynamic. The latter does not necessarily imply a direct anatomical link between agents and may involve indirect connections, feed-forward backward loops, or cascade events.

The key point of analyzing neuronal functional connectivity is to find specific role for neurons, estimate groups of neurons whose activities related to tasks might unveil the emergent properties of neuron firing synchronization patterns embedded in a modular and hierarchical architecture. Understanding how computation modules appear, interact, or evolve in relation to the various types of connectivity could provide new insights of the way information is embedded in the neural code and generalize at a higher scale the discovered principles: the brain [START_REF] Aertsen | Dynamics of functional coupling in the cerebral cortex: an attempt at a model-based interpretation[END_REF] [START_REF] Gerstein | Simultaneously Recorded Trains of Action Potentials: Analysis and Functional Interpretation[END_REF]. (see Appendix A, B)

Population coding

In population coding, groups of neurons appear through physiological co-active events, corresponding to neurons firing within a same meaningful time window. The coactivation of neurons strengthens their functional coupling regarding the Hebbian rule, making their joint activity more likely to occur. Population coding theory claims that the incident activation of such interconnected networks of neurons might encode a mental representation of a cognitive task. Co-active populations would appear as fundamental units of computation and the substrate of sensory, behavioral, or cognitive functions [START_REF] Carrillo-Reid | What Is a Neuronal Ensemble?[END_REF]. The communication between neuronal populations through coincident biophysical events is supported by an optimized distributed system whose topological structure matches the neuron groups (see Appendix A). The topological structure confers to these neuron groups the intrinsic ability to trigger a inter-group self-maintained response. This response induces basal sequential patterns whose existence is proven but exact role is not yet fully understood but could be related to learning through neuron group stabilization.

Neuron groups can be artificially recalled [START_REF] Carrillo-Reid | Encoding Network States by Striatal Cell Assemblies[END_REF] by optogenetics. These groups of neurons can be activated by triggering specific neurons with particular roles, called Completion Pattern Neurons, that seems to be the cores of the synchronicity inside the interconnected groups. The synchronization patterns and consequently neural network topology can be altered by neurological diseases that lead to neural network functioning abnormalities. For example, some evidence seems to describe schizophrenia in mice as a dysregulation of neural synchrony [START_REF] Hamm | Altered Cortical Ensembles in Mouse Models of Schizophrenia[END_REF]. Indeed, it seems that a disorganization of activation dynamics appears during the presentation of a visual stimulation. Locally, the persistent destabilization of synapses and dendritic spines seems to induce a deprogramming of the groups that are the supports of a stable cognition. Conversely, the ablation of completion pattern neurons could equally alter the ensemble activation and deprogrammed behaviors. In Hydra, for instance, the somersaulting behavior activated by the neural microcircuit called Rhythmic Potential 1 can be reduced by an ablation of some core neurons of this group [9].

The existence of neuron groups and synchronization patterns in neural networks lead to a new neural doctrine based on neuronal ensembles [START_REF] Carrillo-Reid | What Is a Neuronal Ensemble?[END_REF].

From neural networks to neuronal ensembles i A new neural doctrine based on neuronal ensembles

The definition of a neuronal ensemble or neural assembly is intrinsically elusive [START_REF] Carrillo-Reid | Cell Assembly Signatures Defined by Short-Term Synaptic Plasticity in Cortical Networks[END_REF] [93] [START_REF] Diana | Bayesian inference of neuronal assemblies[END_REF] in the literature. It mainly depends on experimental data, mathematical models, but above all the prism of our understanding of the neural computation. A commonly adopted manner to define such neuron ensemble is as a population of nervous system cells involved in a particular cognitive computation [START_REF] Papadimitriou | Brain computation by assemblies of neurons[END_REF]. This definition is even more complicated to formulate since neural network activity largely depends on the maturation of the nervous system, the targeted cortical area, or the model animal under the scope. For example, in early development the response of nondifferentiated cells and immature modules do not demonstrate any group activity (A.4), animals can sometimes react to a stimulus without a neuron group activation, or neural circuitry often process information by a large distribution of the computation over multiple entities leading to an overlap of ensemble activity mediating similar cognition. In such systems, a formal definition of a neuronal ensemble becomes even more complicated.

Originally the first neural doctrine from S. Ramón y Cajal (1899) [START_REF]La Textura del Sistema Nerviosa del Hombre y los Vertebrados[END_REF] [START_REF]Estudios sobre la cortexa cerebral humana[END_REF] and C. Sherrington (1906) [START_REF] Sherrington | Observations on the scratch-reflex in the spinal dog[END_REF] was claiming that individual neurons were the structural and functional units of the nervous system. Later studies have nevertheless demonstrated that perceptual representations of the external world and cognitive computation like working memory, decision-making, consciousness, are encoded by neuronal populations [START_REF] Brown | On the nature of the fundamental activity of the nervous centres; together with an analysis of the conditioning of rhythmic activity in progression, and a theory of the evolution of function in the nervous system[END_REF] Here, to provide a mathematical framework able to extract the neuron ensembles, we decided to adopt the neutral definition of neuron ensemble or assembly [106] [92] [107] [START_REF] Carrillo-Reid | What Is a Neuronal Ensemble?[END_REF] as a group of neurons repeatedly firing together. It means that the consistency of a group is based on the statistical co-activation of neuron population vectors [START_REF] Carrillo-Reid | Endogenous Sequential Cortical Activity Evoked by Visual Stimuli[END_REF] that repeat over time.

ii Neuronal ensembles embed randomness

Observed co-activation of neurons is inherently stochastic. Indeed, each time a neuronal ensemble is active, the co-active neurons taking part in the activity of this group are not the exact same. It leads to a huge variability in the group response. Furthermore, neurons are characterized by a firing basal activity, i.e. neurons can fire outside the activity of their respective groups. Some information is carried out by neuronal ensembles' activation but equally by their silences. A non-response of a group or a subgroup might represent inhibitive effects, modular and hierarchical responses, or functional overlapping. Finally, the mathematical definition of a neuronal ensemble could be : "a neuronal ensemble is a particular set of neurons whose co-activation is statistically significant compared to chance (for a given basal activity of individual neurons)".

Another major issue is to relate the neuron ensemble definition with specific cognitive tasks or experimental stimulations. In other words, it questions whether the neuron ensemble exists independently of the task and the stimuli or should its existence be defined in correlation with the input process. All these concepts are primordial, to highlight the properties and limitations of the existing algorithms to identify neuronal ensembles.

iii Neuronal ensembles are functionally overlapping

Due to the modular organization of the nervous system, (sub-)modules communicate to distribute the processing of cognitive tasks or mental representations [109] [110] [111] [112] [START_REF] Markram | Interneurons of the neocortical inhibitory system[END_REF]. In parallel, the reusing of a single module to encode different behaviors or mental tasks would represent a non-negligible gain in neural computation efficacy and complexity.

Therefore, the exact relationship between a neuron ensemble and a cognitive function needs to be clearly defined. If a neuron ensemble is defined as a subset of neurons that encode one task, behavior or state, it provides a functional definition of an ensemble. It means that an ensemble would be in bijection with a function. Functional overlapping, thus, corresponds to the ability of a neuron to be recruited or recalled by several functional ensembles. These intersecting neurons take part in several cognitive tasks such as neurons of the visual cortex activated in response to both specific stimulus shape and motion. These neurons could have major roles in the neural computation.

An overlapping response of neurons is expected, for a lot of in vivo neural networks, especially in mice and zebrafish visual cortex [START_REF] Smith | Parallel processing of visual space by neighboring neurons in mouse visual cortex[END_REF] [3]. For mice, Allen Institute provides an open-source dataset of the visual cortex, called Brain Atlas, whose neuron responses to different visual stimuli (Static, Drifting Gratings, Natural Scenes and movies, Locally Sparse Noise) directly show strong overlapping responses. In zebrafish, J. Mölter, L. Avitan and G. Goodhill [3] highlighted overlapping responses equally. In the Hydra vulgaris case, neuron ensembles are expected to be non-overlapping [2] structurally and functionally. Physiologically, neurons of different ensembles are non-directly connected to each other, functionally they just belong selectively to at most one single neural network. Moreover, neuron ensembles are correlated with behaviors of the animal: longitudinal contractions, elongation, radial contraction, and nodding.

Neuronal ensembles are, thus inherently characterized by activation randomness and overlapping, which combined with the high levels of noise in free-moving animals, hinder the robust and accurate detection of neuronal ensembles from the in vivo imaging of neuronal populations. Standard techniques for the identification of neuronal ensembles are based on theories and experimental observations about dynamics and topology emerging during development through the functional interactions between neurons arranged in interconnected networks. These physics theories of networks topology and structure led to detection algorithms searching for particular patterns like assortative modules. Providing tools able to extract neuronal ensembles from microscopy imagery and correlate these ensembles with stimulations or spontaneous activity would pave the way to a better understanding of modular neural coding.

Chapter III

Chasing neuronal ensembles that encode behaviors 1 Imaging and tracking single neurons in behaving animals i Neural activity monitoring : From patch clamp to calcium and voltage indicators

Monitoring neuron activity is required to understand how neural networks embed information to integrate input stimuli or trigger behaviors. The early techniques using patch clamp used to be the gold-standard techniques to record electric neural activity of action potentials. However, these techniques were invasive and limited to a restricted number of neurons with poor spatial-temporal resolution.

The advancements of fluorescent calcium microscopy in the 1990s and the progresses of genetically engineering in 2000s allowed to perform "multi-cell imaging". It addresses the activity monitoring of individual neurons at a higher scale by recording proxy indicators instead of direct action potentials. The gold-standard method was the protein-based Genetically Encoded Calcium Indicator (GECI), GCamP [START_REF] Nakai | A high signal-to-noise Ca2+ probe composed of a single green fluorescent protein[END_REF]. Each time an action potential is emitted, a calcium influx occurs at axon terminals by the opening of voltage-gated Ca 2+ (see The method requires an accurate quantification of the light emitted by fluorescent indicators received on the camera's detector. Several fluorescent calcium microscopy techniques have been invented and applied on a huge diversity of model animals such as the spinning disk confocal microscopy, the structured illumination microscopy, or the 2-photon microscopy.

Over the last few years, more direct approaches to record membrane potential changes called Genetically Encoded Voltage Indicator (GEVI) have been introduced. It uses protein directly sensitive to membrane potentials, whose molecular conformation and optical properties can change relatively to voltage. GECI offers several advantages for multi-cell imaging compared to GEVI [START_REF] Kim | Fluorescence imaging of large-scale neural ensemble dynamics[END_REF] [START_REF] Zhu | Identifying the pulsed neuron networks' structures by a nonlinear Granger causality method[END_REF]. Its SNR is 8 to 20 times higher than GEVI's, providing a better spatial resolution, since fluorophore molecules are precisely located in neuron cell bodies while potential voltage can contaminate an entire neuropil corresponding to dendrites, axons, synapses, or glial cells located in the cell neighborhood. Furthermore, they are characterized by a better sensitivity and are more stable to photobleaching effect. Finally, they are more versatile and straightforward in terms of use and applications. However, a large temporal discrepancy between calcium and voltage dynamics has been observed [START_REF] Zhu | Population imaging discrepancies between a genetically-encoded calcium indicator (GECI) versus a geneticallyencoded voltage indicator (GEVI)[END_REF]. GECI misses sub-threshold depolarization, is polluted by external cellular processes that produce noticeable changes in calcium ion concentration and suffers from poor temporal resolution that exaggerates temporal summation and mixes bunch of close spikes. GEVI is more robust than GECI to fluorescent baseline and capture region to region propagation latencies.

GECI and GEVI are both fluorescent imagery techniques that suffer from limitations in freely behaving animals. Firstly, their accuracy suffers from noise due to motion artifacts [START_REF] Hamel | Cellular Level Brain Imaging in Behaving Mammals: An Engineering Approach[END_REF] [START_REF] Kerr | Spatial Organization of Neuronal Population Responses in Layer 2/3 of Rat Barrel Cortex[END_REF]. Indeed, motion of the animal itself, heart rate, breathing, non-specific vibrations can blur the image recordings. In addition, calcium microscopy techniques excite fluorescent molecules and must account for back-scattering effects of light in the tissues, and photobleaching effect corresponding to photo-toxicity that leads to the transition to permanent dark state of fluorescent cells over time and results in a non-linear decrease of the fluorescent intensity. For animals without spatial movement constraints, the displacements, the z-axis motions, or deformations can also lead to time-varying fluorescent baseline. Finally, the superimposition of out-of-the-focus plane signal imaging and the activation of surrounding agents in the neuropil can lead to distorted fluorescence intensities.

Hardware and software corrections have been implemented and are described hereafter.

ii Imaging neurons

Hardware solutions correspond to the implementation of experimental protocols and microscopy techniques to limit the previous effects [START_REF] Sawinski | Visually evoked activity in cortical cells imaged in freely moving animals[END_REF] when monitoring the neural activity. To limit motion artifacts, one solution is to physically constrain the motion of the animal.

Mechanical devices to immobilize the animal are used, for example, head fixation under the microscope field-of-view [START_REF] Karagyozov | Recording Neural Activity in Unrestrained Animals with Three-Dimensional Tracking Two-Photon Microscopy[END_REF] [124] [START_REF] Seelig | Two-photon calcium imaging from head-fixed Drosophila during optomotor walking behavior[END_REF] [4] [START_REF] Carrillo-Reid | Identification and Targeting of Cortical Ensembles[END_REF], paralysis via anesthesia [5], surgical surgery or constraining substrates [START_REF] Ahrens | Large-scale imaging in small brains[END_REF] [128] [3].

Then, the basic operating principle of the microscopy methods is to reduce the tissue scattering effects and photobleaching improving spatial and temporal resolution. Two main microscopy techniques are used: the spinning disk confocal microscopy (SDCM) and the two photon-microscopy (2P-M). SDCM (see Figures III.3, III.4), invented by Petráň et al, in 1968 [START_REF] Petráň | Tandem-Scanning Reflected-Light Microscope*[END_REF], scan a sample by rotating arrays of aperture pinholes corresponding to microlenses that act like excitation focal volumes. As output, it provides a high-rate real-time confocal image of the in-focus sample. Noticeable improvements of emission photon collection, detection sensitivity, frame rate to equivalent SNR, photobleaching toxicity reduction were proved using this method compared to previous classical microscopy techniques [START_REF] Oreopoulos | Spinning-disk confocal microscopy[END_REF]. [START_REF] Denk | Two-Photon Laser Scanning Fluorescence Microscopy[END_REF], in 1990. This method allows to image deeper in the tissues by making several photons interfere with each other. While in confocal microscopy, the light selectivity is performed through optical sectioning with pinhole arrays, this method uses longer laser wavelength and the optical theory of up-conversion to do so. The molecule excitation is a non-linear process that involves the absorption of two infrared photons whose combined energy is greater than the energy gap between molecule's ground and excited states allowing a molecular transition. This optical phenomenon allows a better molecule targeting, the ability to image at different depth in the tissue corresponding to different area of the visual cortex for instance and to reduce tissue photo-toxicity.

It is worth noting that other imagery techniques such as light-sheet microscopy allow volumetric imaging with good time resolution (∼ 1 Hz) of small brain volumes or whole animal models [START_REF] Hillman | Light-Sheet Microscopy in Neuroscience[END_REF].

iii Tracking calcium activity Software solutions, inspired from signal processing, aim at extracting single neuron calcium activity from microscopy imaging. These techniques try to account for 2 main issues : a) Either the animal is fixed by hardware solutions and the motion is just residual, or the animal is freely behaving, and more advanced tracking techniques are required. In the first case, residual motion are usually vibrations, obtained for example using 2P-M in mice brain visual cortex [START_REF] Giovannucci | CaImAn an open source tool for scalable calcium imaging data analysis[END_REF]. The frame averaging and the correction of deviations from equilibrium positions can be sufficient to account for such motion artifacts. In the second case, when the animal is displacing and deforming freely, like for Hydra Vulgaris [START_REF] Lagache | Tracking calcium dynamics from individual neurons in behaving animals[END_REF], the robust estimation of the elastic deformations of a free-migrating organism has been previously studied in Caenorhabditis Elegans [START_REF] Nguyen | Whole-brain calcium imaging with cellular resolution in freely behaving Caenorhabditis elegans[END_REF] [136] [START_REF] Ahrens | Large-scale imaging in small brains[END_REF] or in Hydra [START_REF] Chenouard | Multiple Hypothesis Tracking for Cluttered Biological Image Sequences[END_REF] [135] imaged using Spinning Disk Confocal Microscope.

Therefore, the image analysis pipeline is divided into several sequential signal processing steps as follows : 1) Motion estimation and correction , 2) Neuron identification and tracking 3) Source denoising and Spike deconvolution, 4) Neuronal ensembles identification [START_REF] Cantu | EZcalcium: Open-Source Toolbox for Analysis of Calcium Imaging Data[END_REF] [134] [START_REF] Pachitariu | Suite2p: beyond 10,000 neurons with standard two-photon microscopy[END_REF].

a Motion estimation and correction

The robust estimation of motion in the focal plane of a freely-behaving animal is fundamental to target a neuron on an entire video, correct its displacement and efficiently extract its activity in terms of fluorescent intensity. For fixed animals with only residual motion artifacts, classical elastic registration methods are used. For example, the NoRMCorre tech-nique [START_REF] Pnevmatikakis | NoRMCorre: An online algorithm for piecewise rigid motion correction of calcium imaging data[END_REF] based on template matching is a well established semi-rigid corrective method. The idea is to find piece-wise rigid displacements by comparing every frame of the recording to a template averaged image. Each frame is divided into up-sampled overlapping patches that match subareas of the template to estimate the local cross-correlation maxima and thus the local displacement vectors. Then, the template is updated by averaging the frames corrected by interpolated displacement vectors. Another widely used registration algorithm is part of the CaImAn project [START_REF] Giovannucci | CaImAn an open source tool for scalable calcium imaging data analysis[END_REF], but it worth noting that a huge diversity of other techniques exists : intensity based methods [START_REF] Thevenaz | A pyramid approach to subpixel registration based on intensity[END_REF] [142] [START_REF] Pachitariu | Suite2p: beyond 10,000 neurons with standard two-photon microscopy[END_REF], inter-frame matching key points based methods [START_REF] Liu | Fast and Accurate Motion Correction for Two-Photon Ca2+ Imaging in Behaving Mice[END_REF] or deep-learning techniques [START_REF] Haskins | Deep learning in medical image registration: a survey[END_REF] [START_REF] Thompson | The Computational Limits of Deep Learning[END_REF] are part of them.

For freely behaving and deforming animals likeHydra Vulgaris this issue is even more complicated and still appears as a partially open-problem. Indeed, the animal's motion combined with the elastic deformations of the body [START_REF] Lagache | Tracking calcium dynamics from individual neurons in behaving animals[END_REF] make the displacement impossible to perfectly correct. However, some promising efforts have been performed to efficiently track the moving neurons and use the tracks as indicators of the animal deformation.

Following motion correction, an segmentation is performed to extract the fluorescent intensity a neuron or soma, which corresponds to a contiguous set of pixels constituting a Region Of Interest (ROI) within the field of view that are correlated in time [START_REF] Giovannucci | CaImAn an open source tool for scalable calcium imaging data analysis[END_REF] [139] and can be active or silent regarding the neural state. At first glance, the mouse visual cortex appears as a quite simple case since, when the motion correction has been properly carried out, the summation or averaging of all frames already provides a good estimate of each neuron location (see Figure III.7). In this case, ROIs are fixed between frames, matched with a neuron and can be simply segmented (see Figure III.8). However, when the density of ROIs is high the process detection needs to be refined. Indeed, the simple averaging of frames, is no longer precise enough since it misses sparse active cells with low baseline. Goldstandard techniques for neuron ROIs demixing are: a) matrix decomposition with positivity constraints embedded in a convex optimization framework [START_REF] Pnevmatikakis | Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data[END_REF] or, b) matrix factorization with greedily segmentation of high-correlated neighboring pixels [START_REF] Pachitariu | Suite2p: beyond 10,000 neurons with standard two-photon microscopy[END_REF]. Their aim is to refine ROI detection and demix calcium sources.

• Constrained nonnegative matrix factorization (CNMF) [START_REF] Pnevmatikakis | Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data[END_REF] CNMF limits overlapping in ROI detection to avoid "crosstalk" that could impair the monitoring of individual neuron activities. This technique embeds a spatio-temporal model of component location and demixing. The spatial component is handled by the introduction of a non-negative and sparse "spatial footprint" that encodes the spatial calcium concentration profile at each timestep. Concomitantly, a temporal calcium dynamic as an autoregressive model summarizes the evolution of calcium in each pixel/voxel. The method alternatively solves several convex optimization problems with sparse and positive constraints to estimate these spatio-temporal features. Specific initialization procedures, ranking, removing, and merging components are performed. However, the algorithm finally highlights fixed location of ROI corresponding to somas compensating overlapping and handling time-varying baseline.

• Suite2p

Suite2p embeds a spatially localized signal model that handles the activity of the ROI, the neuropil contamination, and the noise. Three alternating steps solve iteratively the optimization problem. 1) Data reduction through support vector decomposition applied on each pixel time series. 2) New source detection through pixel-to-pixel trace correlation in neighborhood. Extraction of the maxima as new candidates. 3) Source and neuropil activity updating by square error minimization. 4) Updating of spatial source locations. Repeat steps 2), 3), 4). Here, an equivalence between ROI and neuronal cells spatially fixed in the image is assumed.

In Hydra Vulgaris' processing pipeline, ROIs appear as bright spots detected over each frame. Detection can be performed with standard algorithms such as wavelet decomposition [START_REF] Olivo-Marin | Extraction of spots in biological images using multiscale products[END_REF] and thresholding or deep learning methods (e.g. Stardist [START_REF] Weigert | Star-convex Polyhedra for 3D Object Detection and Segmentation in Microscopy[END_REF] [START_REF] Schmidt | Cell Detection with Star-Convex Polygons[END_REF]) Since animal motion and deformation induces moving ROIs between images, the association in a succession of ROI detections to estimate the activity of a specific neuron is thus required, especially, for neurons with long time-lapse silent states (no firing).

b Neuron identification and tracking

ROIs detected on consecutive frames can be associated in partial tracks called tracklets using standard tracking algorithms such as the probabilistic eMHT [START_REF] Chenouard | Multiple Hypothesis Tracking for Cluttered Biological Image Sequences[END_REF]. However, the prolonged undetectability of non-firing neurons (without calcium signal) requires a tracklet to tracklet stitching to merge the tracklets corresponding to the same neuron. Thus, a first promising technique is EMC2 [START_REF] Lagache | Tracking calcium dynamics from individual neurons in behaving animals[END_REF] (see Figure III.9). It associates the most relevant tracklets to reconstruct the entire tracks. To do so, the algorithm estimates the elastic deformation of the animal by interpolating via a thin-plate-spline function the position of non-detected neuron ROIs due to long silent neural states from the detected ones. Finally, a global tracklet-to-tracklet distance is optimized in the elastic corrected space to reconstruct the tracks. The elastic motion of the animal is extracted on the fly. Some recent efforts from [START_REF] Reme | Tracking Intermittent Particles with Self-Learned Visual Features[END_REF] The diversity of spike inference techniques is important and can be divided into several working principles [START_REF] Diana | High frequency spike inference with particle Gibbs sampling[END_REF] From the vast number of spike inference techniques three main categories of spike inference techniques emerge: 1) Deterministic techniques using optimization, 2) statistical techniques using inference based on a model and 3) machine learning techniques. Generally, the methods formulate analytically a model that describes the dependency between fluorescence intensity, calcium concentration, and spiking activity. The model characterizes the relationship between the three time-series by embedding the fluorescent noise. To fulfill that task, a transient kernel that summarizes the calcium dynamics is required and embedded in an auto-regressive process whose parameters are assumed to be fixed. A deterministic or stochastic optimization algorithm is, then, run on the individual neuron fluorescence traces to infer jointly the most relevant spiking sequences, and calcium time-series. The level of noise, the kernel dynamics and more rarely the baseline are estimated directly.

Such techniques, usually assume the flatness of the fluorescence traces and spikes distributed according to a homogeneous Poisson distribution [START_REF] Pnevmatikakis | Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data[END_REF] [173] [START_REF] Vogelstein | Fast Nonnegative Deconvolution for Spike Train Inference From Population Calcium Imaging[END_REF], to provide fast, accurate and robust estimates of the action potential emission dates that forms the neuron spiking matrix. Some efforts [START_REF] Theis | Benchmarking Spike Rate Inference in Population Calcium Imaging[END_REF] [153] [START_REF] Pachitariu | Robustness of Spike Deconvolution for Neuronal Calcium Imaging[END_REF] have been carried out to benchmark few methods on synthetic data but remain insufficient to identify a gold-standard technique and assert the efficiency of the techniques on real neuronal experimental recordings.

The neuronal spiking matrix uncovers the convergence of individual neural activity through the occurrence of co-active biophysical events corresponding to action potential emissions. These are supposed to highlight the activity of neuronal ensembles that are the substrate of cognitive states or mental representations [START_REF] Carrillo-Reid | What Is a Neuronal Ensemble?[END_REF]. The robust extraction of neuronal ensembles from the neuron spiking matrix is supposed to provide insights about information integration.

iv Detection of neuronal ensembles

Community detection techniques are generally theorized in more general contexts (social networks, economy, politics...) than in vivo neural networks. They are based on topological observations, on the detection of recurrent activation patterns or a more exhaustive modeling of the studied phenomenon [192] [193] [194].

The existing methods to detect communities are extremely numerous and generally divided into three main categories (see Figure III.15):

• Graph theory techniques. The fundamental idea of node partitioning in com-munities is to identify sub-graphs of densely interconnected nodes called modules, clusters,communities, groups or assemblies. Regardless of their structure, community nodes are assumed to be "assortative", meaning strongly connected to other nodes in their own community and weakly connected to neurons from other ones. This supports a modular organization of a small-world network (see Appendix A), relevant with the observations made in in vivo neural networks allowing autonomous, specialized, or distributed functions optimizing the node wiring, communication, evolution by being robust to failures. It requires a pairwise node adjacency matrix obtained from a Functional Connectivity (FC) metric between pairwise neuron spiking activities. The connectivity matrix is, then, thresholded and a clustering algorithm on the resulting graph is applied. The algorithm choice is fundamental since it operates according to an own definition of a modular community [195] [196] [197] [198] requiring sometimes a projection of the graph rectifying or discarding information. It equally embeds a set of assumptions on the way the graph has been generated by the functional connectivity metric which have a large impact on the community detection.

• Spectral techniques. They correspond to the application of algebraic methods that extract coherent subgroups from activity matrix [199] [200] [201] [108] [126] [START_REF] Peyrache | Replay of rule-learning related neural patterns in the prefrontal cortex during sleep[END_REF] . Each time step is interpreted as a population vector activation on which machine learning techniques are applied to reduce data dimensionality, projected data in some chosen subspaces to find heterogeneous repeated components and to extract common underlying patterns. The task of extracting relevant repeating populations is usually performed using spectral techniques like principal component analysis, independent component analysis or some improved versions.

• Biologically model based techniques. The idea is to explicitly and analytically model what biologically a neuronal ensemble is and try to uncover neuronal assemblies based on an explicit interpretable definition [START_REF] Diana | Bayesian inference of neuronal assemblies[END_REF]. The main interest of the techniques is to provide meaningful results that handle the specificity of the underlying problem providing statistical evidence. [START_REF] Jebabli | Overlapping Community Detection Versus Ground-Truth in AMAZON Co-Purchasing Network[END_REF] on open-source datasets using the commonly used evaluation metrics : 1) (O)-NMI [START_REF] Danon | Community Structure Identification[END_REF] [212] and, 2) Omega Index [START_REF] Collins | Omega: A General Formulation of the Rand Index of Cluster Recovery Suitable for Non-disjoint Solutions[END_REF]. A gold-standard technique does not exist and generally good estimates on experimental data are obtained by combining several community detection techniques and "aggregating" the clustering results [START_REF] Fortunato | Community detection in networks: A user guide[END_REF]. Some efforts have also been carried out to face promising methods on in vivo neural network for the specific problem of neuronal ensemble detection [3]. The techniques have been benchmarked, firstly, on synthetic data and then, on zebrafish experimental data.

The existing community detection techniques have been benchmarked [203] [204] [205] [206] [207] [208] [209]

Finally, the previous benchmarkings of community detection techniques reveal their efficiency to recover communities on simulated data. However, generalizing results on real-world networks especially for the neural networks in vivo case is hard in the absence of groundtruth. Indeed, the exact ability of the community detection algorithms to account for the specific issues of neuronal ensembles such as stochastic activation or ensemble overlapping has been poorly studied yet. The development of methods able to estimate additional evaluation features could also be relevant. For example, estimating neuronal ensemble activity in relation with behaviors or stimulations in free-living animals is a promising goal. Problem inversion using the model

Model of neuron ensemble interactions

s(t) t s(t) t s(t) t Θ 1 Θ 2 Θ 3

Detection neural assemblies

Neuronal ensembles encode behaviors

Understanding how neuronal ensembles encode behaviors requires to quantify precisely the link between motion of a free-behaving animal and its brain computations. The improvements in techniques to track animal motion and to statistically couple stimuli or spontaneous behaviors to neural circuitry activity take advantage of deep learning techniques. Understanding the extent in which neuron ensembles better explain behaviors than individual neuron cells would reinforce the neural doctrine that ensemble is the cornerstone of neural code. Questioning how to precisely define a behavior, quantify it, correlate it with a neural activity or measure the ability to predict such actions are recent open questions essential to pursue the neural code deciphering.

i Behavior tracking

Detecting the occurrence of behaviors in a free-moving animal requires either to perform a manual labeling or to algorithmically detect the behaviors. Motion of the animal needs, thus, to be quantified in order to identify an exhaustive list of all possible behaviors called the behavioral repertoire of the animal. The quantification of motion often implies to fit control points on the animal to extract position and kinetics properties in relation with behaviors.

The robust tracking of the centroid of the animal is generally the most straightforward way to quantify its spatial navigation and locomotion, complexifying the tracking shape can allow to uncover richer information about behavior [START_REF] Pereira | Fast animal pose estimation using deep neural networks[END_REF] [214]. To identify complex behaviors, pose estimation is usually required (see Figure ii Behavior prediction Linking neural activity with behaviors is at stake to decipher the neural code [START_REF] Pereira | Quantifying behavior to understand the brain[END_REF]. The correlation between neural activity and stimulations or spontaneous behaviors would represent and intrinsic ability to encode, and then to predict (see Figure III.17). Long-Short Term Memory (LSTM) is a deep-learning technique based on recursive neural network (RNN) that can be use to decode information about an animal behavior [START_REF] Tampuu | Efficient neural decoding of selflocation with a deep recurrent network[END_REF]. Animal spatial location, locomotion or deformation can be estimated by inputting neural activity in such artificial neural networks. Its accuracy level makes this tool more accurate than classical Bayesian inference to infer behavior and promising to assess the ability of the input to predict the output. Such tool could be used to test if neural ensembles are more relevant to predict behaviors than individual cells. It has been successfully used to decode animal's position from neuron population activity when applied to hippocampal neural data of free moving rats [START_REF] O'keefe | Hippocampal place units in the freely moving rat: Why they fire where they fire[END_REF].

Deciphering the neural code requires to unveil the intrinsic mechanisms of how interconnected neurons convey information and encode animal state and behavior. This understanding usually requires to go from the in vivo imaging of neuronal populations in behaving animals to the extraction of the most relevant neuronal substrates. This is a multistep complex process that requires the continuous development or improvement of image analysis and signal processing techniques. The main contributions of this thesis, described hereafter, are dedicated to the development of a robust statistical pipeline to decipher the neural code of model animals based on the in vivo imaging of their neuronal activity. The behavioral model is defined by 3 states : elongation, contraction, and rest. Any kind of classifier could be used to make the distinguish between the states. For Hydra, a naive hard-threshold is applied on the distance between control points. c) A predictive algorithm is implemented to predict the behavior. For Hydra, a LSTM is used. The provided input is the average pixel intensity of parts of the image. Variations of the input parts are performed to estimate the level of predictability of the pixel intensity of the different parts of the animal. The idea is to unveil if a specific part is involved in a specific behavior prediction.

Chapter IV

Publications 1 Extracting neuronal ensembles in calcium imaging i Pipeline reminder

We remind the reader that the contributions are structured as follows around the different steps of sequential processing pipeline (see The simulator aims to assess the robustness and efficiency of the methods applied to the Hydra vulgaris case, and to point out potential limitations.

2) Step 2 & 3 -Inferring spikes from fluorescence traces & The estimation of a network topology:

The second contribution is the development of a statistical tool to infer neuron-to-neuron statistical coupling in order to extract the network topology while being robust to false positives. Its use as a functional connectivity method facilitates spike inference by filtering out neuron activation artifacts.

3) Step 2 -Inferring spikes from fluorescence traces:

The third contribution of this work is the application of spatial statistics framework to filter out spike artifacts in the presence of a decorrelated noise. 4) Directly from step 2 to step 4 -Neuronal ensemble detection: The fourth contribution aims to extend existing tools for estimating functional neuronal ensembles by using a Bayesian Inference framework to analyze emergent synchronization properties at network scale from the activation of individual cells. This framework uses synchronicity as a biological interpretable neuronal ensemble definition handling overlapping and embeds the point process statistical coupling as an evaluation of assemblies relevance in response to stimulations. The spatial statistics framework is adapted to measure the correlation between stimulations and neuronal ensemble activity.

ii Difficulties of going from single cell activity to neuronal ensemble estimation

In this thesis, we focus on the four first steps of the decoding process. We detail hereafter the main technical issues encountered in these steps and the statistical tools that we developed to improve the overall accuracy and robustness of the process.

a Tracking individual neurons and extracting calcium fluorescence

Tracking errors

Tracking individual neurons in free-behaving animals are usually two-steps processes with, first, a detection of the neurons on successive frames, followed by the association of detections in coherent tracks using for example, a Kalman filter. This filter corresponds to a model of displacements of the regions of interest based on the estimation of their respective speed given previous frames and anticipating their following positions. When, an animal's motion occurs too fastly regarding the time resolution, the constant speed hypothesis of the model is violated and two different somas can be associated wrongly with the same single track providing a concatenation of two non-relevant dynamics.

This problem is particularly important in Hydra vulgaris that can deform very rapidly and importantly during contraction phases (see Figures III.11, III.12), providing the concatenation of inconsistent dynamics (see Figure V.2).

Out-of-focus motion and cross-talk contamination

Animal's fast motions and deformations can blur the image recording leading to spurious intensity variations. The animal's z-axis motions can equally generate artifact intensity variation that needs to be estimated and corrected.

The high density of somas on images can lead to overlapping and a mixing-up of their activity in a phenomenon called "cross-talks". This phenomenon is observed at mice visual cortex and equally in Hydra vulgaris that gathers during contraction phases lots of somas in a cluttered environment [START_REF] Lagache | Tracking calcium dynamics from individual neurons in behaving animals[END_REF] providing once again a mixing-up between inconsistent dynamics (see Figure V.3). This contamination is even more important that contraction of the animal increases the global fluorescence intensity in its whole body by diffusive effect generating common artifact dynamics in all somas [2] (see Figure V.4). Activity contamination can also be caused by neighboring agents. Indeed, supporting cells allow the ion channels opening in living animals or the neural information to spread in the neural networks. Thus, the recorded calcium response can be partly contaminated by other biological local agents than somas. Especially, glial cells or nematocytes can be spots where intensity calcium response is high or time-varying and if recorded can contaminate the soma's fluorescence trace and mislead in the extracted dynamics.

b Inferring spikes from fluorescence traces

Non-linear baseline fluorescence

Neuropil contamination, uncorrected motion and deformation of the animal, or superimposition of neurons spatially close on the frames pollute the neuron spiking activity with false action potentials. Even if some preprocessing steps can be implemented such as detrending or smoothing, false positive artifacts remain and advanced motion correction, source demixing or global contamination correction are required.

Concretely, spike inference techniques are generally based on the deconvolution of flat stationary signals with spikes distributed according to an Homogeneous Poisson distribution whose calcium dynamics is assumed to be constant. Non-corrected artifacts generate timevarying baseline at different frequencies. The underlying baseline induces non-flat signals whose non-corrected dynamics are assumed to represent meaningful information by the spike inference algorithms. Each time a dynamic occurs, a bunch of spikes is located at that specific moment. The main risk is to match underlying dynamics, corresponding to noise or motion, instead of correlating pairwise neuronal activities.

Non-homogeneous firing patterns

Being robust to time delays is equally at stake in spike trains. Since spike inference techniques generally perform the spike estimation using a fixed estimated deconvolution kernel, often represented by auto-regressive coefficients [START_REF] Pnevmatikakis | Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data[END_REF] [174], the deconvolved sequence can demonstrate artificial lags. The problem is exacerbated in Hydra because the kernel dynamics are expected to be neuron firing regime dependent. A kernel that would not fit properly the time varying calcium dynamics could introduce time delays leading to bias. It is worth noticing that the limitation of non-adaptive deconvolution kernel was highlighted in 2020 by our article [START_REF] Kubler | A Robust and Versatile Framework to Compare Spike Detection Methods in Calcium Imaging of Neuronal Activity[END_REF]. The use of Markov-Switching Autoregressive Model (MS-AR) [START_REF] Ailliot | Markov-switching autoregressive models for wind time series[END_REF] to account for the dynamic changes in different firing regimes was a track to consider. No long afterwards, some efforts of Bayesian Inference techniques have been carried out to handle time-varying baseline and dynamic changes of kernel in the deconvolution [START_REF] Diana | High frequency spike inference with particle Gibbs sampling[END_REF] providing insights that our conclusions were relevant.

c Estimating the functional connectivity between individual neurons

Choosing a functional connectivity metric

It is important to bear in mind that most functional connectivity analyses in Neurosciences focus on large scale ensemble-level neuronal signals often assumed to be stationary Gaussian [START_REF] Sheikhattar | Extracting neuronal functional network dynamics via adaptive Granger causality analysis[END_REF]. It is still unclear how well each metric can be adapted to address single-cell level descriptions, such as population calcium imaging data. However, other communities exist and handle specific problem about spiking pattern through a mathematical field called Point process theory with some promising but non-straightforward tools like Granger Causality (see Appendix C), Generalized Linear Model (GLM) or Renewable and Non-Rewable Hawkes process... This field addresses the statistical coupling of point process through the works of Wilson Truccolo, Emery Brown, Liam Paninski or Jonathan Pillow to only cite few of the main contributors [START_REF] Truccolo | Collective dynamics in human and monkey sensorimotor cortex: predicting single neuron spikes[END_REF] [START_REF] Brown | Donald O. Hebb and the Organization of Behavior: 17 years in the writing[END_REF]. The direct application of their work seems, to our knowledge, limited by the experimentalists since the mathematical theories are generally quite advanced, hard to implement and non-related directly to biological plausible problems. For this reason, their work is just mentioned here but is not under the scope of this manuscript.

[224] [225] [226] [227] [20] [62] [228]
It is, however, worth to noting that, existing toolboxes implement the main functional connectivity metrics: SIFT toolbox [START_REF] Delorme | ERICA: New Tools for Advanced EEG Processing[END_REF], MVGC toolbox [START_REF] Barnett | The MVGC multivariate Granger causality toolbox: A new approach to Granger-causal inference[END_REF], Hermes toolbox [START_REF] Niso | HERMES: Towards an Integrated Toolbox to Characterize Functional and Effective Brain Connectivity[END_REF], TRENTTOOL toolbox [START_REF] Lindner | TRENTOOL: A Matlab open source toolbox to analyse information flow in time series data with transfer entropy[END_REF], Chronux toolbox [START_REF] Bokil | Chronux: A platform for analyzing neural signals[END_REF].

How to choose efficiently a functional connectivity metric is still elusive in a context where data are not stationary Gaussian time-series signals, with Non-Homogeneous Poisson spike distribution and non-flat fluorescence baseline. All these points would need to be considered to correctly identify the information flow but are usually neglected and user-defined without further explanations.

Graph Theory and Spectral methods require thresholding

Once the functional connectivity metric has been chosen, it is applied on each pair of nodes to generate the connectivity matrix for community detection via graph theory. According to the metric chosen different natures of information are expected. If the metric handles causality, positive and negative coefficients in the matrix are possible. If a simple correlation is applied, coefficients are expected to be positive. It can provide different kind of matrices. We previously mentioned that topology structures are expected to be smallworld with topological features (see Appendix A), characterized by specific node degree, clustering coefficient, average path length or the existence of hubs, bridges, rich-clubs or cliques. A distance is equally required for spectral techniques based on the definition of an euclidean space and sub-space whose scalar product definition allows the comparison between neuronal populations.

Generally, the connectivity matrix is binarized to highlight the topological structure and neuronal population thresholded using a statistical z-score test so as to run downstream a Neuron Clustering Algorithm. The thresholding step deletes edges too close to noise corresponding to spurious connections. Different techniques can be applied to obtain this Thresholded Connectivity Matrix.

• A hard threshold : a statistical threshold is applied using a certain number of standard deviations relatively to the mean. Only adapted if the number of connections follow a gaussian distribution. [235] [236] [237] [START_REF] Zalesky | Connectome sensitivity or specificity: which is more important?[END_REF] • Shuffling techniques to generate surrogate data : it shuffles the data to generate time independent spike patterns to create a H0 null hypothesis and remove all links too close to a random situation. [START_REF] Maccione | Multiscale functional connectivity estimation on low-density neuronal cultures recorded by high-density CMOS Micro Electrode Arrays[END_REF] [240] [START_REF] Grün | Analysis of Parallel Spike Trains[END_REF] The average link path is shorter between high functional connected nodes. The userdefined threshold modulates the structure of the functional connectivity matrix. The lower the binarization threshold, the higher the number of connections and the more homogeneous the link distribution. The higher the threshold, the shorter the link path and more visible the patterned networks. That suggests according to [START_REF] Poli | Functional connectivity in in vitro neuronal assemblies[END_REF] "the importance (and the dependence) of the threshold selection in this kind of measures" since it directly alterates the representation of the topology known to be fundamental for community extraction. Thus, the exact influence of the choice of the user-defined threshold combined with a specific functional connectivity metric and their application regarding a specific neuron clustering algorithm has been little studied and provides a huge diversity of results in the estimation neuronal ensembles.

d Extracting neuronal ensembles A huge diversity of methods

An exhaustive listing of all community detection methods that could be applied on neural network is impossible and going into the exact details of each method is out of the scope of this thesis. The most used or promising techniques applied on in vivo neural networks are provided in Appendix B. A figure is provided to summarize a method classification based on our understanding of the literature and to highlight methods, from the community detection, that could be applied on in vivo neural networks. First and foremost, it is worth noticing that the classification taxonomy largely varies according to the reviews and rarely forms a disjoint partition of methods.

The lack of ground-truth

Benchmarking the community detection techniques is extremely hard since each method operates using its own vision and interpretation of how a neuronal ensemble emerges. Benchmarking methods on synthetic data are often biased since better results are reached for estimation techniques that fit the way data have been generated. However, real world experimental data often provide a higher level of complexity than simulated ones making the technique efficiency hard to estimate without ground-truth, and results obtained on synthetic data hardly generalizable.

Mölter and Goodhill [3] provided the most promising benchmarking study of community detection techniques in the context of in vivo neural network activity firstly on surrogate datasets of calcium imaging data and then on real experimental data of zebrafish optic tectum responses to simple visual stimuli. They benchmarked spectral techniques ICA [START_REF] Lopes-Dos Santos | Neuronal Assembly Detection and Cell Membership Specification by Principal Component Analysis[END_REF] [200], PCA (Promax [START_REF] Romano | Spontaneous Neuronal Network Dynamics Reveal Circuit's Functional Adaptations for Behavior[END_REF], similarity analysis (CORE [START_REF] Miller | Visual stimuli recruit intrinsically generated cortical ensembles[END_REF]), graph theory (SGC [START_REF] Avitan | Spontaneous Activity in the Zebrafish Tectum Reorganizes over Development and Is Influenced by Visual Experience[END_REF]) and frequent item set mining (FIM-X) for different simulated parameters and by comparing the results with Reference Assembly Configuration (RAC) on real data using BestMatch Index [START_REF] Goldberg | Measuring Similarity between Sets of Overlapping Clusters[END_REF].

Finally, SGC was, according to the authors, the best method to extract the neural ensembles since it provides a number of communities more relevant with the expected number and a better BestMatch score compared to RAC. RAC are assemblies whose fluorescence response after each stimulation is statistically positively significative. The use of such ensembles to account for neuronal ensemble estimation efficiency is a user-defined choice that reflects the lack of ground-truth for method accuracy evaluation.

The main limitation of the benchmarking is the inability to correlate the activity of the inferred neuronal ensembles with stimulations. Providing assemblies in response to evoked visual stimulations without measuring a coupling to describe how well these ensembles match and explain the stimulations is a lack. The chosen evaluation metric is equally not a standard metric like NMI or Omega Index usually used in the community detection field. Finally, the main limitation of spectral methods is the overlapping case, which is not, to our knowledge, benchmarked on the simulator or further discussed on real dataset, even if overlapping in visual cortex of model animals has been demonstrated. Some theories have also been developed claiming that it is the combination of results of several clustering methods that allows an efficient and robust detection of neural ensembles. Techniques to merge and average the clustering results have been developed known as consensus clustering [START_REF] Fortunato | Community detection in networks: A user guide[END_REF]. What about biological interpretability ?

Many methods do not account for the specificity of the underlying biological problems. It means that the results are equally hardly explainable since they do not provide interpretable estimates supported by statistical evidence. Implementing a biological-model based technique is at stake to fulfill this task allowing the technique to be largely adopted by a community of neuro-experimentalists.

In the next following sections, we present our four main contributions that partially tackle previous issues. These contributions pave the way to a robust and accurate pipeline to go from calcium imaging to the characterization of neuronal ensembles in vivo and their functional implication.

A versatile simulator for spike inference technique that embeds Hydra vulgaris neural activity properties

Calcium fluorescence imaging enables real-time activity monitoring of single neurons in living animals. A critical inverse problem resides in the precise inference of spike locations from noisy fluorescence traces, especially in the presence of burst spiking and non-linear fluorescence intensity. Several spike extraction algorithms have been proposed in the recent years, but a robust and objective evaluation of their performance still remains elusive due to the unsupervised nature of the problem. Here we propose a biologically-inspired mathematical framework to reproduce synthetic fluorescence traces from a time-series data of spike-trains. The idea is to create a versatile platform to objectively test the state-of-the art spike inference methodologies over a large range of experimental parameters. Our solution appears as a complementary but more exhaustive approach to determine the robustness of existing solutions to different nature of signals, imaging artefacts, sensitivity to hyperparameters and pre-processing steps. We benchmark state-of-the-art algorithms with the proposed simulation platform, and validate the results on an experimental dataset of the Hydra Vulgaris. We hypothesize that, in contrast to the common practice of qualitative evaluation, quantitative measure of algorithm robustness is essential in understanding the suitability of a spike inference algorithm to be used in an automated computational pipeline to decipher the neural code.

Introduction

To understand the emergent computational properties of connected single neurons, it is necessary to monitor the coordinated activity of many single cells with high temporal resolution [1]. Fluorescence calcium imaging [2] remains the gold-standard for such studies. Thanks to the recent advances in calcium probe engineering and optical microscopy, thousands of interconnected neurons can now be imaged in living animals with high temporal resolution [3]. The success of reverse engineering the brain now hinges upon effective signal processing techniques to extract critical information from imaging datasets.

In addition to the automatic tracking of single neuron fluorescence [4], another critical issue in the processing of calcium imaging data is that of analysing neuronal spike trains from the videos of firing neurons (see Fig. 1).

The inverse problem of reconstructing the neuronal spikes from the calcium imaging data is inherently unsupervised. Given a temporal sequence of calcium uptake of the firing neurons, it is difficult, even for a neuroscientist, to reconstruct a precise representation of the spiking pattern. State-of-the-art deconvolution methodologies exhibit significant variability in extracting the spike patterns from same videos [5] [6] [7]. Due to the absence of ground-truth annotation in most experimental data, qualitative evaluation appears to be a complementary norm to measure the efficacy of deconvolution algorithms [6]. However, supplemented by quantitative experimental evaluations, such qualitative assessments are unreliable due to inherent subjective bias. In the absence of gold standards, the proper mechanism to study the behaviour of deconvolution algorithms is via rigorous testing on realistic synthetic data.

In this paper we present a robust generative model to realize time-series distributions of fluorescent calcium traces from a given neuronal spike-train. Existing fluorescent trace simulators [5,9] do not always account for the intricate features which exist in real calcium imaging datasets. In contrast, our proposed methodology enables a realistic simulation by integrating four critical design factors: (a) a model for photobleaching of calcium indicators, (b) imaging noise model, (c) a tunable non-linear baseline dynamics, and finally, (d) inhomogeneous firing patterns, with neurons that alternate between low basal and burst firing, which integrate refractory time periods due to the depletion of the calcium indicator reservoirs. Existing deconvolution techniques may be categorized in three major groups : 1) deterministic methods that compute the instantaneous spiking rates together with the estimated calcium dynamic of the fluorescent indicator [10]; 2) probabilistic approaches that compute the convolved time series of spikes that maximises the likelihood of observed data [5]; and 3) machine learning approaches [6]. We study the robustness of four representative methodologies and evaluate to which extent they can be adapted on heterogeneous simulated signals with strong underlying dynamic. Finally, we show that the results provided by the simulator, and those obtained within experimental dataset from Hydra Vulgaris are correlated.

Method

Our simulator can be decomposed into four main parts, shown in Fig. 2 : 1) the modeling of the spike impulse train using the Firing Instantaneous Rate (FIR) function; 2) the modeling of the calcium trace from dirac impulse signal using a kernel learned from experimental dataset; 3) the modeling of the fluorescence trace from the calcium by handling the photobleaching effect through a Poisson-Gaussian noise; 4) an additive non-linear baseline integration. Each step is detailed hereafter.

To model the firing pattern s(t) of a single neuron, we implement a Poisson process with a FIR describing the probability at each time step that the neuron fires. The FIR is modeled as a piece-wise step function whose high constant steps correspond to burst events, and low steps to the background spontaneous firing activity.

FIR(x) = n i=1 γ i χ [ti,ti+di] (x). ( 1
)
Here χ is an indicator function, equals to 1 when x ∈ [t i , t i + d i ] and 0 otherwise. The firing-rate γ i follows two regimes (bursting event or not), and is modeled with a truncated Gaussian law to keep frequencies positive. The number of steps n is linked to the burst-rate parameter. The burst duration d i follows a Poisson law, and the spike burst locations t i are drawn uniformly on the available time interval. The spike impulse 2 2. A VERSATILE SIMULATOR FOR SPIKE INFERENCE TECHNIQUE THAT EMBEDS HYDRA signal s(t) is then derived from the FIR using an adaptation of inhomogeneous Poisson process simulation by thinning [11]. An exponential course of return to the equilibrium is used after each firing (the observed collapses of the FIR shown in Fig. 2) to model calcium reservoir depletion.

c(t) = s(t) * k(t) with k(t) = Ae -( t τ D ) β 1 -e -t-µ τ R (2)
Here τ D is the time constant of the calcium concentration return to steady-state, β > 0 is a power law, µ is the median time of calcium increase after the electrical spike. τ R is the corresponding rising time constant. In our simulations, we calibrate the convolution kernel with data extracted from Hydra calcium traces. In the third step, we model the photo-bleaching of calcium indicators with a mono-exponential decrease [12]

λ(t) = c(t)e -t
τ , with τ being the photo-bleaching time constant. To account for the Poisson shot noise of microscopes, we further model the recorded signal as a Poisson process such that P (t) ∼ P oisson(αλ(t)), where α is the gain of the microscope [13].

To obtain a realistic calcium trace, we add a Gaussian noise G ∼ N (m, σ 2 ) with constant mean m and standard deviation σ to the recorded Poisson signal (mixed Poisson-Gaussian representation) and we also add a periodic deterministic baseline B(t) = Asin(2πf t) where A and f are tunable amplitude and frequency respectively. Finally, the fluorescence calcium trace y(t) is derived from:

y(t) = P (t) + G + B(t) (3)
The photobleaching time constant and additive Gaussian noise parameters are fitted to experimental datasets using least square method. In our simulations, the SNR of generated calcium traces is modulated using the gain α of the microscope detector [13], an increased gain leading to higher SNR.

Results

Benchmarking Deconvolution Algorithms

Using synthetic fluorescence traces, we benchmark four despiking methods (two of them are issued from the CaImAn library [9]): 1) Deterministic OASIS [14]; 2) Deterministic CDfoopsi (Constrained Foopsi) [10]; 3) Probabilistic MLspike [5] and 4) a Naive method that consists of smoothing the signal with a wavelet thresholding, before computing the first derivative of the signal and estimates the spikes locations with derivatives greater that one standard deviation of the derivative over the entire calcium trace. After having benchmarked the methods using different pre-processings proposed in the literature such as lowfrequency filtering and normalization [15] or polynomial detrending [16], we conclude that the efficiency of a pre-processing method depends on each deconvolution algorithm, and we choose to perform a polynomial detrending and data normalization for OASIS, and no preprocessing for CDfoopsi and Mlspike.

Comparing the accuracy and robustness of methods

Deterministic methods (OASIS and CDfoopsi) do not handle non-linear baseline. Therefore, high baseline amplitude A will lead to false spike burst detections. To tackle this technical issue, we filtered the inferred spikes with respect to their estimated amplitude using a user-defined Decisional Amplitude Threshold (DAT). The DAT is the only parameter we optimize in tested deterministic methods since the others have been exhaustively analysed in previous study [7]. For Naive thresholding method, the only parameter is the threshold (typically a multiplicative factor of the derivative standard deviation). The probabilistic method MLspike that concomitantly estimate the spike locations and the underlying fluorescence baseline presents four main hyper-parameters: the relative amplitude a of the spikes, the decay time constant τ of the calcium fluorescence trace after a spike, the a priori level of Gaussian noise σ and the drift d of the estimated stochastic baseline. To compare the performances of the different deconvolution methods over a large range of simulation parameters, we compute for each method the rate of missing and false positive spike detections and summarize detector performances using the Error Rate (ER) indicator proposed in [5].

3

We first evaluate the accuracy of each method by varying the different parameters of our synthetic simulator (see Fig. 3) : the gain α that modulates the SNR, the rate of burst events, the spiking rate inside bursts and the baseline amplitude A and frequency f (n = 10 simulations per set of simulation parameters). The method hyperparameters are calibrated using a grid-search which minimizes the error rate with the available ground-truth. Each method calibration is performed once for a reference set of simulator parameters: Gain of the detector α = 1, Firerate In Burst f IB = 0.75 Hz, Firerate Out of Burst f OB = 0.15 Hz (firerates correspond to γ i in/out burst regime respectively see Eq. 1), Burstrate B = 0.02 Hz, Baseline Amplitude A LF = 50, Baseline frequency f LF = 2e -3 Hz (see Eq. 3). These values have been chosen to make the simulator relevant with the experimental Hydra Vulgaris dataset. We find that, after hyperparameter calibration, the probabilistic MLspike is overall the most accurate method for each set of parameters. However, we observe that MLSpike performance rapidly degrades as simulation parameters change, especially for the gain α (SNR) and baseline amplitude. The other methods are overall much more robust to parameter variations while being less accurate for the specific set of calibration parameters. We highlight that the naive derivative thresholding produces good results for isolated spikes but can not properly handle bursts. CDfoopsi method [10] seems to provide the best compromise between accuracy and robustness to parameter variations.

Correlating simulation results with Hydra experimental dataset

The trade-off between accuracy and robustness of a despiking method is important for experimental in vivo applications as different neurons can present different firing patterns (isolated or burst) and SNR for example. We therefore compare despiking methods using experimental dataset composed by ∼ 250 single neurons imaged over 2000 frames at 10 Hz inside freely-behaving animal Hydra [17]. This dataset was obtained thanks to the robust monitoring of single neuron activity with tracking algorithm EM C2 [4]. Method parameters (DAT for deterministic methods and hyper-parameters for MLspike) were calibrated by manually labelling spike locations on real fluorescence traces, and by minimizing the distance between the inferred spikes and the manual annotation. Hydra calcium traces are heterogeneous in terms of baseline, noise and spiking dynamics, with neurons spiking quite sparsely and other in bursts (Fig. 4, [17]). As observed with previous simulations, the second deterministic method [10] is the most robust to calcium trace heterogeneity and presents the best average Error Rate compared to manual labelling (OASIS: ER avg = 0.40, CDfoopsi : ER avg = 0.18, Naive : ER avg = 0.28, MLspike : ER avg = 0.42) The average poor performance of MLspike is likely due to its poor robustness to changes in experimental conditions (see MLspike in Fig. 5).

The lack of ground truth for validating the deconvolution methods on experimental dataset imposes to analyse the consistency of the inferred spike locations. This evaluation hinges on biological assumption of existing neuronal ensembles (e.g. RP and CB neural networks in Hydra [17]) that are supposed to provide correlated spiking patterns. Also, a sufficient and balanced average number of spikes per trace should reflect the robustness of each method to heterogeneous baseline variation and noise level. In this regard, the spike distributions over the dataset are summarized in the raster-plots (see Fig. 5-a). We also calculate the average number of spikes per neurons (OASIS : 12.61 ± 4.43, CDfoopsi : 22.62 ± 14.70, Naive : 19.22 ± 4.47, MLspike : 17.51 ± 26.10) and the neuron pair-wise correlation matrices using a Jaccard distance relaxed by a time tolerance (see Fig. 5-b

). d J (n i * , n j * ) = s(ni * )∩(nj * ) s(ni * )∪(nj * ) = TP TP+FP+FN where s(n * i ) is the binarized impulse spike signal inferred for neuron n * i .
The heterogeneity in spike distribution varies between the methods. MLSpike computes many false spike detections (horizontal lines in raster-plot combined with a high standard deviation of the average number of spikes per sequence), while the naive thresholding method misses more spikes, but provides more balanced impulse spike signals in terms of average number of spikes per fluorescence trace. We observe that the highest global average Jaccard distance is obtained for the naive thresholding and CDfoopsi methods (see Fig. 5-b). Finally, we conclude that CDfoopsi appears as the best trade-off solution since it provides the best reconstruction results compared to our manual labelling, a sufficient and balanced number of spikes per trace and one of the highest global Jaccard correlation. The results obtained on experimental Hydra Vulgaris dataset are congruent with the ones derived from the proposed simulator which enforces the conclusion drawn by [7] about the versatility of CDfoopsi, especially within dataset with heterogenous baselines.

Conclusion

In this paper, we provided a mathematical framework for a fluorescence trace simulator to objectively compare and validate spike deconvolution techniques. State-of-the-art deconvolution methods were benchmarked on synthetic and experimental datasets. Such comparative analysis is necessary to account for complex underlying biological processes and the lack of ground-truth in neurosciences. We even argue that a quantitative benchmarking of methods on synthetic data without a qualitative and statistical analysis on experimental dataset make them unsuitable. Future efforts will focus on improving our simulator by modeling neuron ensemble effects, and go further in an automatic pipeline to break the neural code.

SODA -A spatial statistic framework for neuron coupling estimation and behavior prediction

The observation of physical phenomena often goes through the recording of discrete time series of events, that can be represented with marked point processes. The robust estimation of the correlation between two point processes can, therefore, unveil physical mechanisms underlying the observed phenomena. However, the robust estimation of correlation between two, or more, point-processes is hindered by the signal noise (leading to false and missing point detections), the important density of points, and possible time-shift between coupled points. We propose a statistical framework that uses hypothesis testing to estimate coupling between time point-processes. Using simulations, we show that our framework accurately estimates the coupling between two time point-processes even for noisy signal (with false point detections), for high density of points and in the presence of a time shift between coupled points. By applying our statistical framework to the recordings of neuron population activity in mouse visual cortex, we measure the functional coupling between individual neurons, and cluster neurons into functional ensembles.

Introduction

In many fields, the observation of physical events can only be done through discrete time series of events. This is the case, for example, of volcanic eruptions [1] and earthquakes [2] in geology, molecules arrival and departures form specific processes sites in cellular biology (e.g. endocytosis and pathogen entry [3], and neuronal activity through the firing of action potentials from individual neurons [4] which the study case of the experimental part of our article. The statistical characterization of relations (coupling) between two (or more) time series of events can unveil important mechanisms that underly the observed processes. For example, the observation of the sequential arrival of molecules at endocytic sites with fluorescence microscopy helped to unravel the mechanisms of cell entry [3]. Another example is given by the observation of the firing of individual neurons within a population that provides information on neuronal communication and coding [5].

Observed time series of events can often be modeled as marked point processes [6], with the point being the time location of an event and the mark its attributes (e.g. intensity, color, duration...). Therefore, the characterization of the correlation between different time series of events reduces to the estimation of the coupling between time point-processes. In the case of neuronal activity studies, the two main methods used to estimate the coupling between point-processes are either based on the estimation of underlying firing rates (i.e. the intensity of associated point-processes) [7] and the analysis of the correlation between estimated intensities, or coupling estimation is directly performed with colocalization analysis between discrete time point-processes [4]. While the methods using the estimation of firing rates are more robust to the missing and false detections of single point events, they depend on the robust estimation of firing rates and are not well-suited for detecting the synchronization of single point events. On the other hand, the colocalization methods are sensitive to false and missed point detections. Moreover, high point density can lead to fortuitous point colocalization and overestimation of processes' coupling, whereas time shifts between coupled points can lead to coupling underestimation.

To tackle these technical issues in colocalization analysis, we developed a statistical method to robustly evaluate the coupling between time point-processes, even in the presence of a time-shift between individual coupled points. Our method uses the multi-distance Ripley's K function [8] to measure the time-shifted accumulation of points from one time point-process relatively to the other. It corresponds to an adaptation of state-of-the-art statistics of point-processes to 1-dimensional temporal case of spiking events. This method has been developed to account for false inferred spike detections and potential time-shifts that results from spike deconvolution methods in calcium imaging [9]. As our method is based on statistical characterization of the Ripley's K function and hypothesis testing, it is robust to noise (false point detections) and remains accurate even for high point density. We assess the robustness of our method with synthetic simulations, and show that it outperforms state-of-the-art colocalization metrics. Finally, we use our method to compute the functional relations within a population of neurons in the visual cortex of mice, from their individual spiking activity. 1 2 Method

Measuring the coupling between time point-processes with the Ripley's K function

Ripley's K function introduced by Brian Ripley in the 70's [8] remains the gold-standard to measure the coupling between spatial point processes.

The idea of our temporal adaptation is to measure the coupling between two time point-processes s 1 = [t 1 , . . . , t n1 ] and s 2 = [t ′ 1 , ..., t ′ n2 ] by creating a regular mesh grid whose topology is dependent on s 1 point locations in time and by comparing the effective distribution of s 2 points falling in each mesh to the distribution expected under a random distribution assumption through a statistical hypothesis rejection test. Thus, the Ripley's vector G N = [G 0 , . . . , G i , . . . G N -1 ] embeds a 1-dimensional mathematical mesh grid implementation composed by an ensemble of N fixed size rings of radius [r i , r i+1 ] centered around s 1 points. Thus, Ripley's function just corresponds to a statistical effective counter of s 2 points that fall into rings around s 1 objects with a boundary corrective term w which corrects for the potential underestimation of neighbors to points that are close to the starting-and ending-points of the time study period. This correction is inspired by the 2-dimensional corrective term used by Ripley in [8].

G i = |Ω| n 1 n 2 t k ∈s1 t ′ l ∈s2 1 {ri≤|t k -t ′ l |≤ri+1} w(t k , t ′ l ) (1)
with

w(t k , t ′ l ) = 1 + 1 {|t k -t ′ l |>|t k -δΩ(t k )|} ( 2
)
where |Ω| is the length of the time period over which the two time point-processes are observed, 1 is an indicator function such that 1 {A} = 1 if A is True, 0 otherwise. δΩ is the coordinate of the closest boundary of the study domain ie δΩ(t k ) = 0 if t k ≤ |Ω|/2, Ω otherwise. To detect a significant coupling between time series (point-processes) s 1 and s 2 , we design an hypothesis testing approach. We compare how far is the effective counter in rings represented by the K-Ripley function from the expected number of points expected under a null hypothesis H 0 of Complete Spatial Randomness where points are located according to an homogeneous Poisson distribution for s 2 point-process. Analytical mean and standard deviation parameters are derived by calculating intersection of 1-dimensional volumes corresponding to the overlapping of rings. Under H 0 , Ripley's K function tends to a normal distribution [10] in accordance with the central limit theorem. Thus, the distribution of G = [G i ] i=0..N -1 is fully characterized by its mean

M N = [E{G i } = µ i ] i=0..N -1 and its standard deviation Σ N = [E{G 2 i } -µ 2 i = σ 2 i ] i=0..N -1 .
Using the CSR hypothesis for s 2 time points and the boundary correction (Eq. 2), we compute that

µ i = 1 n 1 t k ∈s1 y∈Ω 1 {ri≤|t k -y|≤ri+1} 1 + 1 {|t k -y|>|t k -δΩ|} dy (3)
and

σ 2 i = |Ω| n 2 1 n 2 t k ∈s1 I 21 (t k , r i , r i+1 ) + tj ∈s1,tj ̸ =t k I 22 (t k , t j , r i , r i+1 ) (4)
with

I 21 (t k , r i , r i+1 ) = -µ 2 i |Ω| + y∈Ω 1 {ri≤|t k -y|≤ri+1} 1 + 3 × 1 {|t k -y|>|t k -δΩ|} dy I 22 (t k , t l , r i , r i+1 ) ≈ - µ i |Ω| + 1 {|t k -δΩ|>ri+1} 1 {|t l -δΩ|>ri+1} × y∈Ω 1 {ri≤|t k -y|≤ri+1} 1 {ri≤|t l -y|≤ri+1} dy
The size and number of the rings provide a maximum duration beyond which colocalization can no longer be detected and a temporal resolution to distinguish two close interactions.

Statistical test of time point-processes' coupling

To build a statistical test of time point-processes' coupling, we introduce the reduced statistics

G = A -1 G -M N Σ N (5)
with A a correction matrix for ring's overlapping [11]. Under the null hypothesis of s 2 randomness, G is a standard normal vector (i.e. each of its component Gi ∼ N (0, 1). Therefore, a significantly high value of a vector component Gi would indicate an accumulation of coupled points around reference points corresponding to a positive coupling (while a low value would indicate a depletion of s 2 points corresponding to a negative coupling). Similarly to the statistical test introduced for spatial point-processes [11], we use the maximum component of reduced Ripley vector Gmax = sup 0≤i≤N -1 Gi to test if there is at least one ring [r i ; r i+1 ] where s 2 time points accumulate significantly. To compute the p-value associated with the observed maximum component Gmax , we compute that, ∀x > 0,

Pr{ Gmax ≥ x} = 1 -Pr{∀i ∈ [1..N -1], Gi < x} = 1 -(Pr{N (0, 1) < x}) N = 1 -cdf N (x),
where cdf(x) is the cumulative density function of the standard normal law. Finally we obtain the pvalue

p-value = 1 -cdf N ( Gmax ) (6)

Quantitative characterization of time point-processes' coupling

To further characterize the putative coupling between two time point-processes, we determine the components of Ripley's reduced vector G that are significantly high by using the universal threshold T (N ) = 2 log(N ) [12], which is widely used in image processing to determine the significant component of a signal corrupted with standard white noise. Thus, Gi > T (N ) indicates that there is a significant accumulation of s 2 points at a time shift comprised between r i and r i+1 from s 1 points. This allows the detection of coupling between two time point-processes at different distances and not only the detection of points co-occurence at the same time. Hence, our framework can handle coupling estimation with time shifts and varying delays. Finally, we convert the reduced Ripley vector components into a coupling probability between all the points (t k , t l ) of time point-processes s 1 and s 2

P (t k , t l ) = N -1 i=0 1 {ri≤|t k -t l |<ri+1} σ i Gi 1 { Gi>T (N )} G i (7)
3 and extract a global coupling metric on the entire time series

GC(s 1 , s 2 ) = 1 n 1 t k ∈s1 1 n l t l ∈s2 P (t k ,t l )̸ =0 P (t k , t l ) (8)
where n l is the number of t l coupled points in s 2 , i.e. such that their coupling probability with t k in s 1 is strictly positive.

3 Results

Synthetic simulations

To validate our proposed statistical framework we use simulations where the coupling characteristics between time point-processes are known, and we compare the results of our method to standard measures of signals' correlation.

Robustness to variations of coupling level and time shift

Using simulations of time point-processes with varying coupling level and time shifts, we compare the accuracy of our statistical method with standard correlation metrics (see Appendix for details): 1) the Pearson correlation coefficient, 2) F1 score and 3) the Cluster Core Index (C3I) [13]. For increasing simulated levels of coupling and several time shifts ((µ s , σ s ) = (0, 0), (1, 0), (0, 0.3), and (1, 0.3), we compare our statistical method with other standard methods (Figure 2). For µ s = 0, i.e. when most coupled time points are co-occurring in same time steps Ω t , all methods provide a correct estimation of the effective coupling between time point-processes. We highlight that the non-linear, convex shape of C3I curve is due to the used index of correlation [13]. However, when the mean time shift increases to µ s = 1, our statistical framework and C3I outperform classical correlation coefficients that are not well-equipped for time-shifted coupled point-processes. We conclude that our framework can handle coupling with or without time-shifts between coupled points and correctly estimates the simulated level of coupling.

Robustness to false and missed detections

A common issue when dealing with point-processes is the presence of false and missed detections in time point series. These artefacts are due to the presence of noise in the acquired biological signal. For example, when analyzing the spiking activity of individual neurons with calcium imaging, noise in calcium intensity traces and imperfections in deconvolution algorithms used for extracting neuron discrete spikes (time points) lead to false and missed detections (see section 3.2). In Figure 3, we measure the robustness of our method to increased levels of false and missed detections. We plot the F1 score of the proportion of couples estimated with our method. For any simulated level of false detections (0%, 10%, 30% and 50 %), our method estimates an accurate number of points' couples (> 80% of couples are recovered). Conversely, our method is much more sensitive to missed detections. Indeed, an increased number of missed detections reduces the number of observable couples. Typically, for a percentage 0 ≤ p ≤ 1 of missed points in processes s 1 and s 2 , the expected proportion of missed couples will be equal to 1p 2 . Finally, we also test the sensitivity of our method to the intensity λ 0 of underlying point-processes. Indeed, an increased intensity (i.e. overall number of points over the observed period ∼ λ 0 |Ω|), can lead to an increased number of points' couples due to chance. We observe that, actually, our statistical method is very robust to an increased intensity of underlying point processes from λ 0 = 0.1 Hz to 1 Hz.

From previous simulations, we conclude that our method outperforms other standard correlation methods when there is a time-shift between the coupled points of the two time point-processes, i.e. when coupled points are not necessarily co-occurring in the same time steps. Moreover, we assessed the robustness of our method when false detections (i.e. random points) are added to the original coupled time point-processes, as well as when the intensity of processes is overall increased. The only sensitive parameter is the level of missing points that leads to an expected decreased of detected points' couples.

A

Functional coupling between individual neurons

We apply our method to measure the functional coupling between individual neurons from their monitored spiking activity. We use the online dataset from [14], corresponding to two-photon imaging of neuron activity in mouse visual cortex (File M1d1AS in the dataset). From calcium fluorescence traces, the exact spiking times can be obtained using spike inference techniques with variable accuracy and robustness [9]. A representative rasterplot of neuron spiking activity obtained with the constrained FOOPSI deconvolution algorithm [15] is shown in Figure 1-A (|Ω| = 5 minutes, image acquisition rate = 12.3 Hz). We measure the coupling between individual neurons using the coupling index (eq. 8) (Figure 4-A). Size of used Ripley's vector (eq. 2) is equal to N = 4, with identical time-shifts r i+1r i = 1 frame (∼ 80 ms) for i = 0..3. Proportions of coupled spikes for the different time shifts are respectively equal to [65%, 13%, 12%, 10%], meaning that 2/3 of the coupled spike times are co-occuring in the same time step Ω t , while 1/3 are timeshifted by more than one frame. These latter, time-shifted coupled spikes are completely missed by standard correlation techniques. Using the coupling information between individual neurons, we represent neurons' couples with an undirected network graph, with edges corresponding to strictly positive coupling indexes (Figure 4-B). We identify neuronal ensembles with a Louvain clustering algorithm [16] which maximizes graph modularity, and obtain 3 neuronal ensembles (with n ≥ 10 neurons), in agreement with [14].

We conclude that our statistical framework allows the robust estimation of functional coupling between individual neurons from the calcium imaging of their spiking activity. Contrary to standard correlation metrics, it allows the estimation of points' coupling even in the presence of time-shifts, and does not require any thresholding of correlation coefficients for the network representation and clustering of neuronal connectivity.

4 Conclusion

We have proposed a statistical method for estimating the coupling between time point-processes that use the multi-distance Ripley's K function and hypothesis-testing framework. Our method is able to accurately estimate the coupling between time-shifted correlated point processes, and is robust to high intensity of point processes and false detections. The unique ability of our framework to compute the coupling between time-shifted point-processes is used to quantify the functional coupling between individual neurons imaged with fluorescent calcium indicators in mouse visual cortex [14].

Appendix

Pearson correlation coefficient For a discretized time period Ω = {Ω t } 1≤t≤T , we introduce the indicator functions δ t (s k ) = 1 {∃t l ∈ s k |t l ∈ Ω t } for each time-process {s k } k=1,2 that determines whether at least one point of each time point-process falls into the time step Ω t . The Pearson correlation coefficient is then given by

r(s 1 , s 2 ) = 1≤t≤T δ t (s 1 )δ t (s 2 ) -T s1 s2 1≤t≤T δ t (s 1 ) -T s1 2 1≤t≤T δ t (s 2 ) -T s2 2 , with sk = T -1 1≤t≤T δ t (s k ).

F1 score

To compute the F1 score between time point-processes s 1 and s 2 , we set a tolerance for point-matching (tol = 2 frames) and define true positive (TP) as T P = 1≤t≤T δ t (s 1 )δ t±tol (s 2 ) with

δ t±tol (s 2 ) = 1 ∃t l ∈ s 2 |t l ∈ t-tol≤j≤t+tol {Ω j } .
False positive (FP) and false negative (FN) are respectively given by

F P = 1≤t≤T δ t (s 2 )(1 -δ t±tol (s 1 )) and F N = 1≤t≤T δ t (s 1 )(1 -δ t±tol (s 2 ))
. Finally F1 score is equal to

F 1 = 2 × T P 2 × T P + (F P + F N ) . 7

DAT -Automatic thresholding for spike estimation

Monitoring the activity of multiple neurons with high spatio-temporal precision in living animals remains a gold standard in neuroscience. The recent advances of fluorescence imaging of neurons' spikes using genetically-encoded calcium indicators and calcium probe design require significative improvements in the post-processing of calcium fluorescence traces to accurately extract spikes. Analyzing neuronal spike trains from time-lapse calcium fluorescence sequences appears as an inherently unsupervised inverse problem, compounded by low signal-to-noise ratios, slow indicator conformation changes, and unpredictable baseline fluctuations. State-of-the-art deconvolution methods exhibit variability and parameter sensitivity, necessitating the manual calibration of spike inference techniques to make the distinction between meaningful spikes conveying information and spikes resulting from artifact or noise deconvolution. To address this issue, a statistical framework inspired from point process theory is derived to calibrate the parameters of spike deconvolution techniques robustly and automatically. The proposed functional connectivity metric is validated through simulations of convoluted spiking patterns and apply to in vivo two-photon calcium imaging of single neurons' activity in the mouse visual cortex. Such approach highlights a significant coupling between individual neurons and allows for the identification of correlated neuronal ensembles paving the way for more accurate and objective analyses of neuronal network dynamics in living animals.

BINOE -Bayesian Inference framework for Neuronal Overlapping Ensembles

"Groups of neurons repeatedly firing together" intrinsically constitute the functional unit of neuronal computation. Input stimuli and spontaneous behaviors are encoded by sequential activation patterns of these co-active neural ensembles whose extraction is the cornerstone in understanding and breaking the neural doctrine. Bayesian inference appears as a promising statistical framework to accurately and robustly estimate such overlapping ensembles inherently characterized by randomness. Our non-supervised clustering technique based on synchronous activation patterns handles for the first time overlapping and gets rid of the arbitrary thresholds commonly used by classical graph theory approaches. It provides further statistical evaluators to quantify neuron ensemble relevancy and mutual interactions. The results of our analysis demonstrates its efficacy in improving identification of neural clusters and correlation with with stimuli and suggest its potential in advancing our understanding of neural ensembles.

Introduction

Monitoring the activity of neurons in free-behaving animals is a fundamental prerequisite to decipher the neural code. Recent advances in microscopy and genetic engineering have made calcium fluorescent microscopy techniques the gold-standard methods for monitoring neuron population activity at the individual cell level in small organisms [1]. These techniques allow for the identification of "neuronal ensembles," corresponding to specific groups of neurons that repeatedly fire together [2] in relevant time windows. These neuronal ensembles of coactive neuronal cells are assumed to be the fundamental units of the neural doctrine and the intrinsic substrates of brain computation and neural states. Their robust identification from neural spiking activity and the correlation of the activity of the ensembles with mental states of the animal, external stimulations or spontaneous behaviors could make the understanding of the neural code a step further.

However, going from fluorescent microscopy images to a robust identification of neuronal ensembles requires several important steps. In the in vivo neuroscience context, neuronal ensemble detection methods are generally extracted from the resulting neuronal spiking activity matrix. The robust extraction of the exact spike time emission dates from raw fluorescence traces is a critical step which is still not trivial despite the huge diversity of methods implemented to fulfill that task [3] [4] [5] [6] [7] [8]. Besides, several attempts of benchmarking have been carried out [9] [10] [11] but they do not identify a gold-standard method and demonstrate that method's generalizability on real experimental data is elusive due to the lack of ground truth.

Then, from the estimated spiking activity, various methods have been proposed to detect ensembles inspired from social networks, economy or politics. Many reviews [12] [13] [14] [15] [16] try to exhaustively summarize this abundant literature. These methods are, mainly, divided into three categories: graph theory techniques [17], spectral techniques [18] [19] [20] [21] and model-based statistical techniques [22]. Graph theory techniques involve building a neuron pairwise connectivity matrix based on a functional connectivity metric and then applying on it a projection and a thresholding to unveil a resulting topological network. On this graph, an optimization algorithm based on a specific cost function (Modularity optimization, edge affinity. . . ) is applied. Spectral techniques involves applying algebraic frameworks on spiking activity matrix to extract similar repeating patterns using euclidean distances or similarity metrics. Data are usually projected in mathematical subspaces to reduce data dimensionality and orthogonal components are extracted using various machine learning methods like principal component analysis, independent component analysis or extensions... Model-based statistical techniques formulates an explicit analytic model of the interactions between neuronal ensembles to robustly identify them through inverse problem theories.

Estimating neuronal ensembles is hard since a lot of user-defined choices are performed about functional connectivity metrics, thresholds or neuronal ensemble definitions making the comparison of the methods hard. The intrinsic problem of ensemble detection is not even straightforward since a neuronal ensemble is characterized by (a-)synchronicity, activation randomness, and overlapping. Indeed, each time an ensemble 1 is activated only a subset of its neurons are co-recruited leading to a varying group response. In addition to this, the neuron basal firing induces a spontaneous non-significant activity of each neuron outside the activation of its specific group inducing a-synchronization events. Finally, neuronal ensembles are assumed to be functionally overlapping i.e that a neuron can encode several tasks or mental representation by firing with several groups simultaneously or at deferred moments.

We developed a Bayesian inference framework to robustly extract neuronal ensembles from their estimated spiking activity. This model provides an explicit biologically-inspired definition of neuronal ensembles based on synchronization. In addition, such probabilistic model handles the possible overlapping and the partial asynchrony of neurons with their associated ensembles. The extraction of additional statistical features to correlate ensemble activity with external stimulations, of core neurons with specific roles for the synchronization [18]. The exhaustive benchmarking of the statistical framework is performed, firstly, on synthetic data to highlight its accuracy and robustness over a large range of simulated parameters. Then, we apply our framework to a real dataset of zebrafish larval optic tectum got from [16]. Using statistical spatial analysis [23] [24] [25], we correlated neuronal ensemble activity patterns with stimulations allowing a better interpretation and even a refining of the inferred ensembles.

Overlapping community based on synchronous events 2.1 The Statistical Model

Bayesian inference aims at clustering neurons in functional ensembles based on their co-activity over time. Gold-standard imaging method to monitor the activity of hundred of neurons in living animals with single cell and high temporal resolution is based on calcium fluorescence microscopy techniques [1]. From the recorded movies we extract calcium fluorescent traces using motion correction, detection or even tracking image processing techniques. The raw intensity corresponds to a count of photons got on the camera detector that gives a graded evolution of fluorescence over time. The robust extraction of patterns corresponding to action potential emissions still remains an open question especially for non-flat fluorescence signal, Non-Homogeneous Poisson action potential emissions, time-varying patterns, noisy data... Spike inference techniques overcome or ignore this issues and extract neural activity in a binary matrix S = {0, 1} N ×T .

Then clustering neurons into overlapping functional ensembles using a robust statistically framework is fundamental. Bayesian inference turns out to be well-suited to handle this kind of structured problem that encounter emergent properties from individual behaviors. Following pioneer work of Diana et al [22], we present a framework to cluster N individual neurons into A overlapping ensembles that activate over time in a W binary matrix. Handling ensemble overlapping combined with a random neuron activation makes synchrony based clustering not trivial at all.

Let A be the number of communities, N the number of individual cells, T the number of recording timesteps or events, fixed and known parameters of the model. Let (S) N ×T be the input binary rasterplot matrix of the N neurons over time. Thus, ∀i ∈ [1, . . . , N ], ∀t ∈ [0, . . . T ], S i,t = 1 if neuron i has been activated in timestep t, 0 otherwise.

Each ensembles will recruit neurons with a fixed probability α k , ∀k ∈ [1, . . . , A] such that a neuron can belong to several communities. A latent matrix (Z) N ×A describes the membership of each neuron to each assembly through Bernoulli distributions.

∀i ∈ [1, . . . , N ], ∀k ∈ [1, . . . , A], z i,k ∼ Bern(α k ). Z i,k =
1 if neuron i belongs to ensemble k, 0 otherwise.

2

The probability that an ensemble fires at each timebin is constant over time and fixed for each assembly. It follows a Bernoulli distribution with probability p k . ∀k ∈ [1, . . . , A], w t,k ∼ Bernoulli(p k ). It defines a latent variable matrix (W ) A×T corresponding to the activities of ensembles at each time.

W t,k = 1 if ensemble k fires at timebin t, 0 otherwise.

The firing activity of a neuron at time t is a function of the combination of sets of ensembles to which it belongs and which are activated at that specific moment. We define, thus, neuron activation through a conditional probability given membership and activation which follows a Bernoulli distribution.

G i = {k ∈ [1, . . . , A]/z i,k = 1} is the set of ensembles, neuron i belongs to of dimension K i = card(G i ).
Obviously, the number of combinations of sets of ensembles that fully define the activation of a neuron depends on the number of ensembles K i neuron i belongs to. An ensemble activation pattern z (t) corresponds to a sub-combination of ensembles, activated at time t, of the whole set of ensembles neuron i belongs to. Then, the probability for the neuron i to be activated at t is :

∀z (t) ∈ R Ki , λ Gi (z (t) ) = P (s i,t = 1|w t,g1 = z (t) 1 , . . . , w t,gK i = z (t) Ki ) (1)
Finally, the statistical model is defined by observed and latent variables X = (Z, W, S) and statistical model parameter vectors θ = (α, p, λ) where α = (α 1 , . . . , α A ), p = (p 1 , . . . , p A ), λ = (λ g1 , . . . , λ gN ). The function that summarizes the interaction between the variables and describes the probability of consistency of the statistical model with the binary observations is the true likelihood function :

P (S|Z, W, θ = (α, p, λ)] = N i=1 A k=1 α z i,k k (1 -α k ) 1-z i,k (a)
× A k=1 T t=1 p w t,k k (1 -p k ) 1-w t,k (b)
× N i=1 T t=1 λ Gi (w t,Gi ) si,t (1 -λ Gi (w t,Gi )) 1-si,t (c)
(a) : Membership term. α k : recruitment probabilities, Z membership binary matrix. (b) : Ensemble activity term. p k : ensemble activation probabilities, W ensemble activation binary matrix (c) : λ Gi Conditional spiking neuron probabilities given a membership and an ensemble activation. Given the binary matrix of neural activation, we solve a statistical optimization problem. We seek to maximize the log-likelihood function which corresponds to the probability of observations given the model. The Maximum Likelihood Estimator (MLE) derives partially the log-likelihood according to all the parameters to zero it. This procedure cannot, however, be carried out here when the model is defined in very high dimensions for vectors of parameters and matrices of non-measurable latent random variables. A Bayesian statistical method called Gibbs Sampling algorithm estimate an approximation of the model posterior distribution to estimate the parameters. This algorithm requires the exact calculation of the posterior probability distributions of the statistical parameters of the model and the latent variables in order to sample iteratively by this MCMC method.

Inference of model parameters

Priors

To calculate the posterior distribution, it is necessary to define a prior on the statistical parameters sought. There are many possible choices (Jeffrey's prior, Uniform prior, Haldane's prior...) Since the Gibbs sampling routine relies on our ability to sample posterior distributions, it is therefore necessary to have known and samplable distributions as output. To do this we use conjugate priors. As all the parameters have similar Bernoulli distributions, Beta priors are used in this article because the posterior distribution is then a beta distribution itself, perfectly known and easily samplable.

∀k ∈ [1, . . . , A], α k ∼ Beta(α * k (α), β * k (α)) ∀k ∈ [1, . . . , A], p k ∼ Beta(α * k (p), β * k (p)) ∀k ∈ [1, . . . , A], λ Gi ∼ Beta(α * Gi (λ), β * Gi (λ))
The α * , β * 's are the hyper-parameters of the model i.e. the parameters of the prior distributions. Their values must be initialized so that the distributions favor probabilities consistent with the variables of our problem.

P (α) = A k=1 α α * k (α)-1 k (1 -α k ) β * k (α)-1 B(α * k (α), β * k (α)) P (p) = A k=1 p α * k (p)-1 k (1 -p k ) β * k (p)-1 B(α * k (p), β * k (p)) P (λ) = N i=1 T t=1 λ Gi (w t,Gi ) α * k (λ)-1 (1 -λ Gi (w t,Gi )) β * k (λ)-1 B(α * k (λ), β * k (λ))
To derive the classical Gibbs sampler, we need to calculate the posterior distribution of each parameter conditioned on every other ones and observations P (θ i |X, θ ∀j̸ =i ) and be able to recognize samplable distributions.

To do so, we just can get rid of all terms considered as fixed by the conditional formulation.

α l |Z, W, S, p, λ, α ∀k̸ =l ∼ Beta(G l , Ḡl ) (2)
with

G l = α * l (α) + N i=1 z i,l and Ḡl = β * l (α) + N i=1 (1 -z i,l ) p l |Z, W, S, α, λ, p ∀k̸ =l ∼ Beta(H l , Hl ) (3)
with Then we need to compute the conditional distributions of the latent variables W, Z which are not, strictly speaking, statistical parameters of the model but which will nevertheless have to be sampled iteratively in the same way in the Gibbs algorithm.

H l = α * l (p) + T t=1 w t,l and Hl = β * l (p) + T t=1 (1 -w t,l ) λ z c |Z, W, S, α, p, λ z ′ ̸ =z c ′ ̸ =c ∼ Beta(T z c , T z c ) (4)
P (w tk = 1|Z, W -tk , S, θ) = 1 1 + ρ tk (5)
such that

ρ t,k = ( 1 p k -1)× N i=1 i∈k (λ Gi (w t,Gi )| w k,t =0 ) sit (1 -λ Gi (w t,Gi )| w k,t =0 ) 1-sit (λ Gi (w t,Gi )| w k,t =1 ) sit (1 -λ Gi (w t,Gi )| w k,t =1 )) 1-sit
∀ξ ∈ {0, 1}, λ Gi (w t,Gi )| w k,t =ξ = P (s it = 1|w t,g1 = z 1 , . . . , w t,k = ξ, . . . w t,gK i = z Ki ) is the conditional probability associated with the activation of the groups K i to which neuron i belongs by setting the activity of the k-th group to 0 or 1.

P (Z ik = 1|Z-ik, W, S, θ) = 1 1 + η ik (6) 5 5. BINOE -BAYESIAN INFERENCE FRAMEWORK FOR NEURONAL OVERLAPPING ENSEMBLES η ik = ( 1 α k -1)× T t=1 λ Gi\k (w t,Gi\k ) sit (1 -λ Gi\k (w t,Gi\k )) 1-sit λ Gi (w t,Gi ) sit (1 -λ Gi (w t,Gi ))

Improving convergence

The implementation of a Gibbs sampler aims at exploring by a random walk the space of the parameters in order to estimate the posterior distribution of the model. In a very high dimensional space, we cannot guarantee the convergence to the global maximum. Strategies of multiple random initializations, simulated annealing, burning period of sample rejection try to avoid local maxima. However, higher is the dimensionality of the model, more likely it is to fall into a local maximum. This situation happens when the number of communities A increases. Note that the number of model parameters k is :

k = A × (N + T + 2) + 3 A
Based on Carrillo-Reid's idea [26], the S neuron activity matrix can be interpreted as a matrix of activation of different neuronal populations. There is an analogy between clustering documents into topics based on word occurrence and clustering neuronal populations into functions based on neuron co-activation, idea mentionned by [27] has been implemented by Carillo-Reid [28] [19] through SVD. Using a natural langage processing community technique, we aim at intelligently initializing the ensemble activity matrix W using Spherical KMeans based on cosine similarity [29]. The goal is to find statistically significant neuronal populations that repeat identically over time but are sufficiently different from each other to suggest that they encode different functions and need to be placed in different groups. Once these populations are detected through a Z-score statistical test, they are normalized using an NLP function called TF-IDF and then clustered using Spherical KMeans based on cosine similarity to generate W 0 , the ensemble activation matrix at iteration 0. To do this, the open-source implementation proposed by [30] is used. This allows us to bring the Gibbs Sampling random search closer in the neighborhood of the global maxima. The clustering technique is firstly evaluated on a synchrony based neural activity simulator and then applied on a real data of V1 mice brain visual cortex in response do static grating stimuli.

Results

Results on synthetic dataset

Synthetic simulator based on the statistical model

The simulator part to generate neural activity is divided into 3 different parts Fig. 3. 1-Drawing of an overlapping clustering in A communities using Bernoulli trials using the set recruitment probabilities. 2-Drawing the binary activity matrix of the clusters using Bernoulli trials using the set of activation probabilities. 3-Generating the neural activity using the conditional dictionary of neural activity given the membership and activation of the clusters at the time. For each combination of memberships, for each combination of activations, the probability of neural synchronous response and spontaneous activities are defined and used to produced S. In the synthetic example presented hereafter, A = 3, T = 1000, p = [0.1, 0.1, 0.1], alpha = [0.5, 0.5, 0.5]. Neurons = [350, … , 400] 𝜆 ∶ {0, 0, 0} → 0.05, 1, 0, 0 → 0.6, 1, 1, 0 → 0.95 0, 0, 1 → 0.8, 0, 1, 1 → 0.95, 1, 1 0 → 0.95 0, 1, 0 → 0.7, 1, 0, 1 → 0.95, 1, 1, 1 → 1

Figure 4: The neuronal conditional spiking probability is defined for all possible combination of membership. For neurons belonging to several groups, a more important number of conditional spiking probabilities is expected. Every neurons is defined with the same basal activity at 5 %.

Parameter and hidden variable estimation

The iterative Gibbs sampler estimation algorithm estimates each of the statistical parameters and latent variables as it loops through, so that these quantities or the error on these quantities can be estimated over the course of the iterations.

We notice that in less than 10 iterations the error rate on the estimated activities of each ensemble is close to 0% calculated with the F1-score. The clustering accuracy is obtained using the Overlapping Normalized Mutual Information (ONMI) extracted from Macdaid that reaches 100% in only 3 iterations. The statistical parameters also converge towards the simulated target values. An averaging of the values obtained during 8 the last iterations should allow to accurately estimate their value.

Robustness algorithm analysis on synthetic dataset

We estimate the robustness of the algorithm by testing different range of simulation parameters. We want to prove the robustness for different ensemble sizes, ensemble activity rate, level of synchronization and a-synchronization given membership. to see the ability of the algorithms to robustly recovers targeted simulation parameters. To highlight the result variability, Montecarlo simulation are performed to estimate the trial to trial error variability.

• Variation of the community activation through p. The algorithm seems robust to a variation of ensemble activities from 5% corresponding to the basal activation rate to 25%. Above this value the accuracy collapses since the ensemble coactivity increases. The preprocessing step is not robust to a high degree of ensemble activity overlapping (> 40%) because it preclusters timesteps as single ensemble activation. The levels of ensemble activation the algorithm is robust to, matches perfectly with sparse real data statistics of neural activities. However, an adaptation would be easily implementable by adding overlapping in the preprocessing for timebins at the border of several clusters corresponding to semitransparent points on fig. 2-e). • Variation of the size of the community through α. The algorithm is completely robust to a variation in the ensemble sizes. • Variation in the ensemble a-synchronicity through λ 0 . The algorithm is robust to a basal activation rate below 30%. When this value becomes higher and reach 50%, accuracy deteriorates since neurons become very active outside the respective activation of their groups and the gap between asynchronization and synchronization firing statistic decreases. It becomes very complicated for the algorithm to make the distinguish between totally desynchronized neurons highly active and synchronized neurons clusters in small different groups explaining the confusion of the algorithm. For real neural activity, basal spontaneous neuron activation rate is around 5% far from this frontier. • Variation in the ensemble synchronicity through λ 1 . The algorithm is robust to level of synchronization above 40 to 50%. Around and below theses values, ensemble activation are not sufficiently coherent to distinguish ensemble coherent activation with random activation. More generally to identify ensembles efficiently their level of synchrony and asynchrony should be enough distant to allow the algorithm to make the distinguish. Otherwise, the ensemble state seems similar to the firing independant case. • Variation of the number A of communities. Modulating the number of communities A can make the algorithm very long. Since this version is mathematically exact, it learns all the possible combinations of ensemble activation given a specific membership for neuron corresponding to 3 * A different synchronization probabilities to learn. But as long as the number of communities remains low enough to limit the ensemble coactivity rate, the algorithm continues to efficiently and robustly estimate the most appropriate statistical parameters given the model. • Variation of N. This algorithm seems perfectly robust to a variation of the number N of neurons.

It almost seems that as the number of neurons increases, the uncertainty on the parameter estimates decreases a little at a fixed recruitment rate as if the increase of neurons outside the groups therefore characterized by random statistics allows to better identify the neurons synchronized with groups. • Variation of T When the number of time steps increases the recording becomes longer and more information about synchronization patterns is extracted from data. It explains why the level of errors is decreasing accordingly.

Additional experiments to assess the durability of the algorithm in response to various spike patterns, given their activity rate have been performed. Although the algorithm's model assumes independent activities 10 of the ensembles, structured activation in the form of sequential peaks does not disrupt its accuracy and robustness as long as the rate of ensemble co-activations is controlled. These findings suggest that the algorithm can reliably extract information even when faced with non-random patterns of ensemble activity

Automatic estimation of the number of ensembles

To automatically estimate the number of ensembles we have set up a dynamic scanning of the number of communities. The idea is to use a statistical criterion of reference such as the log-likelihood or the Akaike Information Criterion (AIC) to automatically increment in an ascending way the number of communities as soon as a convergence regime is reached corresponding to a plateau. As the tested value A of the number of communities gets closer to the target value A*, the plateaus of the log-likelihood criterion increase and the plateaus of the AIC criterion decrease. As soon as A > A * , we notice experimentally that the tendency is reversed, i.e. that the plateaus of log-likelihood and AIC decrease and increase again respectively. This information allows us to define a stopping criterion based on the value of the AIC or the log-likelihood in the convergence regime. On the following example, the statistical properties of the bencharmarking are kept using A * = 3 communities.

Results on experimental dataset -Larval zebrafish optic tectum

Experiment description

We apply our algorithm on a stimulus-evoked calcium imaging opensource dataset from the larval zebrafish optic tectum provided by [16]. The aim is to display light points at different angular positions in the animal's visual field. The recording of the zebrafish's neuronal response should enable us to identify the groups of neurons integrating the response to these stimuli, measure their statistical coupling to the stimuli, measure their intrinsic coherence rate and extract statistical properties about the underlying neuronal ensembles. For this purpose, we have reduced ourselves to the 5 dominant stimuli of the experiment 8, 4, 5, 9 and 7, corresponding respectively to angular positions at 120°, 60°, 75°, 135°and 75°. In order to limit the pollution of the photobleaching effect, we restrict ourselves to the study of 1184 frames located at the beginning of the recording between frames 390 and 1571 and corresponding to the sequential repetitions of 6 stimuli.

Data Pre-processing

Our algorithm works by finding repeated synchronizations in neuron activity to cluster them into ensembles. Overlap in clustering occurs when a sufficient number of asynchronizations in a statistically similar group of neurons is detected. In this case, repeated and consistent asynchronizations of a subset of neurons imposes the existence of another underlying group necessary to explain this activity. Another group must therefore be introduced, and neurons that are too often asynchronized must be placed at the intersection of two different groups. For non-overlapping clustering algorithms, the simple co-activation of neurons, often obtained using a pairwise similarity metric, is sufficient to model neuronal interactions. Introducing overlap into clustering requires the introduction of asynchrony, or the integration of statistically 'silent' ensembles to justify neurons belonging to several groups.

When neuronal activity data are polluted by activation noise, the number of synchronizations decreases, leading to a loss of coherence in synchronous neuronal ensembles and even to false asynchronizations, resulting in artefact overlaps.

To limit these issues, several pre-processing steps are applied to filter out noise in the data: 1) a convolution with a Gaussian kernel is applied to each neuron activity (row) and summed over all neurons (column). Then, a threshold is applied to keep neural activities corresponding to common activation patterns of neurons, called events or wave activations in response to the stimuli, and to remove too singular neuronal activation. 2) the algorithm works on the presence or absence of each neuron at each time step, but a temporal dimension of activation in terms of duration and delay has not been yet integrated into the algorithm's statistical model. To compensate for this, the asynchronizations measured must correspond to inter-and not intra-event asynchronizations. In other words, if a neuron is turned on during an activation wave, it is necessary to normalize its appearance in relation to the appearance of all other neurons during this event, so that all neurons are active together at every single instants of the wave. This limits the potential effects of advance or delay in their activation. Concretely, this means to filtering out insufficiently long activities of a neuron during an activation wave, detecting the beginning and the end of the activation wave and then setting the activity of all neurons to 1 during the total duration of a wave. If desired, wave widths can also be normalized, to avoid making one activation wave more important than another.

Overlapping clustering estimation and synchronization ensemble relevancy

Level of synchronization

Ensemble overlapping (Positional information)

Figure 9: Overlapping clustering visualization with the synchronization levels. The representation provides information of neuron overlapping memberships. Neuron spatial locations inform about the combination of neuronal ensembles the neuron belongs to. The size of the neuron describes the level of synchronization with its groups. In neural coding, comparing neuron clustering algorithms on experimental dataset is tough since the neural ground truth is missing systematically. Generally, the analysis make the assumption that the accuracy results got by each method on synthetic data are perfectly generalizable to experimental dataset. This working hypothesis is false because, firstly, the way synthetic data are generated largely bias the results and then the model of real world dataset is always far more complicated than simulation's and the generalization property is no longer guaranteed. A clustering comparison method on experimental data would need to account for two fundamental criteria so as to answer to the question: which algorithm for neuron clustering is the most accurate ? Firstly, the statistical coupling between inferred ensemble activities and stimuli would be crucial. Indeed, 15 once the clustering step has been performed the statistical correlation between ensemble activity and stimuli would reflects how and if ensembles encode stimuli. As an example, ensemble spike distribution close from a random uniform spike distribution regarding stimuli location would carry no information. In most studies, such analysis does not simply exists.

Supplementary statistical feature extraction

The second criterion hinge on neuron ensemble coherence. Putting any kind of neurons in the same group can provide ensemble activities correlated to stimuli. However, it does not necessarily means that these neurons are similar and that gathering them in a group makes really sense. It is, thus, necessary, that a neuron ensemble is relevant regarding a similarity metric such that co-activity for non-overlapping method or synchronization levels for overlapping ones.

It is worth noticing that the Mölter and Goodhill's article provides a performance evaluation of clustering techniques on synthetic data based on the "Best Match" match clustering criterion and the number of predicted ensembles. This metric expresses how far are the estimated clusters from the ground-truth ones.

A formula that describes the group relevancy based jointly on the global level of synchronization and the statistical coupling between neuronal ensemble activity and stimulation could be at stake to assess neuronal ensemble algorithm efficiency.

Conclusion

We proposed a Bayesian statistical framework to cluster neuron in overlapping ensembles based on their synchronicity. The robust extraction of statistic parameter such as ensemble activity rate, ensemble activity recruitment rate and levels of synchronicity between neurons and ensembles given ensembles states provides an efficient and robust tool to go further in the accurate analysis of the correlation between neuron and stimuli in neuroscience. Our belief is that the mathematical effort to use Bayesian statistics to ensemble activation evaluation is the cornerstone to understand emergent neural network properties that embedded stimuli's information. To do so, we benchmarked as exhaustively as possible this algorithm, providing clear derivation of formulas in appendix and a free open-source software implementation to make this kind of statistical tool as easy to use and accessible as possible. This tool could allow neuroscientists to identify neurons having particular roles in the ignition of the ensembles. Furthermore, it could quantify the interaction of the neuronal ensembles with each other, measure the overlapping rate of neuronal ensembles and even measure the co-activation rate. These questions are of primary importance to many neurobiologists who generally seek to quantify the emergence of synchronized interactions between neurons. Possible improvements lie in the simplification of the model to allow analyses with a much higher number of communities or even to integrate more advanced methods of statistical temporal sampling in order to quantify the temporal evolution of synchronization changes between neurons and ensembles allowing a dynamic temporal clustering, thing totally missing for the time being. Finally, the use of estimated ensemble activity could finally reveal some sequential activation pattern that could be decisive in the understanding and prediction of behavior at living animals. 16

102

A Posterior Parameter Derivation

• α|Z, W, S, p, λ P (α|Z, W, S, p, λ)

α ∝ N i=1 A k=1 α z i,k k (1 -α k ) 1-z i,k A k=1 α α * k (α)-1 k (1 -α k ) β * k (α)-1 B(α * k (α), β * k (α)) ∝ A k=1 α N i=1 z i,k k (1 -α k ) N i=1 (1-z i,k ) α α * k (α)-1 k (1 -α k ) β * k (α)-1 ∝ A k=1 α N i=1 z i,k +α * k (α)-1 k (1 -α k ) N i=1 (1-z i,k )+β * k (α)-1 P (α l |Z, W, S, p, λ, α ∀k̸ =l ) α l ∝ α N i=1 z i,l +α * l (α)-1 l (1 -α l ) N i=1 (1-z i,l )+β * l (α)-1 α l |Z, W, S, p, λ, α ∀k̸ =l ∼ Beta(G l , Ḡl ) with G l = α * l (α) + N i=1 z i,l and Ḡl = β * l (α) + N i=1 (1 -z i,l ) • p|Z, W, S, α, λ P (p|Z, W, S, α, λ) p ∝ A k=1 T t=1 p w t,k k (1 -p k ) 1-w t,k × A k=1 p α * k (p)-1 k (1 -p k ) β * k (p)-1 B(α * k (p), β * k (p)) ∝ A k=1 p T t=1 w t,k k (1 -p k ) T t=1 (1-w t,k ) p α * k (p)-1 k (1 -p k ) β * k (p)-1 ∝ A k=1 p T t=1 w t,k +α * k (p)-1 k (1 -p k ) T t=1 (1-w t,k )+β * k (p)-1 P (p l |Z, W, S, α, λ, p ∀k̸ =l ) p l ∝ p T t=1 w t,l +α * l (p)-1 l (1 -p l ) T t=1 (1-w t,l )+β * l (p)-1 p l |Z, W, S, α, λ, p ∀k̸ =l ∼ Beta(H l , Hl ) with H l = α * l (p) + T t=1 w t,l and Hl = β * l (p) + T t=1 (1 -w t,l ) • λ|Z, W, S, α, p P (λ|Z, W, S, α, p) ∝ λ N i=1 T t=1 λ Gi (w t,Gi ) sit (1 -λ Gi (w t,Gi )) 1-sit × c∈C(A) z∈C(Z) λ c (z) α ( * ) (λc(z))-1 (1 -λ c (z)) β ( * ) (λc(z))-1 ∝ λ c∈C(A) T t=1 λ c (w t,Gi ) N i=1 sitδZ i =c × (1 -λ c (w t,Gi )) N i=1 (1-sit)δZ i =c × [. . . ] ∝ λ c∈C(A) z∈C(Z) λ c (z) T t=1 ( N i=1 sitδZ i =c )δZ i W t =z × (1 -λ c (z)) T t=1 ( N i=1 (1-sit)δZ i =c )δZ i W t =z × [. . . ] ∝ λ z c λ c (z) α ( * ) (λc(z))+ T t=1 ( N i=1 sitδZ i =c )δZ i W t =z -1 × (1 -λ c (z)) β ( * ) (λc(z))+ T t=1 ( N i=1 (1-sit)δZ i =c )δZ i W t =z -1 17
5. BINOE -BAYESIAN INFERENCE FRAMEWORK FOR NEURONAL OVERLAPPING ENSEMBLES λ z c |Z, W, S, α, p, λ z ′ ̸ =z c ′ ̸ =c ∼ Beta(T z c , T z c ) with T z c = α ( * ) (λ c (z))+ T t=1 ( N i=1 s it δ Zi=c )δ ZiWt=z and T z c = β ( * ) (λ c (z))+ T t=1 ( N i=1 (1-s it )δ Zi=c )δ ZiWt=z P (w tk = 1|Z, W -tk , S, θ) = P (w tk = 1|Z, W -tk , S, θ) ξ∈{0,1}
P (w tk = ξ|Z, W -tk , S, θ)

= 1 1 + P (Z,W tk =0,W -tk ,S|θ) P (Z,W tk =1,W -tk ,S|θ) = 1 1 + P (Z,W,S|θ)w t,k =0 P (Z,W,S|θ)w t,k =1 = 1 1 + ρ tk ρ tk = A l=1 l̸ =k T t ′ =1 t ′ ̸ =t p w l,t ′ l (1 -p l ) 1-w l,t ′ A l=1 l̸ =k T t ′ =1 t ′ ̸ =t p w l,t ′ l (1 -p l ) 1-w l,t ′ × p w k,t k (1 -p k ) 1-w k,t | w k,t =0 p w k,t k (1 -p k ) 1-w k,t | w k,t =1 × N i=1 T t=1 λ Gi (w t,Gi ) sit (1 -λ Gi (w t,Gi )) 1-sit | w k,t =0 N i=1 T t=1 λ Gi (w t,Gi ) sit (1 -λ Gi (w t,Gi )) 1-sit | w k,t =1 = 1 -p k p k × N i=1 i / ∈k T t ′ =1 t ′ ̸ =t λ Gi (w t,Gi ) sit (1 -λ Gi (w t,Gi )) 1-sit | w k,t =0 λ Gi (w t,Gi ) sit (1 -λ Gi (w t,Gi )) 1-sit | w k,t =1 × N i=1 i∈k (. . . )| w k,t =0 (. . . )| w k,t =1 ρ t,k = ( 1 p k -1) N i=1 i∈k λ Gi (       z 1 . . . z k = 0 . . . z Ki       ) sit (1 -λ Gi (       z 1 . . . z k = 0 . . . z Ki       )) 1-sit λ Gi (       z 1 . . . z k = 1 . . . z Ki       ) sit (1 -λ Gi (       z 1 . . . z k = 1 . . . z Ki       )) 1-sit P (Z ik = 1|Z-ik, W, S, θ) = 1 1 + η ik η ik = P (Z ik = 0|Z -ik , W, S, θ) P (Z ik = 1|Z -ik , W, S, θ) = P (Z, W, S|θ) Z ik =0 P (Z, W, S|θ) Z ik =1 = N j=1 j̸ =i A l=1 l̸ =k α z jl l (1 -α l ) 1-z jl × α z ik k (1 -α k ) 1-z ik Z ik =0 N j=1 j̸ =i A l=1 l̸ =k α z jl l (1 -α l ) 1-z jl × α z ik k (1 -α k ) 1-z ik Z ik =1 × N i=1 T t=1 λ Gi (w t,Gi ) sit (1 -λ Gi (w t,Gi )) 1-sit Z ik =0 i / ∈k λ Gi (w t,Gi ) sit (1 -λ Gi (w t,Gi )) 1-sit Z ik =1 i∈k = 1 -α k α k × N j=1 j̸ =i T t=1 (. . . )| Z ik =0 i / ∈k (. . . )| Z ik =1 i∈k × T t=1 (. . . )| Z ik =0 i / ∈k (. . . )| Z ik =1 i / ∈k 18 If i ∈ k, λ Gi (w t,Gi ) = P (s it = 1|w t,g1 = z 1 , . . . , w t,g k = z k , . . . , w t,gK i = z Ki ) If i / ∈ k, λ Gi (w t,Gi ) = P (s it = 1|w t,g1 = z 1 , . . . , ( ( ( ( ( w t,g k = z k , . . . , w t,gK i = z Ki )
The activity of the kth group is removed at each moment. In the activity dictionary this is equivalent to taking the kth component systematically equals to zero.

And so finally : Despite the development of advanced spike inference techniques (SIT) and the benchmarking of a few of them, the robust estimation of spikes from in vivo imaging remains difficult, especially in Hydra vulgaris. Indeed, fluorescent traces are noisy due to animal motions and deformations, artifact dynamics, "cross-talk" events or neuron tracking imperfections that largely contaminate the recordings. In this extent, asserting precisely which SIT is the most suited for Hydra vulgaris, demonstrating if a simple preprocessing step could account for properties not handled by the deconvolution techniques or the need to adapt or rethink the existing tools appears as an open question. Our first article [START_REF] Kubler | A Robust and Versatile Framework to Compare Spike Detection Methods in Calcium Imaging of Neuronal Activity[END_REF] created a versatile calcium fluorescent simulator with heterogeneous and more complex properties to benchmark the most promising SIT methods. The modeling of non-linear calcium dynamics, alternating between high-and low-frequency firing regimes, the variations of calcium dynamics or the different noise levels are assumed to provide further insights about SIT's generalizability to deformable free-moving animal case with burst activity events. The difficulties, conclusions and perspectives are detailed hereafter.

η ik = ( 1 α k -1) T t=1 λ Gi\k (w t,Gi\k ) sit (1 -λ Gi\k (w t,Gi\k )) 1-sit λ Gi (w t,Gi ) sit (1 -λ Gi (w t,Gi ))
Firstly, baseline is often neglected and therefore requires advanced pre-processing techniques. Indeed, most of SIT are applied to flat, stationary signals with no time-dependent baseline. This property is often encountered at mice visual cortex using a two-photon microscope. It is under these specific assumptions that the inverse problem of estimating spike locations is often performed. The damages on the estimates caused by the violation of these assumptions has been studied, making the choice of a method to deconvolve Hydra vulgaris' data elusive. Photobleaching is another phenomenon that requires to detrend the signal [START_REF] Nolan | Algorithms for the Correction of Photobleaching[END_REF] making it inherently non-stationary. Indeed, fluorophores switch to dark states leading to a global decrease of fluorescent intensity. Two methods have been tested: 1) Applying a detrending step before the SIT ; 2) Performing no detrending step before the SIT. The tested detrended techniques were model-free. On the data, interpolation (Fixed Degree Polynomial and Cubic Spline interpolation), convolution (Sliding window, wavelet transformation) and filtering (Frequency filtering in the Fourier domain) have been applied. Overall, the spike estimation is better with such detrending steps (frequency filtering and spline interpolation yielding the best results) leading to satisfying results in many cases (see Figure V.1). By doing so, larger dynamics are filtered out before the deconvolution step.

Deconvolution of neuron 400

Increasing Spike Amplitude Figure V.1: Deconvolution of a low noise signal : Neuron 400. The signal is detrended using a sliding window of size 20 and denoise using a gaussian filtering. Spike Amplitude variation is highlighted.

If the spike frequency is too close from the underlying baseline or noise frequency, meaningful information can be filtered out during the detrending step and not recovered afterwards. The level of noise seems to largely vary from one fluorescent trace to another making some data completely impossible to interpret in the absence of ground-truth annotations. This is also true for underlying non-linearities extremely heterogeneous between fluorescence traces.

If the user-defined parameters are not properly tuned for detrending, artifacts appear and cause spurious dynamics. For example, too high a degree for a polynomial interpolation, a wrong choice for the cut-off frequency for frequency filtering or too long a size window for the convolution methods can generate wrong low-frequency dynamics. Furthermore, it seems that a single user-defined parameter set does not exist for all data, especially when fluorescent traces are extremely heterogeneous as it is the case for Hydra vulgaris. This phenomenon can be directly observed on Hydra's deconvolved matrix [START_REF] Kubler | A Robust and Versatile Framework to Compare Spike Detection Methods in Calcium Imaging of Neuronal Activity[END_REF] where horizontal lines suggest that some neuron deconvolution completely failed while others with regular pattern seem to correspond correctly to Rhythmic Potential or Contraction burst activities already highlighted in the animal [2] [START_REF] Lagache | Tracking calcium dynamics from individual neurons in behaving animals[END_REF]. An adaptive way to estimate the user-defined parameters that allow a correct detrending for each individual fluorescent signal combined with a more efficient upstream correction of pollution sources would be at nice extension of this work.

Inhomogeneities of the dataset induce a strong variation of individual neuron signal characteristics that raises a trade-off between accuracy and robustness. Indeed, on synthetic data, the most accurate technique was one handling a time-varying baseline. Indeed, MLspike from [START_REF] Deneux | Accurate spike estimation from noisy calcium signals for ultrafast three-dimensional imaging of large neuronal populations in vivo[END_REF] was proved to be the most accurate method on synthetic data according to the proposed simulator. This method is extremely powerful on a single fluorescent trace when user-defined parameters have been carefully tuned for [START_REF] Kubler | A Robust and Versatile Framework to Compare Spike Detection Methods in Calcium Imaging of Neuronal Activity[END_REF]. However, once we change the fluorescent trace properties, parameters are not well-suited anymore, and accuracy collapses quickly making this method the worse on a Hydra vulgaris dataset. It highlights the trade-off between accuracy and robustness to choose a SIT. The more homogeneous the dataset is, the more similar the method parameters will be. Conversely, if data are ex-

Hard CB/RP

No spike RP CB Figure V.2: Variety of neural activity with different dynamics. Fluorescence traces show simultaneously a neural response to CB events and RP events, observation inconsistent with the non-overlapping nature of the neural networks highlighted by [2] and corresponding to cross-activity contamination. Some traces do not show spiking dynamic at all while others highlight spiking dynamics completely hidden in the noise. tremely heterogeneous like for Hydra, parameters will not be adapted to the entire dataset and a more complex method would fail. This paradox is demonstrated qualitatively by the regularity of the inferred spiking matrix got from a simple naive SIT based on smoothed gradient calculations on experimental data ( [220] Fig. 5.a) and quantitatively using the Global Jaccard Index on the deconvolved spiking matrix ( [220] -Fig. 5.b)) finding a higher degree of inter-spike sequence correlation.

An adaptive estimation of user-defined parameters for each fluorescence trace could be insufficient to deconvolve data. Indeed, surprisingly, some obtained fluorescent traces positively respond to both RP and CB events which is not relevant with previous biological observations made at Hydra vulgaris in [2] claiming that its neural networks are non-overlapping (See Figure I.4). A neuron should only response positively to a single specific network at a time. Two main factors could account for this problem: Burst Pollution and Tracking Errors. Burst pollution is linked to the fact that during a burst event Hydra contracts rapidly. Lots of cells are, then, gathered in a cluttered environment and the global intensity activity seems to increase in every cell, even RP neurons. Spurious information, thus, emerges. It means that, experimentally, about 50 percent or more of the detections could be lost at each fast and strong contraction of the animal and be wrongly reassigned to another cell in the tracklets. In that situation, the beginning of a fluorescence recording, and its end could belong to two completely different somas with different time dynamics.

The Poisson homogeneous assumption of spike inference techniques combined with a fixed calcium transient cannot handle firing switching regimes. Indeed, if we seek to apply a SIT on a fluorescence trace where calcium transients can change according to firing switching regimes, the method is not adapted since the kernel is often estimated as fixed through, for instance, auto-regressive coefficients and artificial time-lags can be induced in the inferred spike locations. An adaptive dynamic deconvolution kernel that could switch from one firing regime to another would be a key solution to handle such time-varying dynamics. It is worth noting that in Hydra vulgaris, the coexistence of several dynamics are experimentally observed in the data (see Figure V.2) even if it seems inconsistent with the literature. Some experiments have also been carried out without detrending preprocessing step before running SIT. In this situation where the baseline is not handled anymore, each positive variation of the signal is assumed to be a meaningful information by the algorithm, and it locates spikes basically everywhere. However, the output spike amplitudes vary according to a level of confidence. The higher the spike amplitude is, the more confident the spike existence at that location will be. The extraction of the highest amplitude spikes could be an attempt to find the meaningful information. The fundamental question is about thresholding the spike amplitudes given that no specific distribution seems to emerge and that a hard threshold that works for all traces of the dataset seems not to exist. Finding similar neuron traces and jointly extracting their DAT by a time-point to time-point statistical coupling robust to false detection was an attempt to fulfill that task. It assumed though that noise would be uncorrelated, which is not necessarily true in practice since noise is partly brought by a motion that could be shared between several neighboring cells.

As a conclusion, the first part of this thesis was aiming to find and justify the use of an existing spike estimation technique to apply on Hydra vulgaris dataset. Even if SIT are often very advanced algorithms applied on low noise data, a gold-standard method has not been proved yet. The fundamental heterogeneous nature of the data unveils that Hydra vulgaris dataset are far more complex compared to other animal models (mice/zebrafish) due to a lot of additive sources of pollution and a perfect SIT adapted for Hydra seems not to exist yet. The track that has been followed, then, was to firstly apply one the most promising technique in the benchmarking CDfoopsi [START_REF] Pnevmatikakis | Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data[END_REF] with a frequency filtering detrending as preprocessing step ; and then, knowing that neural spiking sequences were corrupted by errors (False Detections/ Miss Detections / Time-Lags), to implement a functional connectivity metric based on point-process theory to correct for the errors and estimate the network topology. Two fundamental prospects need to be highlighted. Firstly, the noise cannot be properly removed by just applying a model-free denoising technique. It seems that estimating nonlinearities via a model is at stake. The estimation of the field of view motion of the animal, the z-axis motion of the animal, the correction of neuropile contamination and the filtering of global contaminating bursting activity should be considered. Some solutions emerge. To correct motion, deformation of the animal and improved the soma tracking: A model of movement should be embedded. Kalman filter is one possibility but multi-spectral recording via a second nuclear calcium indicator called TdTomato, that does not vary with neural activity, is currently under the scope. This other solution records the fluorescent intensity of nuclei to efficiently estimate the motion and deformation of the animal in the field of view or perpendicularly to the z-axis. Secondly, the use of an independent component analysis on cell fluorescent intensity could reveal common components. Applied to neighboring cells, it could estimate the z-motion of parts of the animal since neuron cells in tentacles, for example, might have the same intensity decrease if an animal's tentacle leaves the focal plane at a given moment. Applied globally, it could estimate the global intensity increase observed during burst events and be used to correct the pollution above all neuron activities. Neuropile contamination could be reduced by averaging intensity outside soma in its neighboring area. An active contour applied to segment Hydra, combined with Gaussian kernel fitting on soma's neighboring could seem relevant to fulfill that task. Another perspective does not concern the denoising but the spike inference technique itself. Implementing a SIT that could embed a time-varying baseline with an adaptive kernel that switches statistically between several firing regimes would be promising to get rid of the Poisson Homogeneous spiking distribution hypothesis. At the time these observations were made in 2020, the use of Auto-Regressive Markov Switching Model [START_REF] Ailliot | Markov-switching autoregressive models for wind time series[END_REF] was considered. It is worth noting that, in the meantime, a promising article was published highlighting these shortcomings and attempting to provide solutions [START_REF] Diana | High frequency spike inference with particle Gibbs sampling[END_REF]. Finally, the SIT usually works on a single fluorescent trace at a time, but we expect these fluorescence traces to be mutually correlated. A simultaneous multi-fluorescence trace deconvolution would seem relevant. It means that information would need to be exchanged during the deconvolution to increase the confidence in the spike locations regarding several traces in the same time. In other words, if a spike was infered at the same time in several fluorescent traces it would make it more likely to exist.

Time [F]

Fluorescence

Correction of the CB pollution

Correction of the neuropil contamination

Correction of the z-axis motion

Statistical coupling through Spatial Point Process Theory

Functional connectivity tries to infer the co-activity coupling between individual neuronal cells. It allows to create a connectivity map that demonstrates the existence of functional neural networks with specific topological structures. The coupling metric is data-dependent. If neural signals are time spiking sequences corresponding to binary sequences of action potential emissions, point process theory is more adapted than regular correlation metrics for time-series. We highlighted that spiking sequences were polluted by spike false detections, miss-detections and time-delays due to additive sources of noise not properly handled. Implementing a functional connectivity metric based on a spatial statistics framework able to estimate robustly such coupling was at stake in the second publication of this thesis [START_REF] Kubler | Statistical Coupling Between time Point-Processes[END_REF].

i Using spatial point processes for characterizing the temporal coupling of biological processes

Given a point process density, the number of expected points, randomly located in any subspace area, is known to be directly proportional to its volume. Such statistical consideration is known as Complete Spatial Randomness hypothesis. It allows to analytically derive a H0 null hypothesis of the expected number of points in any subspace got by chance. Comparing the observed number of points in subspaces against the expected one under CSR, via a K-Ripley difference vector [START_REF] Ripley | The second-order analysis of stationary point processes[END_REF], demonstrates aggregation or dispersion effects corresponding respectively to positive or negative statistical coupling in spatial point location distribution. When the value is around 0, independence can be assumed. Such analysis have been developed for biological applications to quantify immune response for instance [247] [248].

a Connectivity between individual neurons

Our second article proposed for the first time an adaptation of this framework in a 1dimensional case between time point-processes to infer the functional connectivity between neuron activities being robust to false detections, miss-detections and time-delays. The basic idea is to measure the statistical coupling between a first time point sequence relative to a reference second one. If enough points aggregate systematically at the same distance from the reference points, a distant statistical coupling can be highlighted. This distance corresponds to time delay in that context. By averaging the statistical coupling got for different time-delays, we can recover the global statistical coupling between the two date sequences (see Figure V.5).

On synthetic data, the method demonstrates its ability to cope with false detection events. Indeed, if noise adds randomly false positives in neural data binary sequences, the algorithm will estimate that these points correspond to a random activity. Only an aggregation of false detection points systematically located with the same time lag could bias results by being estimated as positive coupling. The method appears to be less robust to miss-detections since pair of coupled spikes cannot be recovered anymore when information is missing in the recording. Concretely, we can see that results are largely damaged by miss-detections given that for a 50% miss-detection rate, more than 80% percent couples are not matched.

The strength of our statistical coupling framework in comparison with more general frameworks like Pearson correlation or Jaccard index for example is its intrinsic ability to detect the time delay of the statistical coupling. Especially, if neuronal spiking activity has been deconvolved using fixed kernels, time-shifts in the spike locations can occur and bias the classical functional correlations.

b Coupling between stimuli and ensemble activation

Importantly, the provided statistical tool from [START_REF] Kubler | Statistical Coupling Between time Point-Processes[END_REF] can be adapted to achieve other critical neuroscience tasks such at the neuronal statistical coupling to stimuli.

Estimating neural ensembles means that neurons are defined by common synchronous activity patterns. In the neural doctrine, these patterns are the intrinsic representations of a cognitive task, the integration of a stimulation or the activation of micro-circuitry in response to a behavior and would represent the elementary bricks of the neural code.

It means that, first, estimating neural ensembles using a community detection technique necessarily needs to extract ensemble activities. Then, these activation patterns need to be correlated with stimuli submission. Obtaining a neuron clustering where neuronal ensemble activity are not related to the stimulations would be extremely complicated to analyze in terms of consistency and to compare with other clustering techniques. The best neural ensemble detection technique would be a method such that: 1) it makes sense to locate consistent neurons in the same ensemble regarding a specific criterion or metric in relationship to a group definition ; 2) The inferred neural ensembles have an activity that does "explain" the submitted stimuli or the observed spontaneous behaviors.

In the same way that GSODA [START_REF] Mukherjee | Generalizing the Statistical Analysis of Objects' Spatial Coupling in Bioimaging[END_REF] was the adaptation of SODA [START_REF] Lagache | Mapping molecular assemblies with fluorescence microscopy and object-based spatial statistics[END_REF] by transforming the statistical point-to-point coupling into a point-to-object coupling problem, it is possible to modify the developed statistical framework to infer the correlation between a neural activity and a set of stimulations. The research spaces, called Cluster Cores, initially corresponding to concentric rings potentially overlapping and centered around points in a reference sequence, can be matched with the time-space domain corresponding to the submission of stimuli or the execution of a spontaneous behavior. The coupling between a set of points and a time-domain subspace, can thus, be performed to highlight the correlation between neural ensemble activities and the triggering event (see Figure V.5).

ii Using point processes to calibrate SIT

When no detrending step is applied on calcium fluorescence data, the output spikes have different amplitudes. The spike amplitude seems correlated with its probability to occur. At a given point, when spike amplitude decreases, the resulting spikes just correspond to noise. Let's suppose that we can estimate two coupled neuron fluorescence traces with such noisy spike sequences. Finding pairs of coupled neurons could use the statistical coupling framework to find the optimal thresholds DAT (1) and DAT (2) able to separate meaningful spikes from spikes associated with noise. The search-grid research, using the statistical coupling metric between spiking sequences demonstrates a 2-dimensional curve-bell shape whose maximum value, seems to provide an appropriate estimated threshold to differentiate meaningful information from noise components. This analysis is under a writing and submission process in the article "Statistical calibration of deconvolution methods for extracting neuronal spikes from calcium imaging", Kubler et al, not published, and is equally supported by the same mathematical framework.

The method questions equally the correlation of the noise. If the noise is uncorrelated, the framework seems well-suited since we expect to have false spike detections equitably spread over the sequence length that can be filtered out by the statistical test. However, the non-linear dynamics observed in fluorescence traces could disturb the signals at same time locations. This could be relevant if different neurons undergo the same motion at the same time for instance in a free-behaving animal. In this situation, false spike detection could appear at the same time locations in different fluorescent traces creating a spurious aggregation of points and a residual artificial coupling. Such situations could illustrate the limits of the proposed generic algorithm and going backward to the correction of additive noise sources would be at stake.

Statistical inference for overlapping neural ensemble detection

Like previously highlighted in Section I.3.i, a formal consensual definition of neural ensemble does not exist. The definition, explicitly expressed by [START_REF] Carrillo-Reid | What Is a Neuronal Ensemble?[END_REF] as "a group of neuron repeatedly firing together", is adopted for the fourth article of this thesis under writing and submission process of "Kubler et al, Bayesian inference of overlapping neuronal ensembles, Nature Computational Sciences, under-review".

The idea is to provide a statistical framework able to cluster neurons in overlapping ensembles based on their synchronicity. Estimating biological-inspired statistical evaluators that could identify specific neuron roles, extract ensemble neural activity, measure the clusters relevancy and correlation with stimuli would be a novelty that could robustly improved in vivo neural network analysis. This work is an extension of the work initiated by [START_REF] Diana | Bayesian inference of neuronal assemblies[END_REF] but handles overlapping and provides further visualizations and considerations.

Synchronicity is the cornerstone of such statistical analysis. Its basic principle is that, at each time step, a neural population activates, corresponding to all neurons that emitted at least one spike during this time recording window. When a neuron population is repeatedly reactivated over time and during a cognitive task, it is likely that this population constitutes a neuronal ensemble that encodes a mental function. These activations are intrinsically characterized by randomness since the exact same group is never perfectly recalled in practice. When a group activates, only a subset of its neurons can contribute to its activity, but other neurons, which are not a part of it, can equally spontaneously fire meanwhile. This inherent stochastic nature of the problem makes the use of statistical inference a promising solution to consider estimating emergent firing properties.

Overlapping is a clustering feature that makes this open-problem even more complex. It corresponds to neuron capability to be recalled by several different neural ensembles simultaneously or at deferred moments.

The provided algorithm has been applied, firstly, on much simpler model animals than Hydra vulgaris like zebrafish data from [3] and mice visual cortex data from Yuste's laboratory. The aim was to prove that the algorithm works on gold-standard model animals before applying it to much noisier data.

i Synchronicity unveils modular and hierarchical organization of neuronal ensembles

Modularity corresponds to the existence of neuronal populations repeatedly reactivated over time. If a neuronal population is sufficiently distant from other repeating neuronal populations, it will form a module of interconnected neurons. Each neuron will be defined by its rate of synchronization to the group. For a synchronized module to emerge, it is essential that its activation rate is much higher than the spontaneous basal activation of neurons, and that the synchronization rate of its neurons is much higher than their de-synchronization rate. It is also necessary that the number of neurons making up the module be sufficiently large for a group activity to emerge and not be drowned out by activation noise.

Detecting overlapping neuronal ensembles based on synchronicity has to integrate the following principle: statistical significant silences also convey information. Let a neuronal ensemble encodes a cognitive function. If a subset of neurons in this group is silent a significant number of times (= desynchronization effect), this might mean that two regimes of synchronization can emerge in this set. Within the same group, the distribution of neuronal synchronizations could be multimodal, corresponding to several layers of hierarchy.

The most synchronized neurons during ensemble activation can play a specific role for this ensemble. They could, for example, represent the "Completion Pattern Neuron" (CPN) previously highlighted in the literature. The most desynchronized neurons are the closest neurons from neurons outside the group whose activity could be completely random for instance. Overall, if sometimes a silent subgroup of neurons also exhibits external coordinated activations, i.e. out synchronizations from its main ensemble, this might correspond to neurons involved in another behavior or cognitive function. In this case, these neurons may be located at the intersection of several groups, forming a functional overlapping. An a-synchronization, corresponding to a synchronized activation of a neuron outside its group, could therefore be the major pillar explaining the overlapping phenomenon based on synchronicity preventing sometimes two distinct groups to merge.

It is worth mentioning that on experimental data, fluorescence trace intensity of neurons that respond to a stimulation increases. Sometimes, spikes are even detected but the spiking neuronal response can be delayed or can last different durations. In this situation, the time dynamics of the response is not handled by the statistical model and artificial desynchronization events can even occur and generate artificial overlappings. Without accounting for the time dynamic spiking response of the neurons, in the statistical model, it is fundamental to preprocess data to avoid this effect. It is the neuronal desynchronizations between neuronal wave responses that should provide the overlapping effects in our analysis and not the intra-wave desynchronizations (see Contribution 4 -Fig. 7).

ii Evaluating the ensemble detection algorithm

The main advantage of this framework is its ability to estimate many parameters supported by statistical inference that are biologically meaningful related to synchronicity. The size of the group, their activation rate, their average level of synchronization, the basal neuron activation, the levels of synchronization, the (de-)synchronization rates or the conditional firing statistics given membership provide promising indicators to check that the output clustering makes sense and can be directly related to biological features.

The adaptation of the statistical coupling tool introduced in [START_REF] Kubler | Statistical Coupling Between time Point-Processes[END_REF], provides another way to quantify the relevance of the clustering results. Obtaining neuron groups whose activities do not explain stimuli or behaviors would be quite hard to analyze or interpret. Obtaining a correct neuronal ensemble detection algorithm means having consistent synchronized neurons gathered within ensembles whose activity efficiently explains the cognitive tasks. Deriving a quantitative metric to do so, is then, straightforward.

Using this framework, we can represent schemes never highlighted before, to our knowledge. Firstly, the synchronization diagram allows to represent at the same time every combination of overlapping membership, their relative proportion, and the relative importance of each neuron within its ensemble. This information is provided by the level of synchronization corresponding to an activation synchronization rate weighted by the relative contribution for the neuron itself in its own activation. This diagram directly targets completion pattern neurons providing a tool to identify ideal neurons to be manipulated using optogenetics.

Another helpful representation to illustrate the ensemble detection is to directly map neuronal ensembles in the in vivo neural networks of model animals (see Figure V.7). By doing so, it is possible to track over frames, the neuronal activations and visually explain why a neuron has been activated at a specific time frame and whether this activation has been triggered by an ensemble or just corresponds to spontaneous neural activity.

iii Convergence of Bayesian inference

The exploration of the random parameter space is performed using a Bayesian Inference algorithm called Gibbs Sampling that iteratively samples each variable assuming all others to be fixed from the posterior conditional distribution. In a high-dimensional space with a lot of variables, a problem occurs : the convergence toward local maxima. To solve this issue, we implemented diverse statistical techniques to systematically converge toward the correct global maximum during the statistical optimization process.

The successive random initialization strategy was the first appealing technique implemented (see Figure V.8). Largely used in machine learning with clustering algorithms like K-means for its simpleness, it aims at initializing the optimization routine at different random locations to get closer to the global maximum and converge at least once toward the correct solution. The results were promising on the simulator since, for simulation parameters relevant for Hydra's statistics, in a dozen of random initializations the global maximum was systematically reached. More importantly, it demonstrated that the correct clustering result was always reached for a maximum value of likelihood that, thus, could be used as a cost function (see Figure V.8). Unfortunately, these results failed on experimental data like mice brain visual cortex probably due to a higher data complexity.

Another implemented technique was to randomly initialize variables and run the Metropolis-Hasting algorithm embedding simulated annealing (see Figure V.9). The general idea of simulated annealing is to embed noise in the sampling to voluntary move the samples. By allowing samples to be shaken at each iteration you obtain a non-null probability to get out from a local maximum. The MH algorithm looks at whether the sampling deviation improves or worsens the result in terms of likelihood ratio and accepts or rejects it accordingly. The temperature corresponding to the probability to shift the sample decreases over the iterations to avoid final oscillations and allow for convergence. This idea was implemented by adding random or targeted neuron subsets mixing each time a log-likelihood convergence regime was reached, the mixing proportion decreasing over the iterations, like the temperature of the simulating annealing process.

This solution applied alone did not produce convincing results. A possible interpretation is that when initialization is random, the groups at the first iteration are random and group activity is not necessarily consistent. If we reach a local maxima regime at that moment, ensemble activity could have any sense. Trying to group certain subsets of neurons into nonsense groups may fail, as the coherence of a group may not have been sufficiently achieved. In other words, putting non-consistent neurons in the same group provides a meaningless group whose activity that averages the activity of its dissimilar neurons provides a dissimilar clustering. This tool may, though, be used with the following technique.

Another strategy to avoid local maxima is a smart initialization inspired from spectral techniques. The smart initialization avoids performing a random initialization of the clusters at iteration 0. The idea is to find directly a smart initial guess either of the neuron clustering or of the neuronal ensemble activity. To do so, a technique inspired by the work from [START_REF] Carrillo-Reid | Endogenous Sequential Cortical Activity Evoked by Visual Stimuli[END_REF] [START_REF] Carrillo-Reid | Identification and Targeting of Cortical Ensembles[END_REF] has been carried out. If we see each time bin as a population vector, a first initial guess is to estimate heterogeneous repeated population vectors. By performing a spherical K-means algorithm to cluster the statistical significant time steps, we can unveil such populations and initialize the ensemble activity matrix. The initial clustering step used to estimate ensemble activity at iteration 0 was performed using K-Means with cosine similarity on statistical significant population vectors extracted by a z-score test using the framework implemented in [START_REF] Dhillon | Concept Decompositions for Large Sparse Text Data Using Clustering[END_REF] [251].

More advanced sampling statistical algorithms could be considered to avoid local maxima jointly with a smart initialization. Indeed, some attempts have been implemented but were not providing convincing results using more advanced sampling routines like Collapsed Gibbs sampler. The idea is to sample variables using marginalized posterior conditional distributions to improve the speed of convergence [START_REF] Van Dyk | Partially Collapsed Gibbs Samplers: Theory and Methods[END_REF] [253] [START_REF] Robert | The Bayesian Choice: From Decision-Theoretic Foundations to Computational Implementation[END_REF]. The marginalized distributions are derivable in our case and provide a framework more similar to [START_REF] Diana | Bayesian inference of neuronal assemblies[END_REF]. However, to our opinion the issue of convergence was not related to speed but location. In only few iterations, a convergence regime was systematically reached (good or bad) but waiting a huge number of iterations did not change the convergence location meaning that the problem was not related to speed but to local maxima according to the experimental simulations. Results were not convincing. It could, however, appear as an improvement to check whether it increases the accuracy or not.

In brief, smart initialization demonstrated its ability to avoid local maxima on simulated data. On experimental data, in the absence of ground-truth, it is not possible to assert the same. Results are promising though and the combination of smart initialization, simulated annealing and more efficient sampling routine could provide even better experimental estimates.

iv Estimating the number of communities

An implementation of the algorithm has been carried out to automatically estimate the number of communities. This task has been performed by dynamically increasing the number of communities A to scan the model with a varying number of neuronal ensembles. Stopping criteria were assumed to be able to detect the correct number of communities. Indeed, on simulated data the framework demonstrates its ability to uncover the exact number of neuronal ensembles using the loglikelihood and the statistical evaluator Akaike Criterion (AIC) as stopping criteria. Loglikehood measures the fitting quality between the statistical model and the observation data. AIC does the same but penalizing the dimensionality of the model. Experimental results on synthetic data demonstrated the ability of both criteria to respectively provide maximum and minimum values for a correct number of ensembles. This number of ensembles is changed dynamically each time a convergence regime is reached. Such analysis can be quite long and the algorithm is, as implemented, under effective for too high a number of ensembles (>8). Indeed, the formulated statistical model is mathematically exact and tries to infer all statistical synchronization probabilities for each possible combinations of ensemble activation given each possible combinations of membership. This computation allows to exhaustively characterize the statistical interactions between ensembles via statistical synchronization excitatory and inhibitory effects, but appears as computationally very heavy.

The model could, thus, be simplified, since we can legitimately estimate that the probability of having more than 2 co-active communities is low in practice for simple real-world experiments. Simplifying the model of interaction between neurons and ensembles could, therefore, be a promising perspective.

The automatic detection of the number of neuronal ensembles failed on real experimental data. Indeed, it seemed that the higher was the number of communities, the better were the evaluation criteria. This experimental observation made on the in vivo real data seems to correspond to an overfitting phenomenon related to a higher level of complexity of realworld data compared to synthetic data. In this situation, the higher is the number of communities in the model, the more relevant seem to be the synchronization patterns since outliers gathered in non-relevant groups making the remaining ones more consistent.

A promising track to infer the correct number of communities on experimental data is to create quality statistical evaluators that summarize two quantities : the global neuron to ensemble synchronization level and the statistical coupling between neuronal ensembles and stimulations. The first corresponds to an evaluation metric that describes the consistency of the providing clustering in terms of synchronization. The idea is to measure the synchronization rate of all neurons to their respective ensembles by weighting their synchronization with their own synchronization rate. The second is a metric that accounts for how explicable the neuronal ensembles are about stimuli so as to avoid adding non-interpretable ensembles in the model regarding their activities. In Hydra vulgaris, stimulation does not necessarily exist and the statistical coupling should be carried out automatically with detected spontaneous behaviors.

v Perspectives

To summarize, statistical inference based on biology inspired models are promising techniques to analyze emergent synchronization patterns from neuronal individual cells. It largely demonstrates the ability of such frameworks to estimate further meaningful and interpretable statistical features of neural micro-circuitry that pave the way to a better understanding of the neural code. Various animal models demonstrated the integration of neural information through ensembles whose co-activity explain how in vivo neural networks embed a stimulation, a cognitive task, or a spontaneous behavior via a synchronized neural representation. However, this promising framework could be extended further. Time dynamics is currently missing in the statistical model and could be added as a progressive temporal neuron activation. This corresponds to more complex but more biologically relevant activation patterns. Instead of seeing ensemble activation as an instantaneous event where a significant statistical subset of neurons is co-recruited, we can see the phenomenon as a progressive neuron activation and extinction in an avalanche context (see Appendix A). Each ensemble could be characterized by activations spread according to a Poisson distribution to model the occurrence of events that will trigger each time a progressive activation and extinction of neurons during a firing event whose duration could be modeled through a Gaussian distribution. By doing so, the existing normalization preprocessing steps applied on experimental data could be removed, and more accurate results be obtained. Such a model of temporal activation has already been implemented as a generative simulator in this thesis, but raises fundamental questions about statistical estimation. Complexifying the statistical model means increasing the model dimensionality. Supplementary statistical parameters and hidden variables should be estimated while the number of ob-servations would remain the same. Local maxima issue will necessarily emerge and will require advanced techniques too bypass it. Questions would be raised about how to embed external estimation and initializations to avoid such issues, how to implement more powerful sampling techniques, add some regularization in the iterative process and be finally certain at the end of the day that the algorithm converges correctly and systematically toward the good estimate in a reasonable amount of time without ground-truth.

In Zebrafish a very surprising observation has been made. If the time window on which the algorithm is run on is too long, consistency of the results decreases (> 1000-time bins). This seems a paradox since on synthetic data, longer was the simulation duration, higher was the result accuracy. This observation raises another fundamental question: Are the neuronal co-activation patterns the same over time or is there a time evolution (a shift) of the synchronization patterns and consequently of the clustering itself? This phenomenon could be explained by photobleaching and, thus, the change in the neuronal ensemble response. It hence seems that targeting a time window to run the proposed algorithm where consistent synchronous firing patterns are present is necessary. Otherwise, we could attempt to average inconsistent statistics and results could provide spurious information. It raises the fundamental question of the time scale definition of a neuronal ensemble.

The neuronal ensemble could, thus, be estimated as a time-varying phenomenon. The clustering matrix and all other statistical parameters of the proposed framework are assumed to be fixed in the algorithm. However, they might also be represented as functions of time changing the sampling routine. It is fundamental to mention that these questions also appear here as a promising perspective to go further in such Bayesian analysis. Discussions with neurobiologists confirm this need since in in vitro and in vivo neural networks. Evolution is fundamental to understand neural network properties through time dynamics (Appendix A). Having a statistical inference tool that demonstrates explicitly how the synchronization patterns, the expected number of ensembles, the modular and hierarchical structure of neural networks evolve over time, would be, a breakthrough in neurosciences. In Bayesian inference, some advanced statistics have been implemented to analyze time-varying statistics. The general idea is to sample variables using time-dependent conditional posteriors on previous samples. One of the most relevant tool that could be used is called Particle Sampler (PS). It already demonstrates its ability to infer spike locations as a SIT [START_REF] Diana | High frequency spike inference with particle Gibbs sampling[END_REF] and its use for neuron ensemble detection would be a promising perspective that, to our knowledge, has never been explored yet.

Chapter VI

Conclusion

Thesis summary

This thesis work has addressed the different steps from calcium imagery to neuronal ensembles. It, firstly, benchmarked spike inference algorithms to estimate neural activity. Our work highlights that a gold-standard technique does not exist, and the choice of the most adapted one is data dependent. Furthermore, the noisy fluorescence traces extracted, point out the limitations of the existing techniques directly applied to a free-behaving animal. Uncorrected motions and deformations, neuropil contamination, photobleaching effects, tracking errors, non-homogeneous spiking activity and inter-pollution of micro-circuitry make the algorithms insufficient to robustly extract neural activity. In the absence of gold-standard methods, a functional analysis of in vivo neural networks needs to be supported by statistical inference and to account for residual errors in the neural activities. Secondly, this work demonstrates that point process statistical coupling can be used to infer topology and correct this spurious information. Uncorrelated errors can be, thus, compensated to provide a neural network connectivity matrix but its ability in neuroscience goes beyond this simple application. Its use is, equally, promising to measure the statistical coupling between neural ensemble activity and stimuli and could be directly embedded in the definition of what a neuronal ensemble is. This part equally demonstrates that a functional connectivity metric projects observations in a space where data are hardly interpretable biologically speaking and user-defined choices make the clear understanding of neural networks and method comparison even more complicated. The final part of this work invites the reader to question the definition of neuronal ensembles. The more recent neural doctrine accounts for the neural representation as synchronization patterns emerging during in vivo neural networks development. It seems that neuronal ensembles of co-active neuron cells embed a representation of the world, the cognitive computation, or the behaviors through co-activation patterns which involve neurons with specific roles and whose optogenetic manipulation can alter microcircuitry functioning and, consequently, mental faculties, behavior, and decision-making. Based on this neural doctrine, our work finally provides a statistical inference framework to estimate at the same time the neuronal ensemble code, the neuron clustering into functional ensemble and neurons with specific roles. This method, supported by statistical evidence, adopts a biologically inspired definition of a neuronal ensemble making the understanding of the inferred quantities easily interpretable. It is, to our knowledge, the first time, that an overlapping neuron clustering algorithm based on Bayesian Inference provides such analysis and discusses the output clustering relevancy by embedding a statistical coupling framework with the applied stimulations.

Discussion of limitations and prospects

The recent progress of calcium fluorescence imagery techniques, genetical engineering and tracking algorithms of intermittently visible objects makes now possible the functional analysis of neurons in in vivo networks, at the scale of a whole animal like Hydra vulgaris. The use of such simple model animal from an aggregation of cells to an entire free-behaving organism allows for the first time to image neural networks birth, evolution, stabilization, and plasticity. During the course of this thesis, the implementation of a robust fully statistical pipeline predicting Hydra vulgaris behaviors could not be carried out. However, each step of the pipeline has been studied, providing efficient tools, promising results, and further prospects. It is worth noting that each step requires the implementation of tools since lots of existing techniques are too basic to be directly applied on a free-moving animal case such as Hydra which is highly polluted by noise. Three major axes are, according to this work, promising tracks to provide evidence of the neural doctrine in Hydra.

Improving the denoising step by handling the different sources of noise using an explicit model would be fundamental. Indeed, we have finally highlighted that a robust denoising step for free-moving animal does not exist yet. Generally, this step is directly performed on the movies, using convex optimization requiring a very advanced mathematical background. Such adaptation to the free-moving animal case could sufficiently improve the data by its own. It would allow our statistical Bayesian framework to work directly on Hydra like it already does on mice and zebrafish visual cortex.

The neural doctrine claims that neuronal ensembles encode more robustly information than individual cells. Hence, comparing the level of accuracy in the prediction of a LSTM, by inputting the ensemble activities instead of neuron cell activities or pixel intensities could provide further quantitative insights that this neural doctrine is correct.

Finally, Bayesian Inference largely applied in statistical physics has proved its ability to unveil emergent dynamics from individual cells. Complexifying our framework to estimate the avalanche dynamics jointly with the neuron clustering would be a great improvement. A breakthrough would be to provide a framework where neuronal ensemble is related to synchronous activation patterns that evolve as a function of time. Such an evolving definition of a neural ensemble would be promising for neurobiologists to quantify organization emergence in in vivo neural networks through the extraction of time-varying statistics. The combined use of such statistical frameworks with optogenetics to stimulate, disturb or kill neurons in cultures could provide further insights about how neural network embeds information and how we can interact with it. The alteration of neurons with specific roles could even have an impact on behaviors.

Putting it into perspective

The recent interdisciplinary works reviewed in the manuscript highlight that research on functional analysis of neuronal connectivity speeds up. An exponential increase of the publications dealing with functional connectivity has been observed in relation to the emergence of social networks. Bottom-up approaches that try to uncover neural doctrine at the scale of in vitro and in vivo neural networks exploded over the past few years by applying ever more advanced statistical methods. Communities of experimentalists, optical physicists, geneticists, and data scientists work together to decipher the neural code. The main difficulty encountered today in neuroscience is to handle a huge diversity of techniques that are rarely compared exhaustively, drawing parallels between works carried out at different scales, on different model animals, with different imaging methods, different stimuli and different based on a measure of entropy that shows the degree of coordination between neurons.

Self-organized criticality is a cornerstone of the analysis on neural networks dynamics [START_REF] Hesse | Self-organized criticality as a fundamental property of neural systems[END_REF] [START_REF] Bak | Self-organized criticality: An explanation of the 1/ f noise[END_REF]. The sandpile model of self-organized criticality [START_REF] Massobrio | Criticality as a signature of healthy neural systems[END_REF] illustrates the emergence of the neural code in neural networks. Criticality, largely used in Physics, provides a model that links the statistical properties of the avalanches (size, duration, interval) with scale-free power-law distributions [START_REF] Beggs | Neuronal Avalanches Are Diverse and Precise Activity Patterns That Are Stable for Many Hours in Cortical Slice Cultures[END_REF]. It allows a classification of the activation of a neural network in different regimes close or far from chaos: [sub-]critical. It modulates the phase-transition between different stimulus representations.

Topology theory

During its development, the network structure evolves jointly with dynamics. The network needs to adapt depending on the way the information flows, is integrated or distributed. Neural connectivity shares common topological patterns with other large-scale complex networks from politics, economics, or social networks. The fine analysis of the emergent structures unveils statistical distribution, patterns, connectivity that share common key organizational principles. Modularity, hierarchy, centrality, and hub existence result from the evolution of networks and are the first observables required to uncover its intrinsic mechanisms of neural code and computation.

The existing network topologies are generally classified in several categories that result in several generative model of graph structure:

• Random Network (Erdòs-Renyi) [START_REF] Erdős | On random graphs. I[END_REF] : Every pair of nodes have the same probability of being connected. Node degree distribution of large networks follow a Gaussian distribution.

• Scale-free Network (Barabàsi and Albert) [START_REF] Albert | Statistical mechanics of complex networks[END_REF] The node degree distribution follows a power law distribution ∀I ∈ [1, . . . , N ], P (d i = k) ∼ k -lambda /λ ∈ [2,3].A very gradual "heavy-tail" power law decay of the degree distribution implies that the network lacks a characteristic scale. Preferential attachment phenomenon is observed. These networks are unlikely to emerge for network with physical link constraints.

• Regular lattice Network: Extremely regular network defined by a high-level of clustering and a long path length.

• Small-World Network (Watts-Strogatz) [START_REF] Watts | Collective dynamics of 'small-world' networks[END_REF]: Structure largely observed in neural networks. It combines high levels of local clustering among nodes (cliques/families) with a low average short path. On average, only few direct connections remain but information can spread easily due to the interconnected modules.

Graph theory measures from [76]

These topological structures can be supported by some topological features.

• Modularity: The modules correspond to densely interconnected sub-graphs relatively to the connection density between nodes of different modules. Thus, modularity can be estimated by comparing the number of connections insides modules compared to the number of connections expected in a random graph. • Hierarchy: Hierarchy is the successive nesting of modules with different node densities, reflecting the organizational structure of a multi-layer or multi-level network.

• Centrality, and robustness: The centrality of a node measures the number of shortest paths which go through nodes. Robustness refers to the structural integrity of the network after the deletion of nodes.

All that topological features can be quantified in a network represented as a graph using topological feature evaluators:

• Node degree, degree distribution and assortativity: The node degree is the number of connections that link the node with the rest of the graph. The histogram of all node degrees in a graph forms the degree distribution. For a long time, scientists expected to have a Gaussian distribution symmetrically centered to model node degree in complex graph. However, for a lot of real-world networks, like in vivo and in vitro neural networks, this distribution needs to be discarded. For example, the number of sexual interactions between human beings represented in a graph follows a power law distribution like neural networks does. It means that the number of people with an important number of sexual partners is too high to be matched with a gaussian distribution. This is a phenomenon called "rich-gets-richer" or "preferential attachment" [START_REF] Albert | Statistical mechanics of complex networks[END_REF]. It means that people are more attracted by people that already have lots of established links. Similar interactions can be found at neuron scale. This statistic mediates the propagation of a sexual disease in the population like the information flow in a neural network. Preferential attachment is related to development of the power-law scaling in the scale-free and small-world topologies. Assortativity is the correlation between the degrees of connected nodes. Positive assortativity indicates that high-degree nodes are interconnected.

• Average Path length: The path length is the number of edges to go from one node to another. Random and complex networks have short path length unveiling a high global efficiency of information transfer and robustness compared to regular lattice networks. Efficiency is inversely related to path length but easier to estimate.

• Clustering Coefficient [START_REF] Watts | Collective dynamics of 'small-world' networks[END_REF] and motifs: a cluster of a nodes corresponds to all nearest neighbor nodes connected to each other. The clustering coefficient is the number of connections between the nearest neighbor nodes and the maximum number of possible connections. Random networks have low clustering coefficients whereas complex networks have high. It encodes a high local efficiency of information transformation transfer and robustness. The detection and classification of small motifs of interconnected nodes can, equally, be estimated to identify the possible local interaction supported by the network e.g., the cliques or the k-clubs.

• Connection density or cost: The connection density is the number of effective links as a proportion of the total number of possible connections related to the energy of the network.

• Betweenness, Closeness: The betweenness-centrality is the number of shortest paths going through that node divide by the total number of shortest paths between pair of nodes. This metric is used to detect nodes making the bridge in the network. Closeness-centrality measure the speed to reach all other nodes of the network from one node. It corresponds to the average of the shortest paths between the node and every other one.

At the neural network scale, a small-world topology is observed. It is costly for the network to create long and robust connection between distant sites. Thus, the spatial layout of the neural networks tries to minimize the axonal volume by creating densely connected subgraphs corresponding to modules interconnected with each other. It reveals a high local efficiency, robust to failures and associated with a short average path length, allowing each module to communicate efficiently. Lots of neuronal microcircuits have also been identified, leading to the formulation of probabilistic connection rules. The high modularity of neural networks is associated with a hierarchical organization of a functional network that generates diverse and persistent dynamic patterns. Modularity is supported by subgraph of interconnected nodes existence in networks [START_REF] Betzel | Multi-scale brain networks[END_REF] [START_REF] Boccaletti | Complex networks: Structure and dynamics[END_REF] . For in vitro neural networks, modules have small number of neurons. Their size is positively correlated with neuron spatial dispersion. Nodes belonging to the same module are more synchronized. Hub nodes with higher degree connections exist in the structure, having a huge influence on the network. They are at the core of the modules being extremely synchronized with their surrounded environment.

Structural organization of in vivo neural network is closely related to functional connectivity that share common topological features such as modules and hubs. Structural topology abnormalities have been related to neural diseases at the brain level e.g., schizophrenia linked to a reduction of hierarchy and longer average path compatible with inefficient axonal wiring. At neural network scale, the deletion of high-degree nodes (hubs) could be promising to study how dynamic of the system evolves during a failure and how does the topology restructure.

Biological Evidence of dynamics and topology emergence in in vitro and in vivo experiments

Lots of attempts to study emerging neuron network properties during development are now under the scope. Especially, over the last decades, a lot of work has been done to design in vitro interconnected neuronal assemblies [264] [265] [266] [START_REF] Kanagasabapathi | Functional connectivity and dynamics of cortical-thalamic networks co-cultured in a dual compartment device[END_REF] [268] [START_REF] Pan | An in vitro method to manipulate the direction and functional strength between neural populations[END_REF] and visualize in vitro the dynamical and topological property emergence. Levels of activity, neuronal avalanches, network organization, connectivity topology [START_REF] Pajevic | Efficient Network Reconstruction from Dynamical Cascades Identifies Small-World Topology of Neuronal Avalanches[END_REF] [271] [START_REF] Massobrio | Criticality as a signature of healthy neural systems[END_REF] are determinant features inherent to different neural states that will shape different dynamics of the system with an inter-contextual and inter-individual variability.

At early stages, dissociated neuronal assemblies are expected, characterized by a lot of small bursts leading to a persisting activity during the entire network life time [START_REF] Marom | Development, learning and memory in large random networks of cortical neurons: lessons beyond anatomy[END_REF]. At that step, information processing is optimized to work using a restrict number of notfully connected computation units [START_REF] Shew | Information Capacity and Transmission Are Maximized in Balanced Cortical Networks with Neuronal Avalanches[END_REF]. The Self-Organized Criticality (SOC) [START_REF] Bak | Self-organized criticality: An explanation of the 1/ f noise[END_REF] [103] [START_REF] Heiney | Neuronal avalanche dynamics and functional connectivity elucidate information propagation in vitro[END_REF] [274] connects the microscopic and macroscopic levels of interpretation of neuronal system [START_REF] Hesse | Self-organized criticality as a fundamental property of neural systems[END_REF] [START_REF] Poli | Functional connectivity in in vitro neuronal assemblies[END_REF]. The immature network topology is almost completely random with only few interactions that might be due to the absence of a sufficient number of connections to form a more complex structure and the non-synchronous events could pollute the Spike-Timing Dependent Plasticity (STDP) preventing synapses from reinforcing.

During maturating phases, periods of activity increase corresponding to more stable neuronal avalanches and disclose dissociated neuronal network existence. Networks are structured with an edge density naturally increasing as dynamics emerge. A convergence in topology and dynamic appears and communication scale changes going from individual neurons to neuron groups spread in the networks at longer distance. An increase of the degree of small-worldness is observed unveiling deeper and deeper levels of hierarchical arrangement [275] [276]. Downes et al, 2012 argued that the small-world topology emerges from a random network via "synchronized bursts" whereas Antonello et al [START_REF] Antonello | Self-organization of in vitro neuronal assemblies drives to complex network topology[END_REF] claims that it does from the promotion of neurite outgrowth. The successive establishment of new preferential links during the development could be the factor that explains the modification of the observed topology. These links are, first, established locally and the synchronization patterns between nearest neighbor neurons mediate their activity and growth. It means that spontaneous activity, self-maintained by the exchanges between modules would also appear as a stabilizer agent. During maturation phase, connections then develop to communicate at a higher distance, but the average length path does not change significantly. For Antonello, it is the clustering coefficient that seem to contribute to the SW topology emergence that reflects the ability of the network to maximize communication efficiency by minimizing the wiring cost. Watts and Strogatz proposed another theory to account for its emergence. It is based on a random rewiring scenario that could create shortcuts and lead to the observed increase of local clustering.

At longer time scale, neural networks demonstrate multi-level synchronization patterns that seem mediated by balanced excitatory-inhibitory subsystems. Strong structural connection makes functional connection more likely to occur. Finally, functional topology converges from random to small-world with the emergence of associated topological features such as neuron hubs and rich-club.

Neuroscientists shaped different neuron population sizes, type of population, number, and size of (sub-)modules. Hardware material have been made to drive this kind of varying structures and compare the emergent functional connectivity structures. Furthermore, some electrical stimulations can be applied on neural cells during development to modulate emerging neuronal dynamics. For instance, variations on the number of stimulated sites, frequency stimulation, or amplitude of the pulses to change synaptic efficacy were proved to be able to modulate neural network dynamical states. As external electrical stimulation increases, a stabilizing of the integration is observed rather than segregation processes compared to spontaneous activity. External electrical stimulation, regardless of the stimulation sites, affects the number of functional links, as well as the average magnitude of changes. A trial-to-trial variability is observed for a same external input stimulation. The exact influence of stimulation on topology evolution is not fully solved yet, but it seems that initially random networks converge toward a modular structure in the same way, but synaptic strength can be reinforced or altered [START_REF] Wang | Community discovery using nonnegative matrix factorization[END_REF] [START_REF] Wang | Alterations in Brain Network Topology and Structural-Functional Connectome Coupling Relate to Cognitive Impairment[END_REF]. Some evidence have also highlighted similar dynamics and topology emergence in in vivo neural networks.

Hydra Vulgaris is known for its ability to reaggregate from individual cells by recreating an entire animal. This ability of complete reaggregation after dissociation makes this animal a golden-standard study case to analyze the emergent computational properties of neural networks in vivo. From initial random topology structure with a random dynamic, a modular and hierarchical topology emerges, characterized by synchronization patterns. In the reassembly of Hydra Vulgaris' nervous system, transitions between regimes of self-organized criticality would be at stake supported by hierarchical modularity (see A.4).

After 72h, the nervous system of dissociated cells is fully reestablished across the entire animal. Synchronization patterns are appearing with an increase of the functional connectivity strength. Network activity becomes distributed during reaggregation. From 8h to 24h, a hierarchy is emerging directly demonstrated by the relationship between clustering coefficients and node degree corresponding to a scale-free structure. A transition of structure is observed from this hierarchical structure to a distributed one that facilitates this extreme regeneration allowing immature modules to form and resume their function less discriminatory than more highly patterned systems via hub neurons.

The global distance of module node decreases over time. At 24h, modularity becomes more stable with smaller, densely connected ensembles. A synchronization needs to emerge, and it happens through independent sub-modules. It seems to be a strategy to optimize this process by recreating the activity locally first.

Neuronal avalanches whose statistics are in accordance with a phase transition in critical regime is observed. From 24 to 48h, synchronization emerges between modules at the scale of the entire animal. Coactivity of modules lead to module fusion and a loss of modularity at 72h. Some submodules remain isolated corresponding to asynchronous activation. Submodules seem to act as subsystems that can trigger avalanche effects in other modules. Surprisingly, a noticeable decrease of hierarchy is observed in the experience that could account for the existence of hub nodes mediating interconnections between proto-modules and immature ensembles formed on their way to the intermediate state of maximum modularity. These nodes could contribute to Completion Pattern. [START_REF] Carrillo-Reid | Controlling Visually Guided Behavior by Holographic Recalling of Cortical Ensembles[END_REF] Unpublished work about the stimulation of the completion pattern neurons to recall neuronal ensembles via optogenetics at Hydra Vulgaris are, now, under the scope.

To summarize, the connection strength is firstly local, and its distance projection increases during neural network development. This connectivity strength increases during the animal's aggregation leading to a much smaller distance between modules. Second, hierarchical modularity increases revealed by the increase of clustering coefficients. Thirdly, power-law statistics showing that the in vivo network converges toward a small-world critical regime of activity allowing phase transition. Carrillo-Reid [START_REF] Carrillo-Reid | Controlling Visually Guided Behavior by Holographic Recalling of Cortical Ensembles[END_REF] went further in the stimulation of neuronal ensembles in mice visual cortex. The use of 2-photon holographic Optogenetics in mice primary visual cortex can be used to excite Completion Pattern Neurons (CPN) and trigger behaviorally relevant neuronal ensembles. CPN correspond to high interest nodes that are the cores of the synchronization patterns and have the intrinsic ability to trigger the activity of an entire module or neural circuit in correlation with a specific behavior. CPN are basically hub neurons proof of the existence of a modular structure. The article demonstrates the ability to alter behavioral performance. The activation of neurons unrelated to perceptual behavior can degrade behavior while CPN stimulation can trigger the behavior without a real external stimulus. The interpretation is that a neural ensemble encodes an internal representation of a stimulus. The perception of a stimulus could be internally driven using cortical states. Thus, ensembles could be dynamical attractors that implement internal, perceptual or memory states. It seems that the world representation can be altered by interacting with these states by activating specific neurons. Some future works will necessarily emerge to inactivate, via Optogenetics, CPN to, firstly, estimate if a specific behavior could, robustly, be inactivated (s.t Contraction Burst at Hydra Vulgaris )and then if pathological neural diseases could be reduced by killing the potential involved CPN in its neuro-circuitry.

As a conclusion, regular and irregular synchronization patterns seem to shape the network topology regarding neuron ensemble dynamics. Spontaneous activity pattern drives the topology transformation via the self-maintained modular activity. The emergence of Small-world topology via Hebbian learning rule and activity driven plasticity seem to be Step of the Louvain algorithm from [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF] ensembles during reaggregation and development to demonstrate the birth of the emergent computational properties of the neural networks at Hydra Vulgaris.

Stochastic Block Model [128]

Stochastic Block Model is a very promising mathematical tool that combined Statistical Inference and Graph theory to extract statistical features about the connectivity of neural networks. Lots of versions and adaptations exist [START_REF] Jin | A Survey of Community Detection Approaches: From Statistical Modeling to Deep Learning[END_REF].

Avitan [START_REF] Avitan | Spontaneous Activity in the Zebrafish Tectum Reorganizes over Development and Is Influenced by Visual Experience[END_REF] analyzed the evolution of the spontaneous activity response of zebrafish larvae's visual cortex during development using a version of SBM called Degree-Corrected Stochastic Block Model (DG-DBM) [START_REF] Newman | Estimating the Number of Communities in a Network[END_REF].

A statistical model is set on the connectivity graph. Each node will have a probability to generate an edge with every other neuron regarding a specific distribution. The aim will be to inverse the statistical problem by estimating the statistical features that is the best to describe the observed network. Let A be the adjacency matrix. a ij = 1 if node i and j are connected by an edge, 0 otherwise. The membership vector is g where ∀i ∈ [1, . . . , N ], g i is the index of the community neuron i belongs to. Let k be the number of ensembles supposed to be fixed and known. Let γ r be the probability of assignment to community r. (∀r ∈ [1, . . . , k]) It means that each ensemble is characterized by a probability of recruitment fixed and identical for all neurons. We have the normalization condition assuming the disjoint partition of membership: r=1 kγ r = 1 which means that all nodes are assigned to a single group following a categorical distribution. Given a neuron i and a neuron j membership that belong respectively to community r and s, every pair of nodes will have a probability to form an edge. This linking follow a Poisson distribution with mean w rs ∀(r, s) ∈ [1, . . . , k] 2 . The model is summarized by its likelihood function to generate a particular network A, with group assignment g, given the parameters γ, w, and k. Where n r = i σ g i ,r is the number of nodes in group r. σ g i ,r = 1if g i = r, 0otherwise. m rs number of edges between community r and s. m rs = ij a ij δ g i ,r δ g j ,s /r ̸ = s.

The likelihood is divided in three terms. First term summarizes the categorical distribution of the membership corresponding to the membership of the n nodes in the k clusters. Second term represents the creation of the links of the same community. The third term corresponds to the creation of links between neurons belonging to distinct communities. Priors are set on the statistical parameters some of them directly estimated from the observed adjacency matrix known as "Empirical Bayes" techniques while others from uniform distribution. Then, a MCMC sampling routine is performed to estimate the most relevant statistical parameters. This method is an ordinary stochastic block model that can give a poor fit for most real-world network data since it fails to match the broad degree distributions commonly observed in data. Here, it is indeed, assumed, that the number of links is similar for nodes belonging to same ensembles. In real data, huge heterogeneity of node degrees is observed corresponding to scale free topologies highlighted in in vivo neural networks. It means that the statistical model is too simple compared to real-world data complexity. To account for that phenomenon, the degree correction part has been implemented as an extension of the model.

θ is a vector of parameters that describes for each node its ability to create link allowing freer degree distribution. The expected number of edges between a pair of node i and j becomes θ i θ j w rs . (node i and j belong respectively to r and s) A normalization is added to avoid the statistical non-identifiability of model parameters: ∀r ∈ [1, . . . , k], 1 nr i θ i σ g i ,r = 1. Leading to an average number of edges between two ensembles r and s equals to ij θ i θ j w rs σ g i ,r σ g j ,s = n r n s w rs . w rs is still the average probability of an edge between nodes in groups r and s but nodes with released degree can exist. The new likelihood function can be calculated: Where d i = j a ij is the observed degree of node i. Prior distributions are set. A MCMC algorithm is used to iteratively sample the parameters. In some cases, the algorithm can get stuck in "metastable states" where equilibrium is reached extremely fast. Repeated runs with random initialization are then performed keeping the result that provides the highest average likelihood value.

The limitations of the techniques are to be applied using a statistical model of the graph connectivity and not directly a model of the interaction of neurons (obtaining the observed adjacent matrix A still raises the question of the first steps about functional connectivity) and the method does not handle overlapping in the membership part. The mathematical framework remains still very promising highlighting how statistical inference can be used to estimate relationships between nodes given an observable.

Spectral Technique Taxonomy

PCA / ICA [200] [281] [202]

The aim of the component analysis methods is to generate an abstract representation of the data in a new mathematical space where data complexity is diminished allowing easier pattern and relationship identification. Independent component analysis (ICA) and Principal Component Analysis (PCA) are respectively two machine learning techniques able to perform independent source analyses and data dimensionality reduction.

PCA search for the orthogonal directions in the data that maximize its variance as the eigenvectors of the covariance matrix of the spike connectivity matrix. Spike matrix has been normalized using a z-score providing its covariance matrix as the correlation matrix of the spiking activity. This matrix is given by:

C = ZZ T N columns (B.5)
Where Z is the binned spike count matrix got after the statistical z-score. According to the spectral theorem, C is real and symmetric, so diagonalizable.

C = λ i x i x T i (B.6)
Where x i is the i -th eigenvectors of C (the i -th PC of Z) and λ i the corresponding eigenvalue. x i x T i is the projection matrix onto the direction x i and λ i its variance. Under the null hypothesis of a random matrix, some advanced mathematical works has been done to highlight the properties of the spectral distributions [START_REF] Erdős | On random graphs. I[END_REF] [282] [START_REF] Tracy | Level-spacing distributions and the Airy kernel[END_REF] [284] that allow to discard non-correlated interactions. Especially, the eigenvalues of the correlation matrix of a normal random matrix M with statically independent rows follow a probability function described by the Marcenko-Pastur distribution:

P (λ) = q 2πσ 2 (λ max -λ)(λ-λ min ) λ (B.7)
With q = N columns Nrows ≥ 1, σ 2 the variance of the elements of M (equals to 1 for z-score normalization), N columns the number of columns of M, N rows the number of rows of M. λ max/min the maximum and minimum bounds of eigenvalues such that λ max/min = σ 2 (1 + -1 q ) 2 . D = [λ min , λmax] (P (λ / ∈ D) = 0). λ max can be used as a threshold to estimate the number of ensembles. MP distributions have been proved for random matrices whose entries are derived from gaussian distribution and provides a good bound for eigenvalues of matrices composed by independent rows corresponding to uncorrelated neurons. A finite sample bias correction can be applied by using Tracy-Widom distribution instead used for instance by Bickel and Sarkar [START_REF] Bickel | Hypothesis testing for automated community detection in networks[END_REF] do derive a statistical hierarchical ensemble detection algorithm.

Assembly patterns and activity estimation. In linear model, the neuron ensemble activity is a linear summation of the activation of the cells. The co-activation pattern is a vector corresponding to the weights of each neuron in the specific contribution of that pattern. For all time bin b, the activity of cell assembly is then:

R b = i=1 N w i z ib = w T Z b (B.8)
Where N is the number of neurons, Z ib is the Z-score activity of neuron i at time bin b, w i is the weight of neuron i in the assembly, w is the column vector of all neuron weights and Z b is the b-th column of matrix Z. By projecting the neural activity at each time bins regarding the weighted contribution of the neurons in orthogonal decorrelated assemblies, it allows to analyze if different activity patterns emerge. Thus, each cell assembly activity is estimated by this projection of columns of the spiking matrix onto the axis spanned by the corresponding assembly pattern extracted as the PC of the PCA. The projector is P = ww T where w is the unitary vector of assembly patterns that spans the axis. The projection acts like a similarity measure between the activity of the whole population and the assembly patterns. This method has same limitations since in some cases when assemblies overlap for instance estimated ensemble activities can be mixed up [START_REF] Lopes-Dos Santos | Neuronal Assembly Detection and Cell Membership Specification by Principal Component Analysis[END_REF]. In the specific case a neuron can show a large weight contribution for different orthogonal PCs. Thus, if two assemblies share a similar amount of variance in their corresponding axes, it becomes possible that first PC represent the average of two assemblies and not an individual assembly anymore. Then, PCA fails, and further solutions need to be implemented.

Assembly vector estimation is a solution. PC weights carry information about the membership since larger weights encode larger neuron contribution to the axis. The question is to define a membership threshold. AV is a technique to extract both neuron ensemble membership and ensemble activity patterns limiting mixing the ensemble activations. In the AV framework, the number of neurons composing at least one assembly is estimated counting the number of eigenvalues that lie outside Marcenko-Pastur distribution boundaries. The columns of the correlation matrix are projected onto the subspace spanned by the PCs associated with significantly large eigenvalues. Where P C i is the i-th significant PC and P c is a matrix containing all significant PCs (columns). The subspace of projection is the Assembly Space (AS) and the columns of the correlation matrix N i projected in AS are called neuron vectors. The number of eigenvalues that lie outside the theoretical distribution defined the number of significant neurons to handle. The AS gathers neurons with neuron vector corresponding to neuron population [START_REF] Carrillo-Reid | Cell Assembly Signatures Defined by Short-Term Synaptic Plasticity in Cortical Networks[END_REF] [START_REF] Carrillo-Reid | Identification and Targeting of Cortical Ensembles[END_REF] with similar correlation patterns. Neuron with orthogonal correlation patterns will be separated in that space. The neuron pattern similarity matrix called Interaction Matrix is calculated by calculating the inner products of all neuron vectors and those of significant neurons provide the co-activation patterns similarity.

M i,j = N T i N j (B.11)
The distribution of similarity co-activation patterns is expected to be ideally bimodal allowing a binarization separating to noise to significant information. A homemade clustering algorithm is run on the binary matrix to find the cluster of neurons.

Finally, the assembly vectors (AV) are defined as the barycenters of the neuron population N i whose neurons are exclusive to an assembly.

AV a = 1 n a i N i (B.12)
The method is better than a simple PCA but limited by the separation of small and large inner products in the AS. The bi-modal representation cannot always be simply divided by a hard threshold. The method requires the existence of exclusive neurons that could not exist.

The combine use of PCA with ICA is another implementation made to improve the detection of neural assemblies in the absence of exclusive neurons. This work is inspired by [START_REF] Laubach | Independent component analyses for quantifying neuronal ensemble interactions[END_REF]. ICA is a method used to extract non-gaussian components from multivariate signals. Rigorous mathematical formalism of ICA can be easily found but the basic idea uses the central limit theorem. Summing lots of independent random variables provides a gaussian distribution. The more numerous is the number of random variables, the higher the "Gaussianity" of the signal is. ICA projects the data in a space where the data Gaussianity decrease to find independent components that intervene as a linear summation in the resulting signal. A rotation invariance needs to be avoided corresponding to the extraction of gaussian signals that are impossible to recover using this technique. The first part of the algorithm remains the same corresponding to the extraction of the PC and then they are rotated to match ideal assembly patterns that shows significant improvement in the methods. ICA cannot handle gaussian distribution or non-linear spike correlation. Neuronal assemblies are bin size dependent.

SVD

Carrillo-Reid did a huge work to estimate neural ensembles, identify neuron with specific roles (= CPN) and recall the ensembles during behavioral tasks. In his vision, at each time step a neuron population is activated. The aim is to find similar neuron population that are repeatedly activated over time. The method is, thus, based on a spatial and temporal criterion of similarity. His method can be divided into successive steps: 1) Estimation of the high-level synchronous neuron population activities. 2) Normalization using TF-IDF transform to remove non-significative neuron population activation. 3) Mapping of neuron population using Cosine Similarity calculation. 4) Population similarity Binarization. 5) Temporal Profile similarity calculation using Hamming distance calculation. 6) SVD for neuron clustering and ensemble activation estimation. Singular Value Decomposition (SVD) factorization is a generalization of the PCA when we project the data matrix in the highest variance directions instead of its covariance matrix. SVD is more general and does not reduce information dimensionality during the factorization step. For SVD, the number of neuronal ensembles is given by the magnitude of the singular values, so an estimation of the number is not required using a formal mathematical distribution under a null model. A direct mathematical link can be highlighted between SVD and PCA.

Compared to previous work, Carrillo-Reid highlights temporal properties of network activities using a representation of the overall activity as an array of multidimensional population vectors based on time points analysis. The identification of recurrent groups of cells firing together is at the core of the method. Some user-defined thresholds are used in the similarity metrics such as binarizing the matrix to apply the hamming distance on it. The method has not been properly benchmarked against previous spectral techniques and its The matrix is normalized using a TF-IDF transform from the NLP (Natural Language Processing) community. It limits the influence of too activate cells that are not statistically significant. c) A similarity matrix is calculated using a pairwise cosine similarity. d) The matrix is binarized and decomposed using a SVD. It unveils orthogonal components to highlight dissimilar population that repeat over times. e) The decomposition is perfomed in an orthogonal basis. f) An estimation of an ensemble activity matrix is performed and matched with the mice visual stimulations. g) It highlights frame by frame the activated population in response to stimulations. accuracy is unknown in specific experimental conditions such as a high level of ensemble overlapping.

Model based techniques

A model-based technique is a technique that explicitly defined what a biological neuronal ensemble is. It handles the specificity of the problem providing interpretable estimates. These techniques are the most likely to be reused by experimentalist since being meaningful. The main goal of such techniques is to explicitly formulate the model and to try to solve an inverse problem by estimating the most likely model variables given observables.

To our opinion, the most promising technique is the technique from [START_REF] Diana | Bayesian inference of neuronal assemblies[END_REF], that formulates a statistical model with neuronal ensemble definition based on synchrony whose estimates are infer using Bayesian inference to solve the inverse problem. Applied on zebrafish larvae, this technique offers according to our opinion, new perspectives to estimate the emergent properties of synchronization emerging from neural networks and unveiling neuronal ensembles.

This technique is developed around 3 major working visions:

• Neurons are divided into several assemblies. Neuron membership is unique and only defined by a probability n k of recruitment of each neuronal ensemble. Statistically the membership is modeled by a categorical distribution.

• Neuronal ensemble can be active or inactive independently at each time step. It means that a neuronal ensemble firing is defined by the coactivation of a subset of its neurons. This probability activation of each ensemble is model by a Bernoulli law with probability p k assumed to be fixed over a recording.

• Given the state of an ensemble at time bin k, active or inactive, every neurons will have a probability to respond (a-)synchronously. It means that given an ensemble firing, its neurons could fire with it or not. Outside neural ensemble activity every neurons have the same basal probability to be activated corresponding to a spontaneous activity. Where n µ is the probability of recruitment of the community µ, p µ is the probability of activation of community µ, and λ µ (w k,µ ) is the probability of activating the neurons belonging to community µ given the activation state of community µ at time k. (λ t i (z) = P (s i,k = 1|w k,t i = z)) t and W are the hidden variables of the model corresponding respectively to the neuron membership vector and the ensemble matrix probability.

Once the model is explicitly formulated, its inversion and the estimation of the most likely statistical parameters and hidden variables needs to be carried out. It requires an inversion routine. This inversion is based on bayes formula P (Z, θ|data) ∝ P (data, Z|θ) × P (θ) and 2 versions of the algorithm based on two model, and two optimization procedures: the Collapsed Gibbs sampler routine and the Metropolis-Hasting algorithm. Both algorithms are based on the establishment of prior distributions, chosen to be conjugate priors providing sampleable posterior distributions allowing the optimizing procedures. Some posterior distributions are marginalized out to collapse the algorithm and proving a more efficient convergence. The posterior distributions calculated are:

, 0 < u i 1 < u i 2 < • • • < u i J i . N i (t)
is the number of spikes that fall in the time segment [0, . . . , t], ∀t < T . The Conditional Intensity Function (CIF) describes how the spiking past history of neuron i influences its spiking activity at time t. It corresponds to its instantaneous probability to fire at time t given its covariates H i (t).

λ i (t|H i (t)) = lim δ→0 P [N i (t + δ) -N i (t) = 1|H i (t)] δ (C.6)
It corresponds to an instantaneous conditional firing rate. H i (t) denotes the spiking history of all the neurons in the ensemble up to time t for neuron i. H i (t) defined in the interval [t-M i W, t] such that M i (t) non-overlapping time bin windows of width W . For all neuron q = 1, . . . , Q and all timebins m = 1, . . . , M i R q,m (t) denotes the spike count of neuron q in time window of length W convering the time interval [t -mW, t -(m -1)W ]. In the continuous or discrete case, Granger Causality embeds dynamics model through respectively an Auto-Regressive Model (VAR model) and a Global Linear Model (GLM). In discrete time series events, we have :

λ i (t|γ i , H i (t)) = exp   γ i,0 + Q q=1 M i m=1 γ i,q,m R q,m (t)   (C.7)
It means that the probability of neuron i to spike is a linear integration of basal spontaneous activity of the neuron γ I,0 and a weighted contribution of ensemble spiking R q,m(t) . This weighted contribution is defined by a parameter vector γ i γ i = [γ i,0 , γ i,1,1 , . . . , γ i,q,m , . . . , γ i,Q,M i ] (C.8)

Which represents the dependency of neuron i the spiking history of all neurons in the ensemble. This parametric CIF is fit on the data using its point process likelihood function. The likelihood of the neuron spike train i is obtain modeling the process as a Bernoulli distribution of spike in each time bin :

L i (γ i ) = K k=1 [λ i (t k |γ i , H i [k])∆] ∆N i [k] [1 -λ i (t k |γ i , H i [k])∆] 1-∆N i [k] + o(∆ J i ) (C.9)
o(∆ J i ) denotes the probability that neuron i includes more that one spike in the sub interval [t k-1 , t k ]. Measuring the effect of a neuron j to a neuron i means evaluating if its spiking history has a significant statistical impact on the prediction of the probability of neuron i to fire at time t. It means that this notion of "more likely to" is encoded by a ratio calculation between two likelihood functions. One which accounts for neuron j influence, one that does not. The likelihood ratio Γ i,j is thus given by :

Γ ij = log L i (γ j i ) L i (γ i ) (C.10)
This ratio is a coefficient that answer to the question : Is it more likely that neuron j influence the instantaneous conditional probability of neuron i to fire at time t or not. L i (γ j i ) is the likelihood function of the spike train of neuron i where the contribution of neuron j is removed :

log λ j i (t k |γ j i , H j i [k]) = γ j i,0 + Q q=1 q̸ =j M i m=1 γ j i,q,m R q,m [k] (C.11)
Finally the Granger causality coefficients between neuron j and i are given by the following formula providing a Q × Q matrix summarizing all causal influences between neurons.

ϕ ij = -sign   M i m=1 γ i,j,m   Γ ij (C.12)
This coefficient summarizes the interaction between neuron i and neuron j. A positive result, a negative result and a close to zero result reveal respectively an excitory, an inhibitory and a no influence effects. Applying this metric on every pair of node provides a connectivity matrix of size Q×Q. Several versions of this algorithm exist but some recent studies highlight advantages and drawbacks. 7) Tranfert entropy : We focus on the original definition made by [START_REF] Schreiber | Measuring Information Transfer[END_REF] [296] [297] [START_REF] Timme | A Tutorial for Information Theory in Neuroscience[END_REF] and applied by [START_REF] Ito | Extending Transfer Entropy Improves Identification of Effective Connectivity in a Spiking Cortical Network Model[END_REF] in a context of spiking point process. TE is a directed, asymmetric, measure of interaction that can be applied on two time series. TE provides positive coefficients to account for how including information of neuron j's spiking activity improves the prediction on neuron i's activity compared to its own history only. It is an information flow measure. i t denotes the binary status of neuron i at time t. k and l are the order of TE. In this framework, the probability are estimated counting the effective number of situation events occurring.

T E

Figure I. 1 :

 1 Figure I.1: Original Prospective Processing Pipeline. 1-Extraction of single neuron fluorescence traces from imagery using tracking algorithm. 2-Spike inference to estimate the dates of action potential emissions. 3-Extraction of a network topology from neuron spiking interactions. 4-Detection of overlapping neuronal ensembles. 5-Estimation of sequential patterns of neuronal ensemble encode behaviors. 6-Prediction of behavior at Hydra vulgaris based on the inferred computational rules.

Figure I. 2 :

 2 Figure I.2: Mice experimental setup of visual stimulation and recording of evoked neural response in primary visual cortex V1. [4]. a) Experimental setup. Static gratings and gray screen are displayed in front of an immobilized mice. The neural response in V1 is recorded. b) The different stimulation angles displayed and the response of the brain area is recorded. The detected responding neural cells are displayed in red.

Figure I. 3 :

 3 Figure I.3: Zebrafish experimental setup of visual stimulation and recording of evoked neural response in larval zebrafish optic tectum [3]. a) The stimulation is the display of a spot at different spatial location the animal being fixed in an agarose substrate.b) The neural response in the optic tectum is recorded using a 2-photon microscope. A wave activation reponse is observed for each stimuli but the strength of the response and the overlapping largely vary according to the stimulation. c) A very clear spike response can be recorded for some specific stimulation.

Figure I. 4 :Figure I. 5 :

 45 Figure I.4: Non-overlapping neural networks at Hydra vulgaris from [2]. a) Neural networks at Hydra vulgaris are expected to be non-overlapping. There is no specific spatial organization of the networks in the animal except that some networks can be exclusively located in ectoderm, endoderm or in the neighborhood of the tentacles. b) Functionally the neural network are non-overlapping. It means that a very low-level of co-activation is expected between the microcircuits and a specific neuron belongs to a single network at the same time. a) b)

Figure II. 1 :

 1 Figure II.1: The analysis of the neural code can be performed at different scales from whole brain (Macroscopic) to neuron population (Mesoscopic) and single cells (Microscopic)(adapted from [20]).

Figure II. 2 :

 2 Figure II.2: Different animal models are used to image neuronal activity at single cell level. The visual cortex of mice (data from Yuste's lab) or whole animals like Hydra vulgaris (data from[2]) or larval zebrafish (data from[3], image from[START_REF] Diana | Bayesian inference of neuronal assemblies[END_REF]) using GCamP calcium imaging.

Figure II. 3 :Figure II. 4 :

 34 Figure II.3: Sketch of neurons by S. Ramón y Cajal adapted from [20]. a) Neuron of the mammalian cortex observed under microscope. b) Sketch of S. Ramón y Cajal of single neuron structure and signal transmission from a presynaptic neuron to a postsynaptic neuron.

Figure II. 5 :

 5 Figure II.5: Stereotypical representation of action potentials from a voltage recording adapted from [57]. (top) Difference of voltage between a fine tungsten wire placed near a cell in the fly's brain with a reference electrode placed in the body. (middle) Filtering using pass-band of the signal to highlight the stereotypical shape of the action potentials. On the right the superimposition of the action potentials is drawn. (bottom) Timing pulses generated electronically by a threshold discriminator circuit.

Figure II. 6 :

 6 Figure II.6: Scheme of chemical synapse vs gap junction adapted from [58]. a) Chemical synapse with neurotransmitter from presynaptic vesicles to postsynaptic receptors. b) Electrical synapse with direct gap junction connections.

Figure II. 7 :

 7 Figure II.7: Long-term potentiation vs Spike-Timing-Dependent-Plasticity adapted from [20]. a) Schematic drawing of a paradigm of LTP induction. Different strengths of presynaptic pulses evoke different postsynaptic response imprinting a long-term synaptic coupling. b) Spike-Timing-Dependant-Plastiticity. Pairing of neurons is carried out by stimulating a presynaptic neuron and a postsynaptic neuron with a selective delay to induce a short-term synaptic plasticity.

Figure III. 2 )

 2 , and calcium binds to calmodulin (CaM) (see Figure III.1) whose molecular conformation changes emitting a green photon by fluorescence at a peak wavelength around 510 nm. The method monitors calcium fluctuations associated with action potentials in individual cells.

Figure III. 1 : 115 ]Figure III. 2 :

 11152 Figure III.1: Schematic representation of GFP-based Ca 2+ probes from [115]

Figure III. 3 : 130 ]Figure III. 4 :

 31304 Figure III.3: Spinning Disk Confocal Microscope schematic schematic [130]

Figure III. 5 :

 5 Figure III.5: 2-Photon Microscope schematic from [132]

Figure III. 7 :

 7 Figure III.7: Mice brain visual cortex recorded using TPM from Allen Institute open-source dataset. Frame are summed over the recording unveiling soma's locations.

Figure III. 8 :

 8 Figure III.8: Segmentation of Mice brain visual cortex. The segmentation is performed on the sum of all frames of a video.

  have been carried out to refine the tracklet stitching cost function by integrating some state-of-the-art deep learning techniques Discriminative Features using self-supervised trained convolution neural network (CNN) (see Figure III.10).

Figure III. 9 :

 9 Figure III.9: Tracklets from detections using Wavelet transform in EM C 2 from [135]

Figure III. 11 :

 11 Figure III.11: Hydra Vulgaris before Contraction Burst

Figure III. 13 :

 13 Figure III.13: Hydra Vulgaris deformation estimation from [135]

  going from neuron individual fluorescence traces to individual neuron spiking activity : Template matching [156] [157] [158] [159] , Linear deconvolution [160] [161] [162] [163], Finite rate of innovation [164] [165] , Independent component analysis [166], Non model-based Signal Processing [167] , Supervised Learning [168] [151] [169] [170] [171], Constrained non-negative matrix factorization [172] [146], Active set methods [173] [174] , Convex and non-convex optimization methods [175] [176] [177] [178] [179], interior point methods [180], and Model-based bayesian inference techniques [154] [181] [182] [183] [184] [185] [186] [187] [188] [189] [190] [155].

Figure III. 14 :

 14 Figure III.14: Deconvolution of a calcium fluorescence trace from mice V1 cortex with OGB-1 as calcium indicator [151]

Figure III. 15 :

 15 Figure III.15: Processing pipeline for neural ensemble detection. a) Community detection techniques based on pairwise functional connectivity and threshold. (Graph theory) b) Community detection techniques based on data reduction and projection and clustering via Machine Learning techniques. (Spectral Methods) c) Community detection techniques based on an ensemble interaction model using spiking patterns.

Figure III. 16 :

 16 Figure III.16: Pose Estimation on animal models. a) On mice [53] [214] b) On Hydra from E. Martin et al not published.

Figure III. 17 :

 17 Figure III.17: Behavior prediction scheme adapted from [214] and E. Martin et al (not published). a) The behavior is recorded through the interdistance variation between control points and classify to fit a behavioral model. b)The behavioral model is defined by 3 states : elongation, contraction, and rest. Any kind of classifier could be used to make the distinguish between the states. For Hydra, a naive hard-threshold is applied on the distance between control points. c) A predictive algorithm is implemented to predict the behavior. For Hydra, a LSTM is used. The provided input is the average pixel intensity of parts of the image. Variations of the input parts are performed to estimate the level of predictability of the pixel intensity of the different parts of the animal. The idea is to unveil if a specific part is involved in a specific behavior prediction.

Figure 2 -

 2 Inferring spikes from fluorescence traces:The first contribution is the benchmarking of the most promising state-of-the-art spike inference techniques on a calcium fluorescent simulator to estimate the most efficient method on data whose properties match Hydra vulgaris' firing statistics and non-linearities. The proposed versatile simulator embeds non-linear dynamics, non-homogeneous Poisson firing statistics with several firing regimes.

Figure IV. 1 :

 1 Figure IV.1: Benchmarking of ensemble detection methods applied on in vivo neural network of zebrafish larvae from [3]

Figure 1 :

 1 Figure 1: Single Fluorescence trace of a neuron from Hydra Vulgaris. a) Calcium uptake of Hydra neurons displayed using Icy software [8]. b) Temporal response of one Rythmic Potential (RP) neuron with the spike pattern extracted manually.

Figure 2 :

 2 PhotobleachingPoisson-Gaussian noise modeling

Figure 3 :

 3 Figure 3: Accuracy of deconvolution methods over simulation parameters. (a) Variation of the detector gain α (i.e. level of noise). (b) Variation of the firing rate in burst. (c) Variation of the baseline amplitude.

Figure 4 :Figure 5 :

 45 Figure 4: Spike Inference techniques performed on real Hydra Vulgaris' Rythmic potential (RP) and Contraction Burst (CB) neuron fluroscence traces.

Figure 1 :

 1 Figure 1: Statistical analysis of time point-processes A-Stochastic time-point processes (e.g. rasterplot of individual neurons' spiking). B-Multi-time-shift analysis of the coupling between two time point-processes. Depending on the accumulation of points from one point-process around the other for different time-shifts, a coupling probability is assigned to each pair of time points.

Figure 2 :Figure 3 :Figure 4 :

 234 Figure2: Robustness of correlation metrics to different coupling time shifts. Our statistical method (red) is compared with standard correlation metrics (Pearson correlation coefficient (purple), F1 score (blue) and C3I (green)) for increasing simulated coupling level and different time shifts (A-t shif t ∼ N (µs = 0, σs = 0), Bt shif t ∼ N (0, 0.3), C-t shif t ∼ N (1, 0), and D-t shif t ∼ N (1, 0.3).)

𝝀

  -AUTOMATIC THRESHOLDING FOR SPIKE ESTIMATION

2. 2 . 1

 21 Bayes law P (θ|X) = P (X|θ)P (θ) P (X) ∝ P (X|θ)P (θ)

Figure 1 : 90 CHAPTER

 190 Figure 1: Bayesian generative model for neuron clustering

  with T z c = α ( * ) (λ c (z))+ T t=1 ( N i=1 s it δ Zi=c )δ ZiWt=z and T z c = β ( * ) (λ c (z))+ s it )δ Zi=c )δ ZiWt=z c ∈ C(A) represents the combination c of groups among all possible combinations of groups. c ∈ C(Z) represents all possible combinations of activation of groups given the membership.

7 5Figure 3 :

 73 Figure 3: Different steps of synthetic data generation. a) A random clustering is drawn for each neuron in each ensemble using a Bernoulli distribution of parameters α k . b) A random ensemble activation matrix is drawn using Bernouilli distribution at each time step with a fixed activation probability parameter p k . c) Given a conditional probability of neuron activation given all possible combination of memberships, the synthetic neuron spiking matrix is drawn. A Bernoulli distribution is equally used.

Figure 5 :

 5 Figure 5: Estimation of the convergence of statistical estimates on a single synthetic experiment. a) Synthetic neural spiking activity of 400 cells during 1000 timesteps. b) Ensemble activity matrix estimate Ŵ . c) Convergence of the log-likelihood. d) Convergence of the conditional spiking activity of neurons given their combinatory membership. e) Convergence of the reconstructed neuron ensembles spiking matrix using F1-score error rate. f ) Convergence of the neuron membership using ONMI metrix. g) Convergence of the probabilities of recruitment α k of the neuronal ensembles. h) Convergence of the probabilities of neuron ensemble activation p k . i) Neuron membership assignation over iterations.

  The simulation provided under the scope of this article are performed with A = 4, N = 400, T = 1000, α = [0.15, 0.15, 0.15, 0.15], p = [0.1, 0.1, 0.1, 0.1], λ = [λ 0 , λ 1 , λ ≥2 ] = [0.05, 0.8, 1]. For all the simulation parameters variration, the ensemble coactivation rate is fixed at 10%. Here, we can observe the clustering error based on the ONMI metric, ensemble activity error based on the F1-score metric, and relative errors e = |xtarget-x| xtarget

11

 11

Figure 6 :

 6 Figure 6: Automatic ascending scanning of the number of communities based on log-likelihood and AIC criteria. (top) Each time a convergence regime is reached of the loglikelihood function the number of communities A is dynamically increased. (middle)The log-likelihood and AIC reach respectively a maximum and minimum value for the correct number of communities that has been used in the generative process of synthetic data. Here A * = 3 communities. (bottom) The error rates of ensemble activity and neuronal membership reconstructions both provide a perfect estimate for the correct number of communities.

Figure 7 :

 7 Figure 7: Preprocessing step of non-significant spike filtering. a) Raw neural activity for stimuli 8, 4, 5, 9 and 7 between in frames [390, 1571] corresponding to the sequential repetition of 6 trials. b) Filtered and normalized neural activity.

Figure 8 :

 8 Figure 8: Experimental setup of zebrafish larval optic tectum stimulation. Neuronal ensemble activation and statistical coupling with simuli is inferred. a) From [16], we extract the neural activity of the 5 best expressed stimuli (8-4-9-5-7) in the larval zerbrafish optic tectum corresponding to different angle spot locations. b) Superimposition of the pre-processed neural activity with stimuli showing single wave neural responses. c) Estimation of the ensemble activity given stimuli through our framework. d) Statistical coupling between ensemble activity and stimuli derived from 1D spatial statistic frameworks [23] [25] [24].

Figure 10 :

 10 Figure 10: Ensemble overlapping and co-activity extraction. (left) Number of time windows where neuronal ensembles are simultaneously coactive. (right) Number of neurons belonging to several communities.

Figure 11 :

 11 Figure 11: Benchmarking of the overlapping clustering method for different statistical simulation parameters. Default parameters: A = 4, N = 400, T = 1000, p = [.1, .1, .1, .1], α = [.15, .15, .15, .15], λ = [λ 0 , λ 1 , λ ≥2 ] = [0.05, 0.8, 1]. 20

Figure V. 3 :

 3 Figure V.3: Deconvolved neuronal spiking matrix using CD foopsi from [134]

Figure V. 4 :

 4 Figure V.4: Denoising steps to considered in Hydra. a) Correction of Contraction Burst pollution by averaging the intensity over the entire segmented animal or through tracks to remove the common component of CB activity ; b) Neuropil contamination corrected by removing the neighborhood activities outside somas ; c) Correction of z-axis motion and "cross-talk" by extracting via ICA common underlying patterns in neighboring patches.

Figure V. 5 :

 5 Figure V.5: Spatial statistic framework adapted to infer statistical coupling between stimulation and ensemble activity. a) (left) the CSR condition allow to formulate a null hypothesis of the expected number of points in a region through a random spatial distribution. (middle) When the meshgrid is fitted on a first kind of reference points we can estimate the statistical coupling between two sets of points. In the 2D space, a search shape can be used to measure the expected number of points under CSR and the effective number. It provides a point-process statistical coupling called SODA. [247] (right) A simplification of the research can be performed using level-set function instead of concentric rings. It provides another version called GSoda [248] [249] ; b) (left) An adaptation of this point process problem could be formulate in a 1-dimensional case estimating the coupling between two sequences of time dates. (middle) the research shape correspond to 1-dimensional rings i.e pairs of segments (right) a similar simplification can be performed by creating non-overlapping non-overflowing shapes ; c) By fitting the 1-dimensional shapes as time segments it is possible measure the statistical coupling between a time point sequence and a time interval. The time point sequence is the ensemble activity. The time interval corresponding to the time interval of a specific kind of stimulations. The statistical coupling framework based on spatial statistics can be modified to account such point process phenomenon. Adapted from[START_REF] Kubler | Statistical Coupling Between time Point-Processes[END_REF]

6 :

 6 Scanning of the visual receptive field V1 by straight stimuli Estimated ensemble activities on real data Neuronal ensembles detected in mice brain visual cortex. a) Static visual gratings are displayed in front of an immobilized mice. b) Neuronal ensemble activity matrix is estimated and superimposed on stimulation. c) The statistical coupling between neuronal ensemble activity and stimulation is calculated via point process theory.

Figure V. 7 :

 7 Figure V.7: Mapping of neuronal ensembles and extraction of statistical estimates. a) Mapping of the neuronal ensemble in the brain visual cortex. Frame to frame tracking of neuronal activation regarding their ensemble can be performed. b) Number of times neuronal ensembles are coactive. c) Number of times neurons belong to several ensembles.

Figure V. 8 :Figure V. 9 :

 89 Figure V.8: Random initialization allows the Gibbs Sampler to converge on synthetic data.

123

 123

Figure A. 1 :Figure A. 2 :

 12 Figure A.1: Dynamics of a the neuronal system analyze though avalanche statistics. a) (left) The sandpile model analogy summarizes the emergence of the dynamics in a system. Each time a new grain is added there is a probability to move neighborhood grains. When the number of grain is enough high, avalanche cascade effects can happen and displace a huge number of grain in the same time. The grain corresponds to neuron electric firing. (right) A cascade event is group of firing of neurons surrounded by silence. A number of channels is recruited corresponding to a number of diverse elements involved in the process. b). (left) The dynamic system can converge toward different activation regime. It can switch to a balance regime of activation, a high level of activity regime or activity can simply vanish and stay marginal in the dynamical system. (right) This regime convergence can be quantified using the avalanche ratio criterion that summarizes the importance of the avalanches in the dynamics. c). The dynamic theory predicts a power-law evolution of the dynamics of the neural systems. d) These expected power law distribution about avalanche statistical features are demonstrated via synthetic data. e) These theories are confirmed on in vitro neural networks where such power law distribution are experimentally observed. f). The same observations are highlighted in vivo neural network like at Hydra Vulgaris. (Adapted from [255] [81] [259] [235] [49])

Figure A. 3 :

 3 Figure A.3: Small-Worldness emergence for in vitro neural networks adapted from [275]. a) Often observed in in vitro experiments the degree of smallworldness of the networks increase due to an increase of the clustering coefficient given an average link path remaining constant over development. b) Node degree scales as a power law distribution in accordance with [81].

Figure A. 4 :

 4 Figure A.4: Small-Worldness emergence in vivo neural networks at Hydra Vulgaris adapted from [49]

Figure B. 1 :

 1 Figure B.1:Step of the Louvain algorithm from[START_REF] Blondel | Fast unfolding of communities in large networks[END_REF]

Pmrr rr exp -n 2 r w rr 2

 2 (A, g|γ, w, k) ∝ P (g|γ, k)P (A|g, w) r<s w mrs rs exp -nrnswrs (B.3)

R

  b = Z T b P Z b withP = ww T (B.9)

N

  i = P AS C i , such that P AS = i P C i P T C i = P C P T C (B.10)

Figure B. 2 :

 2 Figure B.2: Community detection via SVD from [108] [126]. a) Vectorization of the time bin in the spike matrix. Each time step is seens as a population vector of activated neurons. b)The matrix is normalized using a TF-IDF transform from the NLP (Natural Language Processing) community. It limits the influence of too activate cells that are not statistically significant. c) A similarity matrix is calculated using a pairwise cosine similarity. d) The matrix is binarized and decomposed using a SVD. It unveils orthogonal components to highlight dissimilar population that repeat over times. e) The decomposition is perfomed in an orthogonal basis. f) An estimation of an ensemble activity matrix is performed and matched with the mice visual stimulations. g) It highlights frame by frame the activated population in response to stimulations.

All 3

 3 fundamental model principles lead to a statistical model whose likelihood function can be calculated.P (t, w, s|θ = (n, p, λ)) = ( p µ ) 1-w kµ )( i (w kt i ) s ik (1-λ t i (w kt i )) 1-s ik ) (B.13)

Figure B. 4 : 3 )

 43 Figure B.4: Non-exhaustive literature Tree of Community Detection techniques

  From neural networks to neuronal ensembles . . . . . . . . . . . . . . . . . . i A new neural doctrine based on neuronal ensembles . . . . . . . . . . ii Neuronal ensembles embed randomness . . . . . . . . . . . . . . . . . iii Neuronal ensembles are functionally overlapping . . . . . . . . . . . . From patch clamp to calcium and voltage indicators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Imaging neurons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Tracking calcium activity . . . . . . . . . . . . . . . . . . . . . . . . .

	III Chasing neuronal ensembles that encode behaviors
	1	Imaging and tracking single neurons in behaving animals . . . . . . . . . . . .
		i	Neural activity monitoring :

Introduction 1 Neuronal network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . i From Brain to neuronal network . . . . . . . . . . . . . . . . . . . . . ii Imaging neuronal network in vivo . . . . . . . . . . . . . . . . . . . . . 2 From single neurons to population coding . . . . . . . . . . . . . . . . . . . . i Neurons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Conveying information between neurons . . . . . . . . . . . . . . . . . iii From single neurons to neural networks . . . . . . . . . . . . . . . . . 3 iv Detection of neuronal ensembles . . . . . . . . . . . . . . . . . . . . . 2 Neuronal ensembles encode behaviors . . . . . . . . . . . . . . . . . . . . . . . i Behavior tracking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Behavior prediction . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  𝜔 𝑠 𝑖 = ⋃ Ω 𝑡 𝑠. 𝑡. {∃ 𝑡 𝑙 ∈ 𝑠 𝑖 |𝑡 𝑙 ∈ Ω 𝑡±𝑡𝑜𝑙 } 𝑠 𝑗 ) = (𝐾 𝑠 𝑖 (𝑠 𝑗 ) -𝑚 𝑠 𝑖 ,𝑠 𝑗 ) /𝜎 𝑠 𝑖 ,𝑠 𝑗

	4. DAT -AUTOMATIC THRESHOLDING FOR SPIKE ESTIMATION
	}} 𝑗 = |𝑐 𝑖 Ω 𝑡 𝑇 Ω 𝑢 } Κ 𝑠 𝑖 (𝑠 𝑗 ) 𝑡-𝑡𝑜𝑙≤𝑢≤𝑡+𝑡𝑜𝑙 𝛽 = 𝜎 𝑠 𝑖 ,𝑠 𝑗 [𝐾 ̃𝑠𝑖 𝟏 𝐾 ̃𝑠𝑖 >0 ] (𝑠 𝑖 ) 𝑠 𝑖 𝑠 𝑗 𝛿 𝑡±𝑡𝑜𝑙 (𝑠 𝑖 ) = 𝟏{∃ 𝑡 𝑙 ∈ 𝑠 𝑖 |𝑡 𝑙 ∈ ⋃ 𝑠 𝑖 𝑠 𝑖 𝑡𝑜𝑙 = 2 1≤𝑡≤𝜏 𝑁 𝑖 (𝑠 𝑗 ) 𝑗 𝑁 𝑠 𝑖 (𝑠 𝑗 ) = ∑ 𝛿 𝜔 𝑠 𝑖 (𝑡 𝑗 𝑘 ) 𝑁 𝑗 𝑘=1 𝑁 𝑠 𝑖 (𝑠 𝑗 ) Κ 𝑠 𝑖 (𝑠 𝑗 ) = 𝜏 𝑁 𝑗 𝑁 𝑠 𝑖 (𝑠 𝑗 ) 𝑠 𝑗 Ω 𝜌[𝑠 𝑖 , 𝑠 𝑗 ] = . 𝜎 𝑠 𝑖 ,𝑠 𝑗 𝛽 [𝐾 ̃𝑠𝑖 𝟏 𝐾 ̃𝑠𝑖 >0 ] (𝑠 𝑗 ) , 1 𝐹 𝑠𝑐𝑜𝑟𝑒 (𝑠 𝑖 , 𝑠 𝑖 Ω = {Ω 𝑡 } 1≤ 𝑡≤τ |𝜔 𝑠 𝑖 | = ∑ 𝛿 𝑡±𝑡𝑜𝑙 (𝑠 𝑖 ) 𝟏{∃ 𝑡 𝑙 ∈ 𝑠 𝑖 |𝑡 𝑙 ∈ Ω 𝑡 } 𝑠 𝑖 𝜏 𝑡=1 𝛿 𝜔 𝑠 𝑖 (𝑡 𝑗 𝑘 ) = 𝟏{𝑡 𝑗 𝑘 ∈ 𝜔 𝑠 𝑖 } (𝑡 𝑗 𝑘 ) 1≤𝑘≤𝑁 𝑗 𝑚 𝑠 𝑖 ,𝑠 𝑗 = |𝜔 𝑠 𝑖 | 𝜎 𝑠 𝑖 ,𝑠 𝑗 2 = |𝜔 𝑠 𝑖 |(𝜏-|𝜔 𝑠 𝑖 |) 𝑁 𝑗 𝑠 𝑖 𝑠 𝑗 𝐾 ̃𝑠𝑖 (𝑁(0,1) 𝑠 𝑗 𝜌 0 𝑠 𝑖 𝑠 𝑖 * * ) = 2 × 𝑇𝑃(𝑠 𝑖 , 𝑠 𝑖 𝑠 𝑖 𝑠 𝑗 𝑠 𝑖 𝑠 𝑖 = 𝑐 𝑖 𝑗 + 𝑟 𝑖 𝑐 𝑖 𝑗 𝑠 𝑗 𝑟 𝑖 𝑠 𝑗 𝐶 𝑖 𝑗 | 𝑠 𝑖 * (𝑇) 𝑠 𝑖 * (𝑇) = 𝑐 𝑖 𝑗 (𝑇) + 𝑟 𝑖 (𝑇) + 𝑚 𝑖 (𝑇) 𝑚 𝑖 𝑐 𝑖 𝑗 (𝑇) = {𝑐 𝑖 𝑗 > 𝑇}, 𝑟 𝑖 (𝑇) = 𝛿 𝑡 (𝑠 𝑖 ) = 𝜔 𝑠 𝑖 𝑖 𝑠 𝑖 𝑅 𝑖 = |𝑟 𝑖 𝑗 | {𝑟 𝑖 > 𝑇} and 𝑚 𝑖 (𝑇) = {𝑚 𝑖 > 𝑇} 𝐶 𝑖 𝑗 (0) = 𝐶 𝑖 𝑗 , 𝑅 𝑖 (0) = 𝑅 𝑖 𝑀 𝑖 (0) ≫ 𝐶 𝑖 𝑗 + 𝑅 𝑖 𝑇𝑃(𝑠 𝑖 , 𝑠 𝑖 * ) = (𝐶 𝑖 𝑗 (𝑇) + 𝑅 𝑖 (𝑇)) , 𝐹𝑃(𝑠 𝑖 , 𝑠 𝑖 * ) = 𝑀 𝑖 (𝑇), 𝐹𝑁(𝑠 𝑖 , 𝑠 𝑖 * ) = (𝐶 𝑖 𝑗 (0) + 𝑅 𝑖 (0) -𝐶 𝑖 𝑗 (𝑇) -𝑅 𝑖 (𝑇)) 𝐹 -𝑠𝑐𝑜𝑟𝑒 𝐹 𝑠𝑐𝑜𝑟𝑒 (𝑠 𝑖 , 𝑠 𝑖 * (𝑇)) = 2 (𝐶 𝑖 𝑗 (𝑇) + 𝑅 𝑖 (𝑇)) 2 (𝐶 𝑖 𝑗 (𝑇) + 𝑅 𝑖 (𝑇)) + 𝑀 𝑖 (𝑇) + (𝐶 𝑖 𝑗 (0) + 𝑅 𝑖 (0) -𝐶 𝑖 𝑗 (𝑇) -𝑅 𝑖 (𝑇)) = 2 (𝐶 𝑖 𝑗 (𝑇) + 𝑅 𝑖 (𝑇)) 𝐶 𝑖 𝑗 (𝑇) + 𝑅 𝑖 (𝑇) + 𝑀 𝑖 (𝑇) + 𝐶 𝑖 𝑗 + 𝑅 𝑖 𝐶 𝑖 𝑗 (𝑇) = 𝐶 𝑖 𝑗 × 𝑓 𝑖 𝑐 (𝑇), 𝑅 𝑖 (𝑇) = 𝑅 𝑖 × 𝑓 𝑖 𝑟 (𝑇) and 𝑀 𝑖 (𝑇) = 𝑀 𝑖 × 𝑔 𝑖 (𝑇) 𝑓 𝑖 𝑐 (𝑇), 𝑓 𝑖 𝑟 (𝑇) 𝑔 𝑖 (𝑇) 𝑇 𝑓 𝑖 𝑐 (0), = 𝑓 𝑖 𝑟 (0) = 𝑔 𝑖 (0) = 1 lim 𝑇→∞ 𝑓 𝑖 𝑐 (𝑇), 𝑓 𝑖 𝑟 (𝑇), 𝑔 𝑖 (𝑇) = 0 . 𝑓 𝑖 𝑐 (𝑇) ≈ 𝑓 𝑖 𝑟 (𝑇) = 𝑓 𝑖 (𝑇), ∀𝑇 > 0, 𝑔 𝑖 (𝑇) < 𝑓 𝑖 (𝑇) 𝐹 𝑠𝑐𝑜𝑟𝑒 (𝑠 𝑖 , 𝑠 𝑖 * (𝑇)) ≈ 2(𝐶 𝑖 𝑗 + 𝑅 𝑖 )𝑓 𝑖 (𝑇) (𝐶 𝑖 * (0)) ≈ 2(𝐶 𝑖 𝑗 +𝑅 𝑖 ) 2(𝐶 𝑖 𝑗 +𝑅 𝑖 )+𝑀 𝑖 ≈ 0 𝑀 𝑖 ≫ 𝐶 𝑖 𝑗 , 𝑅 𝑖 lim 𝑇→∞ 𝐹 -𝑠𝑐𝑜𝑟𝑒(𝑠 𝑖 , 𝑠 𝑖 * (𝑇)) = 0 𝑓 𝑖 (𝑇) 𝑔 𝑖 (𝑇) 𝑑𝐹 𝑠𝑐𝑜𝑟𝑒 (𝑠 𝑖 , 𝑠 𝑖 * (𝑇)) 𝑑𝑇 = 2(𝐶 𝑖 𝑗 + 𝑅 𝑖 )[𝑀 𝑖 [𝑓 𝑖 ′ 𝑔 𝑖 -𝑓 𝑖 𝑔 𝑖 ′ ](𝑇) + (𝐶 𝑖 𝑗 + 𝑅 𝑖 )𝑓 𝑖 ′ (𝑇)] ((𝐶 𝑖 𝑗 + 𝑅 𝑖 )(1 + 𝑓 𝑖 (𝑇)) + 𝑀 𝑖 𝑔 𝑖 (𝑇)) 2 𝑇~0 0 𝑓 𝑖 ′ (0) ≈ 0 > 𝑔 𝑖 ′ (0)) 𝑑𝐹 𝑠𝑐𝑜𝑟𝑒 (𝑠 𝑖 , 𝑠 𝑖 * (𝑇)) 𝑑𝑇 (𝑇 = 0) ≈ 2(𝐶 𝑖 𝑗 + 𝑅 𝑖 ) 𝑀 𝑖 [𝑓 𝑖 ′ -𝑔 𝑖 ′ ](0) > 0 𝑇 > 0 𝑇 * > 0 0 𝑇 → ∞ 𝑇 * = argmax 𝑇 𝐹 𝑠𝑐𝑜𝑟𝑒 (𝑠 𝑖 , 𝑠 𝑖 * (𝑇)) 𝑑𝐹 𝑠𝑐𝑜𝑟𝑒 (𝑠 𝑖 ,𝑠 𝑖 * (𝑇 * )) 𝑑𝑇 = 0 𝑀 𝑖 [𝑓 𝑖 ′ 𝑔 𝑖 -𝑓 𝑖 𝑔 𝑖 ′ ](𝑇 * ) + (𝐶 𝑖 𝑗 + 𝑅 𝑖 )𝑓 𝑖 ′ (𝑇 * ) = 0 𝑀 𝑖 ≫ 𝐶 𝑖 𝑗 + 𝑅 𝑖 𝑇 * 𝑇 ̃ * 𝑇 ̃ * [𝑓 𝑖 ′ 𝑔 𝑖 -𝑓 𝑖 𝑔 𝑖 ′ ](𝑇 ̃ * ) = 0 𝜖 ≪ 𝑀 𝑖 𝑇 * 𝜌 𝜖 [𝑠 𝑗 * (𝑇 𝑗 ), 𝑠 𝑖 * (𝑇 * )] 𝑇 𝑗 𝜖 (𝑓 𝑖 ′ 𝑔 𝑖 -𝑓 𝑖 𝑔 𝑖 ′ )[𝑇 * ] = 0 𝑇 𝑖 * * 𝑇 𝑗 𝑠 𝑖 * 𝑠 𝑗 𝒊, 𝒋 𝒔 𝒊 , 𝒔 𝒋 𝒔 𝒋 𝒔 𝒊 𝟎 ≤ 𝜶 ≤ 𝟏 𝒔 𝒊 𝝀 𝟎 𝒔 𝒋 (𝟏 -𝜶)𝝀 𝟎 𝜶𝑵 𝒋 = 𝜶𝝀 𝟎 𝑻 𝒔 𝒊 𝒔(𝒕) 𝒄(𝒕) 𝒌(𝒕) 𝒄(𝒕) = 𝒔(𝒕) * 𝒌(𝒕) 𝒌(𝒕) = 𝑨 𝒆 -( 𝒕 𝝉 𝑫 ) 𝜷 𝟏 -𝒆 -𝒕-𝝁 𝝉 𝑹 𝑗 + 𝑅 𝐹 -𝑠𝑐𝑜𝑟𝑒(𝑠 𝑖 , 𝑠 𝑖 𝑑𝑇 [𝑇] = 𝝉 𝑫 𝜷 > 𝟎 𝝁 𝜏ℎ(𝑇 𝑗 )𝐶 𝑖 𝛽((𝐶 𝑖 𝑗 + 𝑅 𝑖 )𝑓 𝑖 (𝑇) + 𝑀 𝑖 𝑔 𝑖 (𝑇) + 𝜖) 2 [𝑀 𝑖 (𝑓 𝑖 ′ 𝑔 𝑖 -𝑓 𝑖 𝑔 𝑖 ′ ) + 𝜖𝑓 𝑖 ′ ](𝑇) 𝝉 𝑹

𝑖 )(1 + 𝑓 𝑖 (𝑇)) + 𝑀 𝑖 𝑔 𝑖 (𝑇) 0 < 𝑇 * < ∞ 𝑇 * = argmax 𝑇 𝐹 𝑠𝑐𝑜𝑟𝑒 (𝑠 𝑖 , 𝑠 𝑖 * (𝑇)).

  J→I = P (I t+1 , i

	(k) t , j t ) log 2 (k)	P [i t+1 |i (k) t , j t ] (l) P [i t+1 |i (k) t ]	(C.13)

NEURONAL ENSEMBLES ENCODE BEHAVIORS

EXTRACTING NEURONAL ENSEMBLES IN CALCIUM IMAGING

A VERSATILE SIMULATOR FOR SPIKE INFERENCE TECHNIQUE THAT EMBEDS HYDRA VULGARIS NEURAL ACTIVITY PROPERTIES

SODA -A SPATIAL STATISTIC FRAMEWORK FOR NEURON COUPLING ESTIMATION AND BEHAVIOR PREDICTION

3. SODA -A SPATIAL STATISTIC FRAMEWORK FOR NEURON COUPLING ESTIMATION AND BEHAVIOR PREDICTION

DAT -AUTOMATIC THRESHOLDING FOR SPIKE ESTIMATION

CHAPTER IV. PUBLICATIONS

5. BINOE -BAYESIAN INFERENCE FRAMEWORK FOR NEURONAL OVERLAPPING ENSEMBLES

BINOE -BAYESIAN INFERENCE FRAMEWORK FOR NEURONAL OVERLAPPING ENSEMBLES

Acknowledgements

Poisson Process

Kernel Convolution

Poisson-Gaussian Imaging Noise

Non-linear Baseline

taxonomies depending on the communities using them. A transdisciplinary unification of methods and theories is currently lacking and would greatly facilitate future research. Another paradox we encounter today is the existence of extremely advanced analysis methods which are rarely or never used in practice by non-mathematicians due to a lack of clear, illustrated explanations, exhaustive benchmarking, and open-source availability code.

The implementation of statistical methods based on clear model biologically inspired provides interpretable statistical features and appears as the cornerstone of the understanding of neural computation inherently driven by randomness. "As far as he can achieve it, readability is as important for the scientific writer as it is for the novelist". Donald O. Hebb

Appendix A

Dynamics and topology shape neural networks through neuronal ensembles

It is well-known that stimulation of single neurons can affect population activity in in vitro and in vivo experimental conditions. A single neuron can directly influence network dynamics which motivates recent works about optogenetics, heterogeneous receptor expression and retrogradely transported viral vectors. The targeting of neuro-anatomical structure of the microcircuits activated during behavioral tasks is at stake. Such structures evolve during development and highlight common dynamical and topological patterns. For a long time recording 2 neurons with two patch clamps, injecting pre-synaptic current to measure the post-synaptic response was the gold-standard technique to measure the evolution of the coupling efficacy between individual cells. The monitoring of the simultaneous activity of almost all single cells in a living animal or inside an engineered in vitro neuron network provides new perspectives and experiments to highlight how neural network topological configuration can be shaped relatively to the emergence of dynamics.

The study of a neural network is based on two cornerstones inherently linked: Dynamics and Structure.

Dynamics theory

Dynamics is, usually, studied via avalanche size, duration, shape, branching ratio and complexity [START_REF] Heiney | Neuronal avalanche dynamics and functional connectivity elucidate information propagation in vitro[END_REF]. An neural avalanche is a set of recording time bins that shows sequential neural activity surrounded by time bins without neuron spikes. Its size is defined by the number of spikes emitted during an event and the number of different channels (= neurons) involved. The statistical analysis of avalanche size and duration distribution are supposed to provide information about the dynamics of the networks. The quantities follow power law distributions whose parameters define the dynamics of the system and the converge toward sub-and super-critical dynamic states [START_REF] Hesse | Self-organized criticality as a fundamental property of neural systems[END_REF]. The NCC toolbox can be used to extract all these information and classify the dynamic regimes of neuron network recording [START_REF] Marshall | Analysis of Power Laws, Shape Collapses, and Neural Complexity: New Techniques and MATLAB Support via the NCC Toolbox[END_REF]. From an avalanche, the shape can be extracted as a time series corresponding to its size over each time frame. Under specific conditions, the shapes can be rescaled to reveal similar profile. The branching ratio is a mean ratio of the number of "descendant nodes" recalled during avalanche events. It provides a metric showing the relevance in the activation of the neuron and the recall of similar patterns in neural networks. The complexity [START_REF] Tononi | A measure for brain complexity: relating functional segregation and integration in the nervous system[END_REF] of a system is the intrinsic mechanisms of the neural code to allow distant sites to communicate with each other by optimizing the wiring cost. For the last 10 years lots of studies have worked on trying to understand the development of neural networks but only few of them tried to directly interact with it, providing clear insights. Understanding how topological network properties evolve by stimulating or killing neuron cells via Optogenetics would be a promising perspective to unravel the underlying mechanisms of the neural circuitry. Moreover, understanding that topology and dynamics shape jointly during a neural network aggregation provides insights on the tools to detect the neural ensembles or search for some peculiar structures.

Appendix B

Detection of neuronal ensembles techniques

Graph Theory Techniques applied to in vivo neural networks

The adaptation of graph theory approaches to in vivo neural networks relies on the fidelity of a graph representation and the equivalence of neuronal ensembles to network communities in the graph despite all user-defined choices: functional connectivity metric, threshold of connectivity matrix, definition of a community in relation to a graph criterion to optimize for instance. The main issues of the direct application of more general frameworks not directly related to biological features is to lead to erroneous conclusions due to data complexity or the absence of conclusions supported by statistical significance that could handle to the peculiarities of the underlying problem. In the following section, the main significant method applied in the in vivo network assembly analysis is reviewed.

Modularity Optimization using Louvain clustering algorithm

The most widely used technique to detect communities is by maximizing the graph modularity creating clusters of nodes densely connected compared to a random null model. The modularity Q quantifies the quality of any graph partition. The most famous definition of modularity is given by the Newman-Girvan modularity [START_REF] Girvan | Community structure in social and biological networks[END_REF] definition :

A well-known implementation to optimize the modularity is the stochastic two steps greedy implementation Louvain clustering [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF]. Generally, such techniques modularitybased, requires to formulate a graph network null model of random graph. Different null models exist and the comparison between the properties of observed network graph and the expected one under H0 provides insights to remove non-statistically significant edges and detect communities. It is worth noting that a huge effort has been made by the mathematical community to explore the space of null connectivity models. Null models demonstrate statistical distributions at random graph about spectral quantities such as the Tracy-Widom distribution or Marchenko-Pastur distribution but these are advanced mathematical theories barely used in practice and whose implications and bias are not fully mastered for real-world networks.

The algorithm has been applied by J. Lovas 2021 on Hydra Vulgaris to extract neural 15)

In the more advanced model, the number of communities itself can be seen as statistical parameter to estimate and dynamically change in the sampling routine as metropolis-hasting switching regimes. Then, the main algorithm is benchmarked against some state-of-the-art techniques firstly on synthetic data and then, on experimental data.

Appendix C

Graph theories techniques rely on functional connectivity metrics

Functional connectivity highlights the statistical coupling between two distant areas or nodes in their activity. Two sites demonstrate a strong FC if their respective activities are co-dependent, meaning that a behavior or a cognitive function is encoded by a remote interaction between these distant sites. Generally, this measure is performed into 2 steps: the first is to define a coupling metric that will be applied on time-series signals. The second is to apply this metric between all pairs of nodes to generate an adjacency matrix, also called connectivity matrix, that will summarize all existing edges (functional connections) between each pair of nodes (neuron cells). This matrix is the cornerstone of the FC since its structure will reveal the network topology expected (see : previous part) of the neuron network. Another core of community detection is the creation of a null hypothesis of node activity, outside behavior or cognitive task, to detect a statistically significant increase of node cell interactions during an experiment. This step is also a core of the community detection research. For cells, some existing functional connectivity metrics can be directly adapted and applied on fluorescence traces, but it seems that direct application on spike trains are often considered too. It is worth noting that some techniques evaluate the correlation between two signals inside the time domain trying to find linear relationships for instance in the evolution of the time series signals. Other techniques work in the frequency domain, especially for theories where coupling is extracted from oscillatory phase locking. The choice of the FC metric to use depends on the underlying hypothesis that is being tested. The reviews are from [START_REF] Wang | Overlapping community detection in dynamic networks[END_REF] [242] [START_REF] Bastos | A Tutorial Review of Functional Connectivity Analysis Methods and Their Interpretational Pitfalls[END_REF].

Here is a non-exhaustive list of well-known functional connectivity metrics. 1) Cross-correlation: Applied on time series or point process events corresponding to spike trains. It measures the direct influence of a neuron on another reference one [START_REF] Salinas | Correlated neuronal activity and the flow of neural information[END_REF] [START_REF] Heiney | Neuronal avalanche dynamics and functional connectivity elucidate information propagation in vitro[END_REF]. The metric does not handle temporal structure of data nor directionality. 2) Cross-covariance: It is the probability to obtain a spike for a first signal at time s,