Muhammad Tayyab
email: muhammad.tayyab@imag.fr

Yves Usson
email: yves.usson@imag.fr

Thierry Léveillard
email: thierry.leveillard@st-antoine.inserm.fr

Jacques Demongeot
email: jacques.demongeot@imag.fr

Yannick Fouquet

Nicolas Vuillerme

A Lontos

E Promayon

Understanding Physiological and Degenerative Natural Vision Mechanisms to Define Contrast and Contour Operators

Keywords: neural networks, lateral inhibition, retinitis pigmentosa, image processing, Markov random fields cell contouring, cell counting, gastrulation, biomechanical model, streak formation, invagination simulation The living organisms are very complexpart digital and part analogy mechanisms

Advancement in Image Acquisition Equipment and the progress made in the Image Processing Methods have led mathematicians and computer scientists into areas which are critical for physicians and biologists. Early diagnoses of diseases (such as blindness, cancer and digestive problems) are priority areas of interest in medicine and it is largely based on examination of microscopic images. The acquisition of these images was facilitated by the considerable advancement in the development of a number of new equipment, like Two-Photon Laser Microscope, Fluorescent Microscope, etc. which already provide a good idea of the interesting features of the object being viewed. However, some biomedical images are not suitable for fast and easy extraction of relevant information. Image Processing methods have been providing good support to deduce useful information about the objects of interest in these biological images. Fast computational methods allow complete analysis, in a very short time, of a series of images, providing a reasonably good idea about the desired characteristics. The thesis focuses on the application of these methods in 3 independent series of images for 3 different types of diagnosis or inference. First, Images of retina were treated for the detection of cones, in case of total or partial loss of rods, , in the pathology of retinitis pigmentosa (retinitis pigmentosa, or RP). The method used was able to detect and count the number of cones in each frame, and to estimate their average distance, which are important parameters in the diagnosis of degenerative type of RP. Secondly, a gastrulation process in drosophila was studied to observe the movements and cell mitosis, and results were consistent with recent research. Finally, another series of images was treated where biological cells were observed to undergo mitosis. The source was a video from a photon laser microscope. In this video, objects of interest were biological cells. The idea was to track the cells if they undergo mitosis. Cell position, spacing and sometimes contour of the cell membrane are broadly the factors limiting the accuracy during the segmentation process. Appropriate method of image enhancement and segmentation were chosen to develop a computational method to observe this mitosis, while proposing human intervention in some cases to eliminate any false inference.

interesting field of medical imaging and enabling me to work in this challenging domain. It has been a wonderful experience for me. I will always be thankful to him for his moral support whenever I was looking for it. I really acknowledge and appreciate the time, he dedicated to me throughout the course of the thesis. The knowledge I acquired from him is invaluable for me.

I would like to thank Yves Usson of TIMC Laboratory for assisting me in statistical analysis of some results. I also owe thanks to the personnel of AGIM and TIMC laboratories for providing me a pleasant working environment especially Olivier Hansen, Céline Fontant, Nadine Perraud, Caroline Dos Santos and Maribel Chenin for their administrative support. I wish to communicate my thanks to the University J. Fourier of Grenoble (UJF) and to the Centre National de la Recherche Scientifique (CNRS) for allowing me to work in their laboratories.

I wish to extend my thanks to all my family members for their patience (and sometimes lack of it) and prayers (and I knew they were always there for me).

v Résumé Les recherches avancées en acquisition d'image et les progrès effectués dans les méthodes de traitement d'image ont conduit les mathématiciens et informaticiens dans des domaines d'application cruciaux pour les médecins et les biologistes. Le diagnostic précoce de maladies (comme la cécité, le cancer et les problèmes digestifs) sont des domaines d'intérêt prioritaire en médecine et il repose en grande partie sur l'examen d'images microscopiques, dont l'acquisition a été favorisée par le développement considérable de nouveaux équipements, comme le microscope bi-photonique à balayage laser et le microscope de fluorescence par réflexion totale interne, qui fournissent déjà une bonne idée des caractéristiques intéressantes de l'objet observé. Cependant, certaines images biomédicales ne sont pas appropriées { l'extraction rapide et aisée d'une quantité suffisante et pertinente d'information. Des méthodes de traitement d'image adaptées fournissent alors un soutien nécessaire { l'extraction d'informations utiles (en particulier au diagnostic) sur les objets d'intérêt dans ces images biomédicales. Des méthodes de calcul rapides permettent l'analyse complète, dans un temps très court, d'une série d'images, offrant une bonne estimation des caractéristiques souhaitées. La présente thèse porte sur l'application de ces méthodes à trois séries d'images destinées à trois différents types de diagnostic ou d'inférence. Tout d'abord, des images de rétine ont été traités pour la détection des cônes, en cas de disparition partielle ou totale des bâtonnets, dans la pathologie de la rétinite pigmentaire (RP). Le logiciel utilisé a été capable de détecter et de compter le nombre de cônes dans chaque image, puis de donner une estimation de leur distance moyenne, paramètres permettant le diagnostic de dégénérescence de type RP. En second lieu, le processus de gastrulation chez la drosophile a été étudié, afin d'observer les mouvements cellulaires et les mitoses éventuelles et les résultats sont cohérents avec les recherches récentes. Pour finir, une autre série d'images a été traitée, où la source est une imagerie vidéo, à partir d'un microscope bi-photonique à balayage laser. Dans cette vidéo, les objets d'intérêt sont des cellules biologiques. L'objectif est de suivre les cellules et d'observer si elles subissent une mitose. La position de la cellule, l'espacement intercellulaire et parfois le contour de la membrane cellulaire sont globalement des facteurs limitant la précision dans cette vidéo. Une méthode appropriée d'amélioration de l'image et de segmentation a donc été choisie, permettant de Human science, ecological science and many other fields of biology are interlinked and have important information for each other. Bioinformatics emerged to integrate all these fields accelerating the research and development in the field of biology and medicine with the help of computer technology.

After years of research and development, biologists and physicians manage to collect huge amount of valuable data from their experiments all over the world. These biologists and physicians need a tool that can link information from different areas like biology, stati stics, genomics etc…, in order to make their research faster. For instance, they may need some data regarding effects of particular gene or protein on human being and their effect on experimental animal, so that they can interlink and generate some beneficial results or antidote that helps in human development or defence against pathogens. Eventually, bioinformatics provides that help in interlinking information from different fields or species and leads to quick new results.

Digitizing the data is also an important aspect in bioinformatics, so that information could be transmitted wherever required. The format and size of data is an important factor in this regard. Storage and compilation of huge biological data so that it could be used by professional researchers in the field of medicine and biology is the primary objective of bioinformatics. Storage of raw data in the form of databases might not be a difficult task, but extraction of 'only' required information from the raw data is a real challenge for researchers. This challenge led scientists and engineers into development of bioinformatics tools dedicated for biological applications. These tools are designed so as to help the biologists, with little knowledge of computers software, to analyze their raw data without any difficulty in using these tools.

MOTIVATION

A very important application of bioinformatics lies in the early diagnosis of fatal diseases as cancer and other diseases leading to permanent disability of any human organ. Early diagnosis, sometimes, decreases the chances of falling victim of these diseases. This aspect is a prime motivating factor for biologists as well as computer scientists.

In view of above, we discuss in section II, processing of retinal images from an eye suffering from retinitis pigmentosa in which abnormalities of the photoreceptor cells (rods and cones) lead to progressive visual loss. Different phases for processing images of retinal explants are explained from scaling and contrasting to segmentation, inference and diagnosis. Histogram-thresholding technique has been used for segmentation of images in which decision about a pixel to be part of an object, or of its background or of its boundary, is determined from the histograms of the local image intensities in that region. After segmentation of cones, their geometric centres are calculated in order to determine the horizontal intercept distances between the cones. The main objective of this work is to identify the areas of macular degeneration where cones have disappeared. Research publications with supplementary illustration of this work have been included in Appendices A and B.

In section III, we analyze by using image processing techniques the cell deformation and motion in the Drosophila melanogaster embryo searching to delimit the first period of invagination without proliferation of the gastrulation process and then propose a biomechanical model, based only on the consideration of elastic and contractile forces exerted on cell walls and on the centrosome through the combination of myosin contraction and cytoskeleton rigidity. Numerical simulations of this model made during the period of gastrulation without proliferation suggest that the model adequately simulates in-vivo cell behaviour, showing the start of the streak formation at the two extremities of the embryo cylinder, followed by a propagation of the invagination to its central part.

Appendices C and D provide further illustration of the work in this section.

In section IV we explain how this work could be extended to different other fields within the wide domain of biology and medicine.

RETINA AS A SYSTEM

The rationale in choosing retina as a system to explore is based on both biological and medical reasons. Retina is an accessible part of the central nervous system (Fig. 1.1(a), [2]).

It is accessible to biologists because of the laminar organization of different interconnected neurons of the outer and inner retina, and to physiologists that can test the function of the system. Retina is also one of the favorite systems for developmental and molecular biologists as the molecular mechanism of photoreceptor differentiation in the fly Drosophila melanogaster (most commonly used model organism) has been discovered

with the elucidation of the RAS/MAP kinase signaling pathway. The genetic of inherited retinal degenerations (IRD) has been a field of intense research for the last 20 years, with more than 100 genes identified as cause of these diseases ( [3]). For our society, the retina is of interest to medicine since the occurrence of a peculiar form of IRD, the aged-related macular degeneration (AMD) is increasing as our population ages. IRD will be amenable to gene therapy in the near future as dogs have been already effectively treated ( [4]).

From these many aspects of retinal biology, our system will deal with one of the most promising approach for therapy of retinitis pigmentosa (RP) a common form of IRD with sequential loss to rod and cone photoreceptors. The laboratory Inserm U592 led by J. Sahel has pioneered the protection of cones as an effective approach to treat RP ( [6]). A novel tropic factor Rod-derived Cone Viability Factor (RdCVF) was shown to be necessary for cones to survive after rod death in a mouse model of the disease. This therapeutic approach would alleviate a common mechanism leading to blindness in several different conditions that initially affect rods, independent of the nature of the causative mutation. We intend to model the cellular interactions between rods and cones for better defining the biological properties of RdCVF signaling and for treatment of the disease.

Finally, additional reasons to select retina as a system and develop an image processing from the wealth of data generated in the frame of RdCVF research come from the functional genomic approaches: proteomic ( [7]); transcriptomic ( [8]) and images from the expression cloning platforms of the Inserm U592 ( [6]).

RETINITIS PIGMENTOSA AND MACULAR DEGENERATION

Retinitis Pigmentosa (RP) is a group of hereditary diseases of retina of the eye. Retina is a light sensitive delicate layer of cells which picks up the pictures and transmits them to the brain. In humans, there are two types of light sensitive cells in the retina: rod cells and cone cells. RP may be caused by a breakdown in the function of the rods or the cones in some part of the retina. The retina is so complex that breakdowns may occur in a variety of ways and so RP is not a single disorder. The breakdown of cone function may be called Macular Degeneration.

The macula is located roughly in the center of retina, temporal to the optic nerve. It is a small and highly sensitive part of the retina responsible for detailed central vision. The macula allows us to appreciate detail and perform tasks that require central vision such as reading. Macular degeneration is a general term used to describe a number of diseases of the retina. Of these diseases, age-related macular degeneration is the most common, primarily affecting people over the age of 60. With macular degeneration, a spot called a scotoma may appear in the central vision. This scotoma may look light, or dark, or blurred. Fovea is the very center of macula. It measures about 1.5 mm in diameter equivalent to 5 degree field of view represented in Figure 1-2. In Figure 1-3, the retinal fovea (A) measures 1.5 mm in diameter and contains only 50,000 cones on a total of 5 millions cones in the whole retina. The retinal periphery (B) contains 120 millions rods. There is a blind spot in the field of vision where the optic nerve leads back into the brain (Figure 1-2 Left). There are no rods and cones at this point, so a small object becomes invisible if it lies in this field of vision. In Figure 1-4, C and D represent 3D shadow projections of confocal images of the fluorescent cones in fixed retinal whole mounts.

SPATIAL MODELING OF RETINAL DEGENERATION

Our animal model of rod-cone interactions will be the rd1 mouse, a retinitis pigmentosa (RP) model that carries a mutation in a photoreceptor-specific rod, the phosphodiesterase gene (Pde6b), leading to sequential degeneration of rod and cone photoreceptors. In this model, rod photoreceptors degenerate through apoptosis between the first and fifth week post-natal, the cones die subsequently over a period that extends to three months ( [12]) (see Figure 12345). In that study the evolution of both cell populations has been studied at some limited time points. A more precise description of cell population evolution at different post-natal days will be obtained through image analysis involving cone counting together with the estimation of the cone density and homogeneity over the entire surface of the retina.

The recent identification of the Rod-derived Cone Viability Factor (RdCVF) ( [6]) has significantly increased the knowledge about the interaction between rods and cones and has contributed to the understanding of secondary cone degeneration and visual loss in RP.

This trophic factor was identified by expression cloning using a high content screening approach. It is mainly produced by rod cells and is necessary for cone survival. After rod loss, cones die by lack of trophic support. The restoration of RdCVF expression in patients after rods have been lost is a current therapeutic treatment of RP ( [12]).

The aim of the retinal image study is to model the degeneration of cones from a spatial point of view together with images of flat-mounted retina of the rd1 mouse (Figure 12345)

and gene expression data for the same samples by Affymetrix Genechip hybridization (PN8, 9,10, 11, to PN35). The cones are identified by their ability to bind specifically to the peanut agglutinin ([14)]. The U592 INSERM platform of cell counting has been used by T.

Léveillard to generate images of the retinal explants of oriented (nasal-temporal) rd1 mice at post-natal days corresponding to the loss of both rods and cones. The orientation will be used to look for a gradient in the process of cone degeneration as it has been previously

reported ([15]). An objective will also be to complete the image series by analyzing the secondary cone degeneration (weekly from 5 weeks to 3 months). Raw data in the mice retina already represent a huge number of images (about 500 GB)

showing retinas affected by the retinal degeneration. These images, recorded with confocal microscopy at 9 levels of depth inside the neural retina, allow us to detect in the whole retina the positions of cones.

STATISTICAL TEST ON THE UNIFORM DISTRIBUTION OF CONES

To model the spatial degeneration of the retina it is necessary to assess if this degeneration is uniform or oriented in some direction. A gradient of degeneration from the periphery to the center is indeed suspected by some experts of the field. A common reference for the acquisition and representation of the retinal images is defined (based on the naso-temporal axis and on the blind spot position). This reference will permit registration of the images between individuals and for the same individual the in-depth quality control between neuro-retina slices. The determination of a radial grid will allow us to obtain the best framework for counting the cones present in each region and the best scale is the largest one, corresponding to the occurrence of at most one cell in each region. Beyond this uniformity test, a test will be done to detect a possible gradient of the presence of cones in the naso-temporal axis or in the centrifugal direction (due to heterogeneity in the distribution of rod loss, responsible for the secondary cone degeneration).

A mathematical model considers the influence of rods on cones on a semi-lattice Model.

This model is an Ising like statistical structure ( [16][17]). Cones are assumed to be at the nodes of a two-dimensional lattice whereas rods fill the interspaces. Cones are supposed to be either alive or dead. The reason of this modeling is that cones can be observed (see Figure 1-6), whereas rods can not be. Rods will be modeled in a relative way, considering the proportion of living rods compared to the initial rods population. This proportion will be assessed through the rod-specific gene expression (rhodopsin, rod-arresting). For the rd1 mouse retina, the transcriptomic temporal data are available over 11 time points from day 1 to day 35 post-natal.

When the gene expression is maximum, we consider that 100% of rods are present in the interspaces. If the gene expression value is halved, we also consider that the presence of rods is halved. The lattice can be re-indexed as:

𝑥 = (𝑥 1 … 𝑥 𝑛) )
If n cones are present, each variable 𝑥 𝑖 (𝑖 = 1 … 𝑛) is equal to 1 if the corresponding cone i is alive, and to 0 otherwise. The probabilities that a state changes, are given by 𝑃 𝑥 𝑖 𝑡 + 1 = 1|𝑥(𝑡) and 𝑃 𝑥 𝑖 𝑡 + 1 = 0|𝑥(𝑡) . In this model, it is assumed that:

𝑃 𝑥 𝑖 𝑡 + 1 = 1|𝑥(𝑡) = 𝑓(𝑘)
where k(t) is the number of rods alive in a certain range from i (a disk of a certain radius for example) and f a function increasing from 0 to a, with 𝑎 ≤ 1. Different functions can be tried such as the affine function, sigmoid, etc. The model parameters will be chosen so as to maximize its maximum likelihood.

This model is built under the hypothesis that cones and rods are uniformly distributed (hypothesis reinforced if the test above has proved the homogeneity of the occurrence of still alive cones) and exert positive interactions on the cones (through a trophic substance, RdCVF, secreted by rods and diffusing in the neural retina to the cones). However, the model should be changed in case of an orientation of the rods and cones distribution. The model chosen in the present study is an Ising-like or logistic model, with :

𝑓(𝑘) = 𝑒 𝑤𝑘 1 + 𝑒 𝑤𝑘 = 1 1 + 𝑒 -𝑤𝑘 or 𝑤 = log log 1 -𝑓(𝑘) 𝑓(𝑘) 𝑘
where w is the rod-cone interaction parameter. In Figure 1234567; ki is calculated for example in the lattice box i of the SE quadrant by counting the number of rods in its 8 neighbour boxes: this empirical number equals 3. It is easy to

𝑤 𝑆𝐸 = - log 1 -𝑓(𝑘) 𝑓 (
show that the two empirical meanswSE andwNW are significantly different (p=0.01), by using a Student t-test which allows concluding to the heterogeneity of the distribution of cones and rods between the quadrants NW and SE.

SHAPE MODELING DURING CONE DIFFERENTIATION

Transplantation of retinal sheet in the rd1 retina restores visual response attributable to cone function [19] and RdCVF is able to prevent cone degeneration [6]. Nevertheless, RdCVF was identified by expression cloning using cell viability, not cell function, as assay. It is theoretically possible that a secreted differentiation factor not yet identified will complement the viability activity of RdCVF and that both will be required for effective In these retinal cultures from chicken embryos, the morphology of viable cells is not uniform (Figure 1234567). The bipolar cells are cones, the long cell at the top left represents a glial cell and the perfectly round cells are most likely not differentiated. When seeded at low density and in the absence of serum, retinal precursor cells adopt a cone differentiation state by a default differentiation pathway resulting in about 60 to 80% cone-enriched culture ( [20]).

We have developed a statistical method to define the different cell morphologies. The features used to classify the morphotypes are based on classical measurements such as cell perimeter/area ratio, elliptic ratio, etc. Next, a mixture model clustering is performed to identify all the different morphotypes. At least two are present: one for the living, but nonfunctional cones and the other for the differentiated retinal progenitor cells. A further study could involve a more detailed modeling. For instance, for star-convex cells it is possible to represent cells by angular functions, and use distances between these functions to cluster the corresponding cells. In this process one must pay attention to have a method that is invariant under rotation. Even further we can model the cells by using all their surface intensity values descriptions and not only their contours. This would involve using functions with a two-dimensional support.

Having established the list of morphotypes that we will assume to represent different differentiation states, essentially rods and cones in the present study, we have looked for their distribution in the expression cloning experiment. The plan is to look for statistical enrichment of the morphotypes for certain pools of clones. If such aggregation is found, the interpretation will be that in these pools, a secreted factor will be driving the differentiation of retinal precursor to a peculiar differentiation fate, as seen by its morphotype. The selected pools are divided in sub-pools of 10 clones; the assay is repeated by using a limited dilution the cDNA encoding the differentiation factor.

To validate the factor identified, the protein or a viral vector encoding for this factor has 1.6.

IMAGE SEGMENTATION

Enhancement algorithms have been used to reduce image noise and increase the contrast of desired objects. The idea is to extract a wealth of information from an image, which in our case would be to detect the abnormalities in a medical image. Accurate interpretation may become difficult in case noise levels are relatively high.

Human Vision System has the capability to group the image into segments which contain pixels with some common characteristic (same grey level or colour). The goal of image segmentation is to detect and extract the regions which compose an image. Common requirements from a segmentation process suggest that it should be unsupervised, computationally efficient, flexible and robust. In many cases, the segmentation method decides the outcome of the entire analysis, since subsequent process is based on the segmented region(s) ( [21][22][23][24][25][26]).

Simple segmentation techniques include Thresholding and Clustering. In Thresholding technique, thresholds could be static and may be applied globally across the image; or dynamic which could be applied locally so that the threshold varies across the image (thresholds could be chosen using bi-modal and multimodal histograms). Not necessarily such thresholds would exist in an image or it could be difficult to infer an appropriate threshold from the image histogram. In places where noise and texture varies, this technique becomes inadequate for segmentation. In Clustering Technique, a collection of data points that have similar colour or texture are clustered together. This could be done in hierarchical manner or by partitioning the image into various disjoint "imagettes".

Clustering and Thresholding methods are global and do not retain positional information (invariant to spatial rearrangement of pixels). The resulting segments are not connected and could be widely scattered.

The segmentation techniques could also be classified as Region or Boundary-based. In Region-based Segmentation, Region Growing and Split-and-Merge methods exist, whereas in Boundary-based Segmentation, First-Order and Second-Order approaches are common.

In First Order methods, a gradient mask is convolved with the image to obtain the gradient vector. Sobel Operator introduced a weighing local average measure on both ends of the central pixel. Canny proposed an analogue of the lateral inhibition retinal function, the derivative-of-Gaussian (d-o-G) filter as a near optimal filter with respect to three edgefinding criteria, i.e., a good localisation of the edge, one response to one edge, and high probability of detecting true edge points and low probability of falsely detecting non-edge points. Deriche, based on Canny's criteria, used a filter with impulse response similar to that of d-o-G but which lends itself to direct implementation as a recursive filter.

In Second Order methods, optimal edges are found by searching for places where the second derivative is zero, or where zero crossing occurs. Later these boundary elements are connected to form lines or curves (Hough Transform and hierarchical Hough Transform).

In the following chapter, we will focus on a necessary step before the modeling, the image analysis.

CHAPTER 2

IMAGE ANALYSIS

IMAGE ACQUISITION AND FIRST TREATMENT

The images of retina in mice have been obtained after dissection and 2D-mapping from confocal microscopy of retinal fragments. Here we will use only the first of the 9 levels of the confocal slices and we will assume that we observe only a typical imagette at the peripheral level on which we detect a transition to the cone degeneracy (consecutive to the primary rods disease, which causes the disappearance of the trophic factor secreted by rods ensuring the cones survival). Treatment done on the imagette has to be repeated all along retinal mapping, but the methodology is shown on this only one. Figure 2-1 shows the first contrasting step done on the imagette before the stereological study.

STEREOLOGY APPROACH

Statisticians regard stereology as a form of sampling theory for spatial populations.

Stereology exploits the fact that some 3-D quantities can be determined without 3-D reconstruction. For example, the 3-D volume of any object can be determined from the 2-D areas of its plane sections, without reconstructing the object (this means that stereology only works for certain quantities, like volume, and not for other quantities). We assume that the altitudes of the cones, orthogonal to the retinal surfaces, are following Gaussian distribution with respect to a base plane tangent to the exterior of the retina. Hence, the distribution of the z-values of the cone-sections is supposed to be also Gaussian and same is true for the surface areas of the cone-sections.

SAMPLING PRINCIPLES ([1])

In order to extrapolate from a few plane sections to the 3-D material, essentially the sections must be 'typical' or 'representative' ones. In classical stereology, we assume that any plane section is typical. The material is supposed to be homogeneous. A statistical 3-D model of the material is built on these two assumptions. In modern stereology (designbased), we select the plane sections according to a random sampling protocol, by choosing a random position to start cutting the material. This method is effective for nonhomogeneous materials as in biomedical sciences especially for lungs, kidneys, bones and brain.

The intercept method is based on counting the number of segments determined within a set of objects cut by a set of parallel scan lines along a series of directions ( [21][22][23][24]). The advantage of the method is that it still gives a statistical estimation of the grain size distribution even when individual contours are only partially outlined. This is particularly interesting as it matches the real operation conditions. A prerequisite however, is that the subset of contours revealed by etching is a dense and unbiased sample (e.g., no preferential etching due to crystallographic orientation).

THE INTERCEPT METHOD

The intercept method is a stereological tool defined on binary images. It is based on counting segments from objects of interest cut by scan lines along a series of parallel directions ( [21][22][23][24]). This means that the preliminary image-thresholding-operations, to be performed, are totally independent.

Hence, the binary intercept method uses simpler images and simpler algorithms, but it has some disadvantages. The result is not robust with respect to the grey-level threshold especially using absolute thresholding on images with non-uniform luminance. As a result, some boundaries are not detected whereas many spurious ones do appear. Image rotation has to be performed using an adequate interpolation algorithm which may lead to artificial intercept counts. These are required to be filtered out by applying additional processing such as the linear filters [25].

The grey-level intercept method deals directly with the original grey-level image that contains more information than the binary-level image. This method analyzes the greylevel profile or its derivative. Rotation is performed on the original grey -level image using a linear re-sampling algorithm that better preserves the information. The transitions along any direction can be selected with sub-pixel accuracy.

DETECTION OF THE GREY-LEVEL TRANSITIONS

Grey-level transitions can be analyzed at the sub-pixel scale using the tools implemented into the EasyGauge library (Euresys). A scan line is drawn in a portion of image with its associated profile and derivative (Figure 2-2). Peaks along the profile are identified to determine the transitions positions. A peak is the area comprised between the profile (or its derivative) and a horizontal user-defined threshold level. All pixel values within the peak are taken into account to compute the exact transition location. The parameters which help to detect valuable transitions include Threshold level, Amplitude of the transition, Minimal area of the transition, and Thickness (linked to the number of parallel segments used to extract the data profile). The three last parameters reduce the influence of noise. Transition detection can be improved by applying a Gaussian filter along the profile by removing the noise and keeping only real transitions.

Two types of transitions can be relevant. 'White to black to white' transitions correspond to light objects with darker boundaries. In this case, the transitions are detected along the profile. 'Black to white' or 'white to black' transitions correspond to dark objects in a light matrix or vice-versa. They are detected along the derivative profile. The pixels are not treated individually anymore but the relative difference between adjacent pixels is considered.

IMPLEMENTATION OF THE METHOD

A validation program can be implemented by following the three steps given below: The 5x5 mask is obtained using a local grey-histogram thresholding:

• we calculate the local histogram giving the grey level in a pixel and in each of its 24 neighbors,

• we decide a pixel is internal (and respectively external or frontier) if this local histogram is located below (and respectively above or with a part below and the complementary above) a certain grey threshold,

• a pixel is not assigned to one of these 3 types, if its local histogram is not verifying any of the 3 threshold rules. CHAPTER 3

RESULTS

As mentioned earlier, we are considering images of retina in mice obtained after dissection and 2D-mapping from confocal microscopy of retinal fragments. Currently we intend to process only the first of the 9 levels of the confocal slices and we assume to observe only a typical imagette at the peripheral level and on that imagette, we observe a transition to the cone degeneracy (consecutive to the primary rods disease, which causes the disappearance of the trophic factor secreted by rods ensuring the cones survival).

Figure 3-1 shows the image from retinal explants. This image has been divided into imagettes each of dimensions 650x515 pixels. For the current task, these imagettes have been treated to be independent whereas they are actually dependent. The image mainly depicts macular degeneration, and it contains mostly cones in all regions because of the cell pathologic process observed in retinitis pigmentosa. A typical imagette is shown in Figure 3-2.

PROCESSING

IDENTIFICATION OF CONES

Having said that the images chosen are from an eye with retinitis pigmentosa, our objects of interest would be essentially cones. We intend to identify the cones and then explain their different geometries, their horizontal intercept distances, and if possible find their volumes, surface area and perimeters.

CONTRAST ENHANCEMENT

Images having low contrast may result from inadequate illumination, wrong lens aperture settings, or insufficient dynamic range in the image sensor. Contrast enhancement increases the dynamic range of the image. In contrast stretching, the existing grey-levels of the original image are mapped to new levels. The image in Figure 3-2 is a 16-bit TIFF format image corresponding to 65536 grey levels suffering from low contrast. We intend to scale it down to only 256 levels. Further to enhance the total contrast of the image, we set the lowest pixel value in the imagette to zero (black) and the highest pixel value to our new scale of 255. Grey-levels between the upper and lower bounds are set to a linear ramp of values between 0 and 255. Edge detection reduces the amount of data and filters out useless information, while preserving the important structural properties in an image. Gradient and Laplacian are two widely used methods in edge detection. The gradient method detects the edges by looking for the maximum and minimum in the first derivative of the image. The Laplacian method searches for zero crossings in the second derivative of the image to find edges. The derivative represents the maximum change located at the center of an edge in the image.

Sobel filter is a typical gradient filter. A pixel location is declared an edge location if the value of the gradient exceeds some threshold. As mentioned before, edges will have higher pixel intensity values than those surrounding it. So, once a threshold is set, we can compare the gradient value to the threshold value and detect an edge whenever the threshold is exceeded. Furthermore, when the first derivative is at a maximum, the second derivative is zero. As a result, another alternative to finding the location of an edge is to locate the zeros in the second derivative. This is known as the Laplacian method.

SOBEL

The Sobel operator performs a 2-D spatial gradient measurement on an image. Typically, it is used to find the approximate absolute gradient magnitude at each point in an input grayscale image. The simplest Sobel edge detector uses a pair of 3x3 convolution masks, one estimating the gradient in the x-direction (columns) and the other estimating the gradient in the y-direction (rows). A convolution mask is usually much smaller than the actual image. As a result, the mask is slid over the image, manipulating a square of pixels at a time. Sobel Masks are shown in Figure 3-3.

LAPLACIAN

Similarly, Laplacian of the image detects intensity changes in an image and hence is often After the image has been enhanced, the next task is the segmentation of the objects of interest (cones). The resultant image would be binary consisting of existing well-defined cones in the image and the background. The image can be converted to binary by choosing a threshold value after viewing the histogram of the imagette. To accomplish this task, the concept of histogram thresholding is employed.

HISTOGRAM THRESHOLDING

Thresholding is based on the assumption that the objects can be separated from the background according to the grey-level values. Grey-levels within a defined range are selected as belonging to the foreground whereas grey-levels outside that range are considered as the background. As it is clear from histogram of scaled image (Figure 34567), it is difficult to segment this image with a global threshold, whereas selection of an appropriate threshold is also a difficult task. We intend to pick a mask of 5 x 5, find its histogram, check whether all pixel values are low or high or otherwise, following the procedure given below:

• If all the pixel values in the 5 x 5 mask lie in the right half, then the central pixel belongs to the inside of the object;

• If all the pixel values in the 5 x 5 mask lie in the left half, then the central pixel belongs to the background;

• If the histogram of the mask contains both low and high values, then the central pixel is the one corresponding to the edge of the object. -the mean intercone distance (for both barycentric distances on Figure 3-12 and intercept distance on Table 3-1) is about 3.5 times the maximal cone diameter (10 µm, after Figure 1-

3) while it equals between 2 and 5 times this diameter on regular cones (Figures 1-3 and 1-

6),

-the distribution of the objects diameters is about the same as in the slices of Figure 123.

The number of sections of small diameter (less than 5 µm) is equal to 23, that is about twice the number of sections (11) of big diameter (more than 10 µm); these numbers are respectively 15 and 7 in Figure 123.

Hence, we can consider the cones distribution in RI quadrant as practically normal (despite the small number of rods, which can be represented by the small cells).

The situation is opposite in case of quadrant LS in which the mean intercone distance is about 4.5 times (Table . 3-1) the maximal cone diameter and where the big cell objects are absent. Thus quadrant LS presents a pathologic cone distribution typical of retina degeneracy. The situation in the two remaining quadrants RS and LI is an intermediary between LS and RI, which shows that we have chosen to study an imagette of transition between cones degeneracy and cones survival zones. A fuzzy frontier of the degeneracy domain could be obtained by chaining these transition imagettes from the whole retina image of Figure 2-1(a).

After repeating this procedure on all imagettes of the retinal explants, one can identify the regions of macular degeneration, where the cones have disappeared. With the knowledge we have about cone densities in fovea and its periphery, we are in a position to infer the information about the different zones. For instance, we have the image for the retinal explants at day 1, which can be a reference for comparing the results of the image at day 7.

Each imagette comprises a sick region, a healthy region or a combination of two.

Practically, all sick portions in individual imagettes would combine to give a complete picture of the sick region in the complete image. Thus, boundaries of the sick region can be identified by a simple segmentation technique which would be marked along the junctions of the healthy and the sick region. will permit the image registration between individuals and for the same individual, and the quality control in depth between neuro-retina slices, 2. determine a grid scale for data representation with renormalization control and counting. This step will allow us to obtain the best framework for counting, at each vertex of the grid, the present cones or rods and the best scale would be the largest corresponding to the occurrence of at most one cell on each vertex, 3. perform homogeneity test for the occurrence of the still alive cones (test of spatial uniformity for the localizations of the gravity centers of the cones, see . This test will be done to detect a possible gradient for the presence of cones in the nasotemporal axis or in the centrifugal direction (due to heterogeneity in the distribution of rods responsible of the secondary cones degeneracy), 4. test the interaction of cones with a population of dying rods whose activating influence diminishes. This is a test on the significance of strictly positive values of interaction coefficients in an Ising like statistical structure ( [16]), by supposing that interactions are exerted in the framework of an Ising isotropic model, in which the interactions range, sign and intensity can be estimated (using the fact that the Gibbsian distribution of the Ising model with constant interaction weight leads naturally to an exponential structure). This test will be performed with a zero hypothesis in which rods were uniformly distributed (hypothesis reinforced if the step 3 above has proved the homogeneity of the still alive cones occurrence) and exerted positive interaction on the cones (through a trophic substance secreted by rods and diffusing in the neural retina up to the cones). The test will be performed for different values of the density characteristics of the uniform distribution of the rods. This step will also allow the estimation (of the maximum likelihood, convergent and without bias) of the interaction coefficients in the frame of an Ising like exponential structure.

5. test the interaction between cones. The still alive cones can interact between themselves by emitting inhibitors or activators and this possibility will be tested by using the same procedure as in phase 4. CHICKEN 6. characterization of morphotypes of cones from the shape of their cellular membrane.

The features used for classify the morphotypes will be based on classical measurements like nucleo-cytoplasmic index, cell perimeter/area ratio,…, followed by a non-hierarchical clustering procedure (like k-means method) searching at least for two morphotypes, one for the alive but no more functional cones and the other for the still functional cones, with a possible third population representing dedifferentiated neuro-retinal cells.

7. counting the sizes of the morphotypes populations at different days during the retinal maturation for revisiting the classical curve of cones survival.

In future, an application could be possible on the human retina, if a non-invasive procedure (like the fluorescence of non-toxic membrane potential dependent optical dyes) could be done on the still functional retinal cells giving, during the degeneracy, the localization of the still alive and functional cells. gastrula. Some of ectoderm cells (orange) move inward, forming the endoderm (red) [33].

Invagination is the morphogenetic process by which an embryo takes form and is the initial step of gastrulation, the massive reorganization of the embryo from a simple spherical ball of cells, the blastula, into a multi-layered organism, with differentiated germ layers: endoderm, mesoderm, and ectoderm (Figure 1-1). More localized invaginations also occur later in embryonic development, to form coelom, etc. Invagination is the formation of a cleavage furrow during cytokinesis in animal cells.

During the early stages of development, the primitive streak is the structure that will establish bilateral symmetry, determine the site of gastrulation and initiate germ layer formation. To form the streak, reptiles, birds and mammals arrange mesenchymal cells along the prospective midline, establishing the first embryonic axis, as well as the place where cells will ingress and migrate during the process of gastrulation and germ layer formation. The primitive streak extends through this midline (Figure 1-2) and creates the antero-posterior body axis, becoming the second symmetry-breaking event in the embryo (the first being the spherical symmetry breaking due to the fecundation and the second the revolution symmetry breaking due to the streak formation), and marks the beginning of gastrulation. This process involves the ingression they will differentiate into the three germ layers (endoderm, mesoderm, ectoderm) that will give rise to all the tissues of the adult organism. Koller's sickle is a local thickening of cells that acts as a margin separating sheets of cells from posterior margin of blastoderm as they migrate before to push primary hypoblast cells anteriorly to form a secondary hypoblast known as the endoblast.

The Posterior Marginal Zone (PMZ) of Koller's sickle can induce a primitive streak (Figure 1-2) and Hensen's node, acting as an organizer.

STATE OF ART

During the morphogenesis of the gastrula -the second step of the embryogenesis after the blastula stage (Figure 1-1), we observe, in the majority of the animal development dynamics, the following phenomena: i) cell motion is partly guided by chemotaxis, in order to supply their substrate demand, and also to respect the epigenetic architecture ruled by morphogens, ii) cell shape is due to a constriction controlled by cell differentiation, and iii) final gastric tube is obtained from cell proliferation relaxing the forces exerted on the cell plasmic membrane and optimizing the cellular "nutritive Area / inner

Volume" (A/V) ratio.

Concerning the differentiation process in the gastrulation process, some cells of the embryo start to take the shape of a bottle (bottle or flask cells), decreasing the surface at the interface with their nourishment substrate in mammal (like chicken in Figure 1-3) and insect (like Drosophila in Figure 1-4) embryos. At the end of the gastrulation, these bottle cells start to divide and grow, increasing their A/V ratio. In this section, i) using image processing techniques, we try to follow, through pure cell motion, the first phase of the gastrulation, causing the invagination without any proliferation, and ii) we try to observe that the cell constriction results in a streak starting from the two extremities of the embryo and propagating subsequently towards the middle portion.

In order to simulate the initiation of this phenomenon, authors in [37][38] Cells could also move in response to chemo-attractants/repellents. The cells in Koller's sickle could aggregate towards the dorsal midline in response to a chemo-attractant that these cells themselves produce, resulting in their aggregation at the ventral meeting point.

The cells at the tip of the forming streak could then acquire the ability to respond to a gradient of another chemo-attractant, for instance the one that all cells in the epiblast produced, resulting in their movement towards the midline of the embryo.

IMAGE PROCESSING TECHNIQUES: CELL CONTOURING AND COUNTING

Localised or oriented cell divisions could locally expand the cell mass in the epiblast, displacing cells at other locations. However, the initial cell flows preceding streak formation occurred in the absence of cell division ( [34]).

Many cell-based research studies require the counting of cells prior to beginning an experiment. Estimation of cell density in various regions of embryo is thus an integral part of such studies. Profile counts or stereological techniques could be used to have an estimate of the cell density in a particular region. In regions where cell density increases enormously, care must be taken, if they occur in regions of interest for the embryo morphogenesis, i.e., zones with cell deformation and/or cell proliferation. In any case, image enhancement techniques need to be applied on the images to get well-defined cell contours.

The available images ( [3]) depict actin-myosin networks in colour. Here, we are focussing more on cell boundaries in particular regions, in order to have an idea of cell density and shape, as well as to observe whether cells divide in those regions. After evaluating the colour histogram of the image, we filter out the noise by selecting an appropriate colour threshold value. In the absence of noise, cell boundaries are easily visualized and hence cells could be tracked.

Cell boundaries are further enhanced by applying contrast algorithms ( [35]) on the colourfiltered image again by using simple threshold values after studying the gray-histogram.

Cell proliferation results in increase in volume of the mother cell, thereby causing daughter cells to recurrently contract then stabilize. Older cells are more vulnerable to shrink themselves as a result of forces from the neighbouring cells, due to proliferation at the other end. Cell density, as expected, would increase enormously in such regions with respect to other regions (Figure 1-3), which could lead to invagination, closer to the point of maximum cell density.

Observing the colour histogram of the original image (Figure 2-2 (bottom)) suggests we could start enhancement by filtering out the coloured objects which are not contributing to the actual information in the image. Since we are interested only in the cell boundaries which are in red, implying we could filter out all colours other than red, which include mostly the green pixels in the image.

After colour filtering, we observe the intensity image (gray-level) and observe its histogram. Using Sobel edge detection algorithm, we get the boundaries of the cells.

Segmentation sometimes need be done in an iterative manner, where results are not sufficient after one treatment to deduce useful info from the image. Applying the edge filter, the resulting binary image observed might have discontinuities in cell boundaries. These discontinuities can be easily removed by applying a dilation technique, followed by some erosion technique to get well defined contours. Once the boundaries are well defined, all statistical methods could be applied, making it easy to automatically count the number of cells in that frame. Calculating the number of pixels in each cell gives an estimate of the mean size of cells, thus important parameters like cell density, cell shape, etc could be calculated.

RESULTS OBTAINED FROM CELL CONTOURING AND COUNTING

Available series of images were treated to observe the phenomenon of streak formation.

The idea was to observe the cell movement and proliferation in the region where a constriction appeared at a later stage. The source of data was in the form of a video (in AVI format) from which frames (BMP images) have to be extracted at a reasonable frame rate.

Thus a series of data images was obtained (resolution: 510 x 255 pixels). After following the image processing steps explained in 2.1., we were able to have an exact visual idea of how the cells move and how cell density changes in the region where invagination occurs experimentally, and also observe if there is any division of cell in that region before the development of constriction.

We started from the observing the colour histogram of the image (Figure 2

THE EMBRYO ARCHITECTURE

The embryo cylindrical cephalo-caudal structure (Figure 3-1) has two zones of cylindricity breaking, the two hemispheric extremities and the cells at the boundary between the two geometries presents a curvature due to a relaxation of the internal rigid properties of their cytoskeleton. The appearance of the gastrulation streak occurs at their level and later propagates to the central cylindrical part of the embryo. Taking this into account, the differential contractility of these boundary cells constitutes the model core and has been implemented over a 3D mesh representing the external embryo surface ( [37], [38], [42]).

THE MODEL (A. LONTOS)

Each cell evolves following the following mechanical laws:

1. The sum of the orthogonal forces exerted on the cell walls is equal to its mass (proportional to its inner area) times its orthogonal (to the wall) acceleration. The external forces are the resultant of the extracellular (Pext) and intracellular (Pint) pressure forces, i.e., the contact force exerted by the neighbouring cells (sharing a part or the whole of the wall) exerted by other cells and the internal force due to the cell pressure exerted via the cyto-skeleton and the plasmic membrane ( [7]). Each force is equal to a coefficient (e.g., the physical pressure) times the length of the wall on which it is exerted.

2. If we suppose that the initial cell configuration is in an equilibrium state, we calculate an admissible set of parameters values respecting this equilibrium.

3. Then, we leave the cell system evolve depending on the energetic balance ruling the cytoskeleton apical polymerization ( [4][5][6]) controlled by a specific genetic regulatory network comprising essentially concertina (cta), actin, myosin, Rho and RhoGEF genes ( [7][8][9][10]), by choosing a small time step, by updating sequentially each cell and by calculating their displacements respecting the no-overlapping rule. At each step we update the cell common walls by supposing that cell contacts are close, ensured between cells by cadherins and gap junctions ( [11][12]), and with the extracellular matrix by integrins and adhesins. Cell motion involves a change in its inner area: we suppose that growth occurs where internal forces P int are larger than external ones P ext (cell has to be stretched), as for constrained growth in continuous media mechanics.

When the external forces dominate, the cell is supposed to be compressible and can be constricted. The cells divide longitudinally or laterally when their ratio perimeter over area is too small ( [13][14][15][16][17]).

The revolution symmetry breaking consists in making two cylinders inside the embryo cylinder, one coming from the ventral furrow and giving the digestive tube (this phase is called gastrulation) and the other coming from the dorsal furrow and giving the neural tube (this phase is called neurulation). The first gastrulation step consists of an apical concentration over fluctuation of one of the cytoskeleton components (myosin, actin, tubulin,...) or one of the enzymes (ADenylateKinase -ADK -or Nucleoside Diphosphate Kinase -NDK) or one of the energy molecules (ATP, GTP) involved [18][19][20][21][22][23].

This apical overconcentration diffuses and reaches the extremities [11,12] of the cylindrical portion of the embryo at which gap junctions are less important with the cells of the "hemispheric" terminations. Then two first bottle cells can appear at the two extremities of the diffusion line, then this first contraction can propagate until the centre of the cylindrical part as noticed in [11,12] during the 6 first hours of the gastrulation in Drosophila melanogaster. After this first phase of the ventral furrow formation, for regularizing the tensegrity forces, based on a synergy between balanced tension and compression components of cells, first divisions occur from extremities to centre by applying the first Thom's cell law. We can simulate such a process on a transversal slice of the embryo and reproduce the invagination of the gastrulation from both the contraction and consecutively the proliferation phases. In future, spatial and temporal resolution of the data could be increased by sampling in 3D.

Spatial resolution can be improved using confocal or biphotonic microscopy. Temporal resolution can be improved using the high-frequency cine-imaging techniques. Raman or Infrared Spectroscopy could enhance the functional resolution of the data. These techniques will allow to detect, precisely, the frontiers between the successive phases of the morphogenetic processes which cause the symmetry breaking in the embryo, i.e., understand the causes of its disease. Cell tracking can be useful in the treatment of lots of diseases like cardiac cell based therapy [1] or gene therapy [2].

Viability of stem cells used to develop effective cell therapies needs to be evaluated in a noninvasive manner [3]. Although optical imaging methods allow tracking cells in real time [4], there could be a requirement to store the images and analyze later the characteristics of individual cells. In a relatively wide series of images, it is practically impossible to manually analyze all images and infer useful information from them. Automation of a series of thousands of such images could help biologists and physicians to save a lot of time. The limitations of such a method would be decided by the presence of noise in the images, continuous change in cell position, size and shape, and also spacing among the cells. In view of these limitations various algorithms are applied stepwise to address each issue. Noise could be reduced by using contrasting and common noise removal techniques. Once we have considerably removed the noise from the images, next task would be to segment the cells from the background for each frame. After segmentation has been done successfully, the idea is to relate two consecutive frames for all segmented cells. In an ideal scenario, we should be able to find the same cells in the next frame, unless they have divided or disappeared. For this, the position of the centroid of a cell in consecutive frames is one of the criteria on the basis of which we could decide whether we are tracking the same cell. If we see no cell in the vicinity of '𝑑' pixel radius from the previous position of centroid, we could decide that the cell has disappeared. The value of '𝑑' could be chosen by analyzing the speed of cell movement in some frames for a series of images. After segmentation and cell counting, an increase in cell density in a frame would depict an occurrence of cell division in that frame, assuming no cells entered from the exterior of the frame.

PROCESS OF MITOSIS

In Mitosis a replicated genetic material in a parent cell is equally distributed into two daughter cells. The method proposed here involves different steps to track cells individually in a series of frames and also to detect when it divides. Available images were extracted from a recorded video of early mitotic cell divisions in a drosophila embryo.

IMAGE ENHANCEMENT AND NOISE REMOVAL

The primary task in the proposed method is to choose an appropriate method to remove unwanted objects. To separate the objects from background, the task gets difficult where the objects and background have the same grey level, because simple threshold techniques can not be applied directly. Even more difficult it becomes if we have different grey levels in different sections of the images. One way of tackling this issue is to divide the image into blocks (say 100x100 pixels) or imagettes and apply local thresholds instead of one global threshold for the entire frame. A good threshold value could be selected from individual histograms of those blocks.

After we have sufficiently enhanced the image (cf. Section II), we need to remove the noise (and more generally any unnecessary object). From a typical frame, we could observe and decide ourselves which size and shape of cell we are expecting to segment in that frame.

Based on that, any object not following that criteria won't be regarded as a cell. For instant we chose to ignore objects below 16 pixels in area. So again applying this threshold, we remove all objects we don't intend to detect as cells.

SEGMENTATION OF CELLS

After noise removal, we have the image ready to be segmented. We already could observe objects, which are expected to be our objects of interest, cells in our case. Segmentation techniques have been discussed in Sections II and III.

CELL TRACKING

Once cells have been segmented successfully, we can apply the tracking algorithm. Based on the position of any segmented cell, our intention is to follow its movement and observe whether it divides, and then follow the daughter cells. The position of the cell at current time instance is updated based on current and previous observations. Still, errors could be caused by background clutter and multi-object confusion [6]. Observing the data images, we can estimate the speed of motion of cells, and based on that we can estimate the size of the window beyond which the cell won't move from one frame to another. Placing this window at the centroid of the segmented cell, the algorithm need be run only on this window in the next frame to find the same cell. The position of centroid would be updated in each frame and hence position of the window would be updated. The method will detect the biggest segment in this area, unless there a division occurs whereby it will look for two cells. From Figure 1-1, we observe that the position of centroid of the cell changes, by a considerable distance, whenever a division actually occurs. Hence, the method looks for a second object in the vicinity at that point. The two objects could be labeled and tracked Each frame size is 360x480 pixels. One isolated cell is easily segmented, whose movement will be followed in the subsequent frames.

CELL TRACKING

The method depicted in the block diagram of Whenever '𝑑' increases to a value greater than a predefined value, we look for another segment in the vicinity of the current cell.

SOME INTERMEDIATE IMAGES

Although the method tracks the cells in most of the cases, there have been some cases where it was not effective. If the limitations could be addressed, we could have a very fast method to track cells in all frames and we could observe cell division or cell death which could be very useful for biologists and physicians.

I. INTRODUCTION

In the vertebrate retina, cones are hyperpolarized when illuminated by light, but also receive a depolarizing input when receptors some distance away are illuminated. This antagonistic center-surround response is mediated by amacrine horizontal cells (Fig. 1 and Fig. 3), through a signreversing synapse to cones called feed-back synapse and a global mechanism, lateral inhibition [1], involved in edge enhancement and image contrasting [6,35], realizing concretely the Mach (boundary brightness overshoot) and Marr (Laplacian zero-crossing edge enhancement) effects. Multiple contrast illusions (Fig. 2) are based on the lateral inhibition principle. Here, we will study how functions of rods and cones are differentially affected during the retinal degeneration of the retinitis pigmentosa and how this pathologic process can be modelled.

II. PHYSIOLOGICAL AND PATHOLOGICAL RETINA

In physiological retina, lateral inhibition causes illusions like the perception of artefactual stripes or spots. The lateral inhibition causes a reinforcement (or a decline) of brightness in a pixel if its neighbours are black (resp. white). This illusion is easy to simulate by computer and is illustrated in Fig. 2. The pathologies of the retina provoke a progressive death of the rods (like in retinitis pigmentosa in Fig. 3), which leads to cones apoptosis, due to non-secretion (by rods) of growth factor favouring the cones survival, causing the disappearance of the lateral inhibition, hence of the contrasting ability. Fig. 3 (Top-right and Bottom-left) shows a confocal slice of a sick retina where we can observe an important loss of both rods and cones in the left superior quadrant (LS). An analysis of the cells interdistance in the three other quadrants shows that the mean interdistance between cones in the peripheral retina (about 15ߤ) is better conserved than the interdistance between rods (about 3ߤ), proving the primary rod degeneracy.

III. NEAREST NEIGHBOUR MODELS

A formal deterministic neural network ܴ of size ݊ is defined by its state variables ሼݔ ሺݐሻሽ ୀଵ… , where ݔ ሺݐሻ denotes the state of the neuron ݅ at time ݐ (equal to '1' if the neuron fires at this time and '0' otherwise). Then the discrete iterative system ruling the change of states in the network is given by the following equations:

ݔ ሺݐ 1ሻ ൌ 1; if ܪ ሺݐሻ ൌ ݑ ݓ ݔ ሺݐሻ ,אሺሻ
ݒ ݔ ሺݐሻݔ ሺݐሻ ߠ, ൌ 0; otherwise where V(i) is a neighbourhood of ݅; ܪ ሺݐሻ plays the role of the somatic electric potential; ݑ denotes an external field; ݓ designates the synaptic weight resuming the interactions of the neuron ݆ on the neuron ݅; ݒ is a non-linear effect coefficient (due to the presence of a triplet of neighbouring neurons firing together); and ߠ is a firing threshold (Fig. 5). The updating of the neuronal states can be operated [11]:

• either sequentially, after having chosen a certain order for the neurons, • or block-sequentially, by parallel updating of each subnetwork of partition ܴ and then activating these subnetworks sequentially, • or in a massively parallel way. The updating rule can be randomized as follows:

ܲ ቀݔ ሺݐ 1ሻ ൌ 1ቚݔ , ݆ א ܸሺ݅ሻቁ ൌ ݁ ு ሺ௧ሻ ் ⁄ 1 ݁ ு ሺ௧ሻ ் ⁄
We can easily show that this rule is the same as the deterministic rule above, when ܶ ൌ 0. The presence of a cone growth factor secreted by neighbouring rods will be taken into account by putting a '1' in a square of 100ߤ ଶ , containing at least one cone (Fig. 6), and a '0' if there is no cone. The deterministic or stochastic dynamics of such neural networks and their robustness has been extensively studied in [2][3][4][5] and [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22]. The occurrence of the values '1' in a grid made of squares, each of 100ߤ ଶ , will be supposed to be ruled by a nearest neighbour Markov random field [24][25][26][27][28][29][30][31][32] and [36,37], which corresponds exactly to the fixed configuration of the random neural networks defined above. There exists a classical unbiased maximum likelihood inference procedure for estimating all the coefficients ሺݑ , ݑ ଵ , ݑ ଶ ሻ of the neural network (in the hypothesis of a translation invariant isotropic field with temperature '1'), since the corresponding statistical structure is exponential [18].

In the example of Fig. 6, the estimates are significantly different between the quadrants LS and RI, showing the presence of a pathologic process. If we assume the network states determined by this nearest neighbour modelling, we have to check if there is an influence of the boundary neurons (in violet in Fig. 5). It is possible to systematically study this influence by searching phase transition parametric conditions, i.e. values of the coefficients ݑ , ݑ ଵ , ݑ ଶ for which the states on the core depend on the states on the boundary. In the parametric circumstances where: ݑ ݑ2 ଵ ݑ ଶ ൌ 0, (cf. Fig. 8 where ݑ ൌ െ3), we observe a phase transition with the dependency of the core on the boundary, which proves that we have to be very careful in fixing the states on the frontiers of the retina [16,17].

IV. SPATIAL RENEWAL BINARY PROCESSES

There exist different alternatives to the nearest neighbour random neural networks, like the reactiondiffusion process, in which one can identify the diffusion coefficient and also the lateral inhibition parameters of the reaction part. There is also a way to escape the spatial Markovian character of the previous model, by supposing that the occurrence of the state '1' at time ݐ 1 in the neuron ݅, given by ݔ ሺݐ 1ሻ ൌ 1, is depending on the states on the first sphere (for the Manhattan distance ܮ ଵ ) centered at ݅, where we can meet a neuron in state '1'. Fig. 8 -Simulated configuration of the states of a 2-dimensional renewal spatial binary random field [23] A spatial random process verifying such a condition is called a renewal spatial random field [20,23] (analogous to the renewal temporal processes observed, for example, while tossing a coin) and we can estimate its parameters in the same way as for the spatial Markov random fields, that is by considering the associated statistical exponential structure in which estimates are of the maximum likelihood, unbiased and almost surely convergent.

We can see, in Fig. 8, a simulation of such a renewal spatial random field, showing configurations sparser than those observed for a spatial Markovian field, correlating more to the data in Fig. 6. A test of data adequacy can be performed over the spatial Markovian and renewal fields and we can retain the structure which best-fits the data, where the first model favours rod-cone interactions and the second model favours long-range interactions between cones.

V. CONTRAST ENHANCEMENT Another way to detect an abnormality in the rod and cone distribution is to use the observed configuration as an artificial retina in order to contrast an input image. The latter can be made of pixels of different grey levels as in Fig. 9 (top-left). A normal retina treating the image must enhance homogeneous zones, by exploiting the fact that a square of pixels (or a peak) having the same medium (respectively high) grey level is reinforced (respectively undermined) by the presence (respectively absence) of a local activation and its boundaries are contrasted because of the absence of external inhibition (Fig. 10). If the distribution of the weights follows the scheme given in Fig. 2 (bottom-right) and if the neurons are dispatched in the four quadrants as in Fig. 3 (bottom left), then the contrasting is realised only in RI quadrant as in Fig. 9 and Fig. 10.

VI. CONCLUSION

We have shown in this paper that a retinal pathology, as in retinitis pigmentosa, could be studied by identifying from observed data, either a Markovian or a renewal spatial random field, leading to the hypothesis of local rod-cone interactions or of a long-range cone-cone interactions. In the first case, the interactions diminishing due to vanishing of rods depicts existence of an action by rods, like the secretion of a specific cone growth factor, which is in agreement with the genetic studies. The conservation of the contrasting retinal function can be tested for the still healthy zones by using, as input, a reference image having a homogeneous zone inside to enhance (This is being extensively used in image processing applications as in Fig. 10). The absence of contrasting power could constitue a good test of loss of this major functionality related to the integrity of the lateral inhibition rod-cone architecture.

Introduction

''In nova fert animus mutatas dicere formas corpora…'' I want to speak about bodies changed into new forms… (Ovid, Metamorphoses, Book 1 st , 10 AD).

In the vertebrate retina, cones are hyperpolarized when illuminated by light, but also receive a depolarizing input when receptors some distance away are illuminated. This antagonistic center-surround response is mediated by amacrine and horizontal cells (Figure 1), through a sign-reversing synapse to the cones often called feedback synapse, the global mechanism being called lateral inhibition [1][2][3]. This surround response is involved in edge enhancement and image contrasting [4][5][6][7][8][9][10][11][12][13][14][15][16] realizing concretely the Mach (boundary brightness overshoot) and the Marr (Laplacian zero-crossing edge-enhancement) effects, used in many image processing applications [17]. A number of contrast illusions (Figures 2,3,4) have been described [18] based on the lateral inhibition principle. In order to examine how rod and cone functions are differentially affected during retinal degeneration (abolishing the contrast), many studies have been done on the genetic level showing that these two cell types have complementary roles during both development and degenerative processes [19][20][21]. For understanding the retinal physiology as well as this pathology, many models [22][23][24][25][26][27][28][29][30][31][32][33][34] are now available which try to mimic relevant adaptation behaviours of the human visual system, like lightness/colour constancy and contrast enhancement, corresponding to the ability of the visual system to increase the appearance of large-scale light-dark or inter-colour transitions, similar to how sharpening with an ''un-sharp mask'' increases the appearance of small-scale edges.

These models use theoretical developments [35][36][37][38][39][40][41][42][43][44] in dynamical systems, especially the study of their attractors. An attractor represents the ultimate evolution of a dynamical system when time tends to infinity; after perturbations, an attractor recovers its stable dynamical features, like its period and amplitude. That requires a rigorous mathematical framework for defining the continuous flow and its convergence speed to attractors, and after its discrete version, i.e. an iteration process representing the succession of states of the dynamical system. These theoretical advances have permitted the development of fast image processing algorithms used in rapid contrasting methods [45][46][47][48][49][50][51][52][53][START_REF] Brunie | Interpretation of multimodal medical images using connectionist and variational methods[END_REF][START_REF] Mattes | Neural networks for parallel contrast enhancement in medical images[END_REF][START_REF] Bosch | Rate and Temporal Coding with a Neural Oscillator[END_REF][START_REF] Demongeot | Bio-array images processing and genetic networks modelling[END_REF][START_REF] Demongeot | Bio-array images processing and genetic networks modelling[END_REF] implemented in realtime processors [START_REF] Fay | Real-time early vision neurocomputing[END_REF][START_REF] Kobayashi | A layered architecture for regularization vision chips[END_REF][START_REF] Nabet | Implementation of front-end processor neural networks[END_REF][START_REF] Kobayashi | Image processing regularization filters on layered architecture[END_REF][START_REF] Shi | Design of linear cellular neural networks for motion sensitive filtering[END_REF][START_REF] Kobayashi | Fast adaptive contrast enhancement method for the display of gray-tone images[END_REF][START_REF] Kobayashi | A High Fidelity Contrast Improving Model Based on Human Vision Mechanisms[END_REF][START_REF] Moini | Vision chips[END_REF][START_REF] Kwon | Design of regularization filters with linear neural networks (image restoration[END_REF][START_REF] Nakahira | An image processing system using image signal multiprocessors (ISMPs)[END_REF], and the development of contouring methods like snakes, snake-splines, d-snakes, which allow a global definition of the boundaries of objects of interest in an image. These algorithms have emphasized the role played by computer implemented procedures, starting from an initial compact, e.g. a sphere, and ending at the final shape of the object's contours after a certain number of iterations [START_REF] Mumford | Boundary detection by minimizing functional[END_REF][START_REF] Kass | Snakes: active contour models[END_REF][START_REF] Leitner | Dynamic Segmentation: Finding the Edge with Snake-Splines[END_REF][START_REF] Cohen | On active contour models and balloons[END_REF][START_REF] Chou | Estimation of the velocity field of 2-d deformable motion[END_REF][START_REF] Chen | Description of complex objects from multiple range images using inflating balloon model[END_REF][START_REF] Kichenassamy | Gradient flows and geometric active contour models[END_REF][START_REF] Bainville | 3D reconstruction using a new dynamic model: the d-snakes[END_REF][START_REF] Faugeras | Complete Dense Stereovision Using Level Set Methods[END_REF][START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF][START_REF] Osher | Level Set Methods: An Overview and Some Recent Results[END_REF][START_REF] Jiang | Level-Set Curve Particles[END_REF]. The corresponding flow is a compact set valued flow, the simplest deriving from a potential [START_REF] Doyen | Mutational equations for shapes and vision-based control[END_REF][START_REF] Demongeot | Compact set valued flows I: applications in medical imaging[END_REF][START_REF] Lorenz | Compact set valued operator[END_REF][START_REF] Demongeot | Multiscale Modeling and Imaging: The Challenges of Biocomplexity[END_REF][START_REF] Demongeot | Dynamical Systems and Bio-medical Image Processing[END_REF][START_REF] Demongeot | A brief history about the analytic tools in medical imaging: splines, wavelets, singularities and partial differential equations[END_REF]. In general, this methodology allows one to rapidly and automatically obtain 3D contours, which is necessary in medical imaging to perform computer aided medical interventions. If the dynamics are conservative in a neighbourhood of an attractor, the flow becomes Hamiltonian, so we then will define the notion of mixed potential Hamiltonian flow. This flow gives a theoretical support to the Waddington's notion of chreod, particularly relevant in embryonic morphogenesis modeling [START_REF] Forest | What is the biological basis of pattern formation of skin lesions? Viewpoint 4. Skin Patterns belong to three main types, determined at three steps of development[END_REF][START_REF] Forest | Cellular modelling of secondary radial growth in conifer trees: application to Pinus radiata[END_REF][START_REF] Forest | Modelling of auxin transport affected by gravity and differential radial growth[END_REF][START_REF] Forest | A general formalism for tissue morphogenesis based on cellular dynamics and control system interactions[END_REF][START_REF] Michon | BMP-2 and BMP-7 play antagonistic roles in feather induction[END_REF], but also serves in image contouring.

Using the previously introduced theoretical notions, we study an enhancement method for contrasting medical images, using either a discrete neural network approach, or its continuous version, i.e. a reaction-diffusion partial differential system [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF][START_REF] Cottet | Mode `les de re ´action-diffusion pour des re ´seaux de neurones stochastiques et de ´terministes[END_REF][START_REF] Catte | Image selective smoothing and edge detection by nonlinear diffusion[END_REF][START_REF] Cottet | Image processing through reaction combined with nonlinear diffusion[END_REF][START_REF] Cottet | Image segmentation using snake-splines and anisotropic diffusion operators[END_REF][START_REF] Weickert | Coherence-Enhancing Diffusion Filtering[END_REF][START_REF] Tschumperle | Restauration d'images vectorielles par EDP[END_REF][START_REF] Tschumperle | Diffusion PDE's on Vector-Valued Images[END_REF]. Indeed, having the goal of providing for a rapid and efficient action [100-142] in precise surgical robotics as well as in disease diagnosis and satellite control imaging, such pre-treatments are performed for contrasting and then contouring images. The medical community, for example, often uses pre-treated anatomical images coming from imaging devices, like MRI or CT-scanner, whose pre-processing involves two fundamental steps: contrasting and contouring. The natural vision executes these two tasks, the first one being based on the architecture of the retina, which uses lateral inhibition to reinforce the perception of the contours of homogeneous objects in a scene. Because the objects of medical interest are homogeneous with respect to their environment (a tumour or an organ are made of cells coming from the same cellular clone), they are well enhanced by using operators processing as in the natural vision. Therefore, we introduce continuous operators generalizing discrete neuromimetic approaches using lateral inhibition as well as analogs of the Hebbian rule for the evolution of synaptic weights.

Results

The results presented in this Section involve consecutive phases of contrasting and segmenting in order to identify objects of interest in an image. The important features of a scene are the prey, predators, and sexual partners. For the detection of these features, the major characteristics are the ''phaneres'', this word coming from the Greek phaneros: visible. The ''phaneres'' in animals and plants are prominent visible tegumentary formations like feathers, scales, hair, petals, skin spots and stripes of various forms and colours. The role of the contrasting pre-treatment in the retina is to rapidly enhance the characteristics (luminance, colour and texture) on the boundaries of the homogeneous zones in a scene in order to improve their perception and extract the features associated to the vital functions like the nutrition, the survival and the reproduction. This process can trigger very fast actions (like escaping a predator) after a stimulus of about 150 ms [136]. Such fast sensory-motor loops need a very simple and rapid mechanism well encoded in the anatomy and in the physiology of the retina (like the center-surround response of cones and rods [1-3]), early before a semantic recognition and denomination of the prey or of the predator. We will give first some results concerning the natural contrasting process both in a natural and in a simulation context.

Pathologic retina

The lateral inhibition mechanism in the retina is due to the presence of feedback synapses of horizontal cells [1,2], which reverse the sign from activation of the cells surrounding that were illuminated (Figure 1 top left). The retina pathologies provoke a progressive death of rods (as in retinitis pigmentosa) followed by the apoptosis of the cones; then, the non-secretion by rods of a growth factor favouring the cones survival, causes the disappearance of the lateral inhibition, hence of the contrasting ability [4,19,20,21]. As shown in the top right and bottom left of Figure 1 on a confocal slice of a sick retina, we observe an important loss of both rods and cones in the left superior quadrant. An analysis of interdistances among cells in the three other quadrants shows that the mean interdistance between cones in the peripheral retina (about 20 m) is better conserved than the corresponding value between rods (about 3 m), proving the primary rod degeneracy.

Contrast illusions

The perception of artefactual stripes or spots comes from the lateral inhibition effect, which causes a reinforcement (respectively decline) of brightness in a pixel if its neighbours are black (respectively white). This illusion effect is visible on the Figures 2 to 4. In Figure 2 (top-left), the Hermann illusion is provoked by the local organization of inhibition and activation between retinal cells, which is described bottom right. The illusion shows bright squares at the intersection of grey stripes and grey squares at the intersection of white stripes. In Figure 2 (bottom-left), the Mach bands illusion gives an enhancement of the vertical lines separating the different grey zones. In Figure 3 top-left, the tangential vision (which allows to escape the macular vision) gives the illusion of a bright reinforcement at the extremities and middle of the white stripes. On the top right, a progressive change of the vertical bright stripe into bright spots (in false colours) is observed during the feathers morphogenetic process in chicken due to a lateral inhibition effect between morphogens (model and simulation are given in [START_REF] Michon | BMP-2 and BMP-7 play antagonistic roles in feather induction[END_REF]). On the bottom left, we can observe bright and grey activities respectively near the center (vertical black line) and the extremities of the white horizontal diamonds. For explaining these illusions, we can simulate a very simple threshold formal neural network (cf. infra) made of 7 neurons, with a lateral inhibition mechanism defined by the parameter values h = w ii = 2, w ii21 =w ii+1 = 20.5, and a sequential updating from the left to the right hand side. The spots activity appears after 3 iterations as a stable steady configuration, and is the discrete analog that the feature created by simulating the continuous reaction-diffusion operator used for modelling feathers morphogenesis [START_REF] Michon | BMP-2 and BMP-7 play antagonistic roles in feather induction[END_REF]. In Figure 4, the sensation of seeing a 3D pyramid is the generalization of the well known Kanizsa polygon effect. It is due to the artefactual prolongation of the white square extremities as white (respectively black) lines in a black (respectively white) dominant neighbourhood. The illusion effects described above are easy to simulate by computer and can serve as external efficacy criterion when different contrasting methods are benchmarked.

Contrasting and contouring images

The enhancement of the grey level on its maximal gradient lines (identical to the geometric locus formed by all the points where the mean Gaussian curvature on the grey surface vanishes) is due to the retinal processing and causes the sensation of contours. By using an enhancement procedure based on the lateral inhibition effect in an formal neural network receiving as input the grey level of an image, we have obtained a good contrast on the boundaries of homogeneous zones either on simulated or on real images. Figure 5 (respectively 6) shows the result obtained after applying a contrasting algorithm on an artificial image (respectively on the NMR slice of a brain tumour).

The contouring step follows the contrasting one, and we see in Figures 5,6 and 7 contours of homogeneous (in grey level) zones resulting from a snake-spline procedure (i.e. an external snakebased procedure with the constraint to keep a spline closed curve at each step) applied over an artificial isolevel square (Figure 5), a brain tumour (Figure 6) and a forest (Figure 7) made of the same species of elements (pixels, cells and trees respectively).

The two steps of contrasting and contouring are based on classical algorithms of neural networks [24,31,32] and snake spline [START_REF] Mumford | Boundary detection by minimizing functional[END_REF][START_REF] Kass | Snakes: active contour models[END_REF][START_REF] Leitner | Dynamic Segmentation: Finding the Edge with Snake-Splines[END_REF][START_REF] Cohen | On active contour models and balloons[END_REF][START_REF] Chou | Estimation of the velocity field of 2-d deformable motion[END_REF][START_REF] Chen | Description of complex objects from multiple range images using inflating balloon model[END_REF][START_REF] Kichenassamy | Gradient flows and geometric active contour models[END_REF][START_REF] Bainville | 3D reconstruction using a new dynamic model: the d-snakes[END_REF], but they can involve new methods coming from biomimetic procedures. We will describe rapidly four such new methodologies and give examples of their application to real satellite or medical images.

1)

A chemotactic operator. If we denote, at time t and pixel x, g(x,t) as the grey level function, we can consider g as a food or substrate, which living entities (like bacteria) can eat, being attracted from the image boundaries (where they are first located) by a chemical gradient linked to the substrate. Let us denote the bacterial concentration by b(x,t). We can consider the following equations, which constitute a new image processing operator [START_REF] Demongeot | Dynamical Systems and Bio-medical Image Processing[END_REF]142]:

Lb=L t~L b +g j j max Db{x+ b+g ð Þ Lg=Lt~L g Dg{ Kg bze ð Þ= bz1 ð Þ ð1Þ
with Neumann conditions on the image boundary, where +g j j max denotes the maximum value of the g gradient norm, x is the attractive chemotactic constant, L b (respectively L g ) is the diffusion coefficient of the bacterial concentration (respectively grey level), K (respectively Ke) is the maximal (respectively minimal) grey consumption rate of bacteria. These equations imply that the bacteria move towards the concentration of grey considered as a chemo-attractant to consume. They also diffuse as the grey level with respectively the diffusion constants L b and L g . The Figure 7 bottom shows the progressive treatment of the image of a Chilean forest presenting the same characteristics of internal homogeneity as a tumour (the trees replacing the cells); due to the fact that the trees (like cells) belong to the same genetic lineage. After reaching their asymptotic values, the dynamics of contrasting implemented in a discrete scheme of the partial differential equations (PDE) (1), stops and this processing step can be followed by a snake spline contouring step.

2)

A viability contouring operator. If we minimize the following function,

aS K t ð Þ ð ÞzbV K t ð Þ ð Þzc ð LK t ð Þ 1= +g x ð Þ k k ½ dx ð2Þ
we obtain a new snake operator [START_REF] Kichenassamy | Gradient flows and geometric active contour models[END_REF][START_REF] Demongeot | Dynamical Systems and Bio-medical Image Processing[END_REF], where K(t) is a compact object of interest moving toward a limit set K('), whose external surface S as well as its inner volume V are minimized, allowing a contouring with real gloves (precise contour) contrarily to mittens (convex envelop) often observed with the Mumford-Kass-Terzopoulos algorithm in Figure 7 [START_REF] Mumford | Boundary detection by minimizing functional[END_REF][START_REF] Kass | Snakes: active contour models[END_REF]. We see in Figure 7 (top and bottom right) the contouring done by imposing a bicubic spline to the boundary at each time step [START_REF] Leitner | Dynamic Segmentation: Finding the Edge with Snake-Splines[END_REF][START_REF] Cohen | On active contour models and balloons[END_REF], followed by a 3D spline smoothing. Many other approaches can also be used for controlling the active-shape models. This is the case in the level set methods used for computing and analyzing the motion of an interface in two or three dimensions by modelling the velocity vector field through Euler-Lagrange or Hamilton-Jacobi PDE's [START_REF] Faugeras | Complete Dense Stereovision Using Level Set Methods[END_REF][START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF][START_REF] Osher | Level Set Methods: An Overview and Some Recent Results[END_REF][START_REF] Jiang | Level-Set Curve Particles[END_REF]. These PDE's can be used to model the segmentation of a moving 3D object (like the heart) giving a particular status to the pixels having a maximal velocity or acceleration of their grey levels. This procedure has been used for segmenting the pericardium [131].

3)

A non-isotropic reaction-diffusion operator. If we consider the grey level function g(x,0) as the initial image, we can follow the transient behaviour of the non-linear diffusion operator defined in [START_REF] Cottet | Mode `les de re ´action-diffusion pour des re ´seaux de neurones stochastiques et de ´terministes[END_REF]:

Lg=L t~Ldiv 1 0,s ½ +(G Ã g j j :+g À Á À Á ð3Þ
Here G is a Gaussian kernel of fixed variance and with Neumann conditions. Its asymptotes correspond to a constant grey level suppressing the objects of interest inside the image. For that reason, we consider now a non-isotropic reaction-diffusion operator defined in [START_REF] Cottet | Mode `les de re ´action-diffusion pour des re ´seaux de neurones stochastiques et de ´terministes[END_REF][START_REF] Cottet | Image processing through reaction combined with nonlinear diffusion[END_REF][START_REF] Cottet | Image segmentation using snake-splines and anisotropic diffusion operators[END_REF]:

Lg=Lt{div L+g ð Þ~0,dL=dtzL=t~s 2 P +g =t, if +gj j ws ~+gj j 2 P +g z3 s 2 { +gj j 2 Id=2, if +gj j ƒs ð4Þ
where L is a 262 matrix and P =g is the orthogonal projection matrix:

P +g ~1= +gj j 2 Lg=Lx ð Þ 2 -Lg=Lx ð ÞLg=Ly ð Þ -Lg=Ly ð Þ Lg=Lx ð Þ Lg=Lx ð Þ 2 ! ð5Þ
In the equations above, the diffusion constant L becomes variable with the time t and its evolution equation is similar to the Hebbian rule of a discrete neural network operator. Treated images are obtained at the asymptotic state of the PDE dynamics as for neural networks [48,49] with lateral inhibition (Figure 6). A comparison done in [START_REF] Cottet | Image segmentation using snake-splines and anisotropic diffusion operators[END_REF] shows that the asymptotes of this non-isotropic operator are better than for some of the operators described earlier. More generally, we can notice in the other PDE approaches:

a) The application of the pure heat operator [145] quickly leads to a constant grey level b) In the Perona-Malik operator [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF], the viscosity is different within a region and across its boundary in order to encourage smoothing inside the region of interest; this operator can be used transiently for this purpose before the non-isotropic reaction-diffusion operator c) The Catte ´-Lions-Morel-Coll algorithm [START_REF] Catte | Image selective smoothing and edge detection by nonlinear diffusion[END_REF] gives a good contrasting during the transient behaviour of the operator, but has the same asymptotes as for the pure heat algorithm (even it is reached more slowly) d) The non-isotropic reaction-diffusion operator [START_REF] Cottet | Mode `les de re ´action-diffusion pour des re ´seaux de neurones stochastiques et de ´terministes[END_REF][START_REF] Cottet | Image processing through reaction combined with nonlinear diffusion[END_REF][START_REF] Cottet | Image segmentation using snake-splines and anisotropic diffusion operators[END_REF] offers a reasonable asymptotic processing e) The Weickert operator [START_REF] Weickert | Coherence-Enhancing Diffusion Filtering[END_REF] permits the completion of interrupted lines or the enhancement of flow-like structures by choosing the appropriate smoothing direction in anisotropic processes in spirit to the Cottet-Germain filter [START_REF] Cottet | Image processing through reaction combined with nonlinear diffusion[END_REF] f) The Tschumperle ´-Deriche operator [START_REF] Tschumperle | Restauration d'images vectorielles par EDP[END_REF][START_REF] Tschumperle | Diffusion PDE's on Vector-Valued Images[END_REF] allows the regularization of velocity vectors fields in 4D imaging (acquired for example during the motion of a 3D camera).

4)

An attentional focus operator. For focusing on only one region of interest, we have to change the image input on an artificial neural network [START_REF] Bosch | Rate and Temporal Coding with a Neural Oscillator[END_REF]. This input can be constant [24,31,32], stochastic [47][48][49][50][51][52][53][START_REF] Brunie | Interpretation of multimodal medical images using connectionist and variational methods[END_REF] or deterministic periodic [START_REF] Bosch | Rate and Temporal Coding with a Neural Oscillator[END_REF]. This last coding mimics the information storage inside the hippocampus in which the functional unit, made of two neurons in mixed inhibition/activation interaction (Figure 8 top left) has an attractor limit cycle. We can locally synchronize, using an evocation stimulus, and desynchronize, by introducing noise on the inter-unit interactions, the periodic activities corresponding to initially non phaselocked neurons. In this way, we enhance considerably (by forcing the units to add their maximal activities at the same time) the grey level on the zones of local synchronization (Figure 8 E bottom right). Then, by thresholding and segmenting, we get the parts of the initial image (Figure 8 A top right) on which the attentional focus has been exerted (Figure 8D, E, F top right).

Computer assisted interventions

For introducing and driving medical or surgical tools (like needles, electrodes, bistouries) into the human body [118-135], one needs to segment and contour (after contrasting) zones of interest to avoid (as indicated by red zones in Figure 9 left representing tissues of lungs on the top and cardiac muscle on the middle and bottom) or to reach (blue arrow in Figure 9 right indicating a pericardial effusion). This example gives a good illustration of what can be exploited from the contrasting and contouring operators in order to go farther than the descriptive level for diagnosis. That is to really improve some medical procedures, one must automate the process completely, thus replacing the human actor without any loss of speed or precision [119][120][121][122][123][124][125][126][127][128][129][130].

Discussion

Interest of the biomimetic approach

The biomimetic approach used in numerous methods presented in this paper, especially for the contrasting phase, exploits the efficiency of visual data processing procedures that have been selected by natural evolution. These procedures represent an optimum in terms of economy of implementation (small number of living elements involved, like cells, tissues, vessels, etc), speed and precision. They also are based on operations that come after processing by the retina and visual areas, thus providing high level semantic neural networks that define the symptomatology related to the observed medical reality. The extraction of semiotic characteristics of objects of medical interest that have been enhanced and contoured using biomimetic methods allows medical signs and symptoms to be organized in syndromes, thus facilitating the diagnosis process. The concept of biological information encoded in a genetic program that controls development forms a major part of the semiotic metaphor in biology. The development plan is seen as being analogous to a computing program, and ''semiotics of nature'' studies the structural relations as explored by molecular and evolutionary biology [137]. Y.L. Kergosien [138] advocates a semiotics of nature in an epistemological sense for analysing interacting biological systems, in order to increase the precision of terms such as ''signal'' in biology or ''symptom'' in medicine, and to develop new themes of inquiry into the nature of their biological or medical signification. The Kergosien approach indeed allows for a concept of natural signification. The adaptation of an animal to a specific function is seen as the realization of the natural metaphor [137]. This is the case for retinotopically arranged neuronal sets that code for homogeneity features (brightness, colour, texture, etc), oriented contours, and corners of an object. Simultaneous representation by colour neurons, complex model neurons (with oriented receptive fields), and hypercomplex model neurons (responding to corners) makes attention and recognition robust and reliable, in the framework of emergent abilities of optimized complex systems [139][140][141].

The bio-inspired image processing methods also have a tendency to use an information encoding that provides for optimum information storage and query, as done in mnemonic structures like hippocampus. In general these structures possess their own formats of information encoded in periodic temporal neuronal activities that we can mimic to optimize both compression and retrieving procedures [40]. All these neural treatments can induce illusions and artefacts. But the knowledge about their origin can be used for preventing such abnormalities in the low level (contrasting and contouring) as well as in the high level (semantic assignation and recognition) image processing steps. The neural treatments need also to avoid pathologic processing, due to a non-optimal number of their neurons and/ or interactions and to a non robust value of their parameters. To that precise purpose, a deep scientific knowledge about the physiology and the pathology of the retina constitutes an unavoidable inheritance.

Limits of the biomimetic approach

In order to be faster, the methods mimicking the natural process of vision need to be made parallel as in the real neuronal systems. But the attractors of the dynamical systems permitting contrasting and contouring of the images are highly dependent on their modality of implementation, particularly on their updating mode. In general, the fixed configurations obtained by simulating such systems are robust with respect to the mode of updating, but it is not the case for the periodic neural activity we have used in attentional focusing (Figure 8). Hence it is convenient to be very careful until the final step of algorithmic implementation.

The imitation of nature does not push to avoid theoretical studies on the spatio-temporal processes used in artificial vision [142]. Only this fundamental approach is able to finally guide the methodological choice with arguments as fast calculation speed [143], precision, accuracy, and minimal algorithmic complexity. Indeed, these good properties constitute the main criteria for selecting robust, fast and precise image processing tools for reliable procedures of computer aided surgical and medical intervention [118][119][120][121][122][123][124][125][126][127][128][129][130][131][132][133][134][135].

Materials and Methods

Discrete operators

Contrast enhancement. A large number of methods of contrast enhancement have been used in the past to reinforce the grey level gradient on the boundaries of objects of interest. These methodologies can be classified following a typology, based on the mathematical tools underlying algorithms:

-classical filtering (e.g. Gaussian [28]), PDE filtering [START_REF] Faugeras | Complete Dense Stereovision Using Level Set Methods[END_REF][START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF][START_REF] Osher | Level Set Methods: An Overview and Some Recent Results[END_REF]] the simplest being the heat operator analog for the Gaussian filter, grey histogram thresholding [26,30], entropy techniques [25], adaptive filtering [START_REF] Kobayashi | Fast adaptive contrast enhancement method for the display of gray-tone images[END_REF] -multi-scale, in particular wavelets [START_REF] Demongeot | Multiscale Modeling and Imaging: The Challenges of Biocomplexity[END_REF]) approach [27,31] -fuzzy clustering [29] -specific hardware implementation [START_REF] Kobayashi | A layered architecture for regularization vision chips[END_REF][START_REF] Kobayashi | Image processing regularization filters on layered architecture[END_REF] for real-time procedures [22], using either the simplest PC based [23] or the most sophisticated architectures (SIMD [32,52] or MIMD [START_REF] Nakahira | An image processing system using image signal multiprocessors (ISMPs)[END_REF])

-neural networks techniques both discrete [48][49][50][51][52][53][START_REF] Brunie | Interpretation of multimodal medical images using connectionist and variational methods[END_REF][START_REF] Mattes | Neural networks for parallel contrast enhancement in medical images[END_REF] and continuous [START_REF] Cottet | Mode `les de re ´action-diffusion pour des re ´seaux de neurones stochastiques et de ´terministes[END_REF][START_REF] Cottet | Image processing through reaction combined with nonlinear diffusion[END_REF][START_REF] Cottet | Image segmentation using snake-splines and anisotropic diffusion operators[END_REF].

We will focus in the next Section on the neural networks techniques which are the closest to the natural vision processing.

Definition of a formal neural network. A formal deterministic neural network R of size n is defined by its state variables {x i (t)} i = 1, … ,n , where x i (t) denotes the state of the neuron i at time t (equal to 1 if the neuron fires at this time and to 0 if not). Then the discrete iterative system ruling the change of states in the network is given by the following equations:

x i tz1 ð Þ~1, if H i t ð Þ~X w ij x j t ð Þwh, j [ V i ð Þ ~0, if not,
where V(i) is a neighbourhood of i, H i (t) plays the role of the somatic electric potential, w ij designates the synaptic weight representing the influence of the neurone j on the neurone i and h is a firing threshold. The updating of the neuronal states can be operated:

-either sequentially, after having chosen a certain order for the neurones, -or block-sequentially, by operating the updating in parallel in each sub-network of a partition of R and by afterwards activating these sub-networks sequentially, -or in a massively parallel fashion if only one sub-network exists.

Input in a neural network

If an input I i (t) is sent to neuron i at time t, it is merged with the information coming from the neighbourhood V(i) in order to build the somatic potential H i (t):

H i t ð Þ~X w ij x j t ð ÞzI i t ð Þ j[V i ð Þ
A very simple way of generating such inputs is to choose, for each time interval E k (supposed to be independent of the others) between the two consecutive inputs 1, the k th and the (k+1) th , the truncated geometric distribution: Prob({E k #T i }) = 0 and Prob({E k = m . T i }) = p i (12p i ) m2Ti21 , where T i and p i denote respectively the refractory period and the spike occurrence frequency on the afferent fiber i bringing the electric input to the neuron i. The truncated geometric processes are independent or correlated between fibers. In Figure 5, we can see the activity of a formal neural network activated by a non-homogeneous input representing the initial image (top-left), and after iterating the neuronal firing, we obtain as mean asymptotic behaviour (bottomleft). The coding is obtained by taking T i and p i proportional to the grey level of the initial image. Image on the top-right is representing the dynamics of the synaptic weights {w i0j (t)} jMV(i) which follows a Hebbian rule reinforcing the weight w i0j (t) if i 0 and j had the same firing activity at time t:

w i0j t ð Þ~F w i0j t{1 ð Þza x i t ð Þx j t ð Þ{ b 1{x i t ð Þ ð Þx j t ð Þz 1{x j t ð Þ À Á x i t ð Þ Â Ã !
where F is a sigmoidal function of arc-tangent type. The initial distribution {w i0j (0)} jMV(i) is chosen dog-like (i.e. a difference of Gaussian distribution centred at i 0 , the negative Gaussian having the greatest variance as shown in the red dog G in Figure 5 (topright)), for mimicking the lateral inhibition. The image treated is shown in grey level in Figure 5 (bottom-right), from initial to treated asymptotic image. We see that the square having a medial activity is enhanced by the lateral inhibition expressed by the dog function and its final level after iterating the network until it reaches its asymptotic firing regime, has a level clearly augmented (see the orange square on the bottom left and the enhanced ''mesa'' on the bottom right). Such a simulation highly suggests that an analogy between pixels and neurons can be made allowing the transfer of neural filtering techniques in image processing [24,31,48].

Gradient enhancement by a neural network

Image enhancement procedure. We now present in 4 steps, the essentials of a method, easy to parallelize based on the same principles as proposed in [31]:

1)
reduction of a 5126512 NMR image in a 2566256 image by averaging each block of 4 neighbour pixels, in order to obtain the input image (cf. Figure 6 top left).

2)

use of this image as the mean configuration of an input geometric random field transformed by a 2566256 uni-layer neural network implemented in parallel; this network has an internal evolution rule, realizing a treatment of the input signal very close to a cardinal sine convolution, mimicking the lateral inhibition and favouring the occurrence of a very steep gradient on the boundary of homogeneous (in grey level) objects of interest in the processed image. In Figure 6, the object of interest is a brain tumour, its homogeneity coming from the same clonal origin of all its tumour cells. 3) use of the gradient, built by the neural network as the potential part of a mixed potential Hamiltonian differential system, whose Hamiltonian part is given by the initial grey level (before the action of the neural network). 4)

obtaining boundaries of homogeneous objects as limit cycles of the differential system by simulating trajectories of the system in the different attraction basins.

The step 2 consists of defining the input from a geometric random field, i.e. a collection of geometric random processes such that, if p i (t) denotes the probability to generate a spike on the afferent fiber i to the neuron i at time t, we have: p i (t) = 0, if t2s i #T i , where s i is the time of the last 1 on the fiber i before time t and T i denotes the refractory period, chosen as a constant equal to R. p i (t) = a i sin + (v i (t2s i 2R)), if t2s i .R, where sin + denotes the positive part of the sine.

In order to incorporate an adaptation learning effect, a Hebbian evolution of the w ij 's is chosen based on the reinforcement of equal grey activities in the same neighbourhood:

w ij tz1 ð Þ~log X p i s ð Þp j s ð Þ . t , sƒt
where w ij (0) values come from a dog (difference of Gaussian) distribution of j centred at i, for each i, for mimicking the lateral inhibition. This formula corresponds to the fact that w ij (t) is just the non-centred covariance function between the p i (s)'s and the p j (s)'s; if v i2 v j and R are small, w ij (t), when t tends to infinity, tends to log((a i a j /2)sin(v i2 v j )/(v i2 v j )).

Image coding

After normalization of the grey level g(i) in the pixel i between 0 and 1, we take:

a i ~g i ð Þ and v i ~l g i ð Þ
and we start the procedure by iterating the deterministic neural network. It is easy to prove that the probability p i to have 1 as output of the neuron i at time t, just before renormalization, is about proportional to:

p' i ~X X pi s ð Þpj s ð Þ=t h i j[V i ð Þ sƒt
This last formula has been used to make the gradient enhancement visible in Figure 6 (top-middle). The behaviour of the function p' i is similar to a convolution by a cardinal sine function, because of the approximate asymptotic formula:

p 0 i ~Xj[V i ð Þ log a i a j 2 À Á sin v i -v j À Á v i -v j À Á À Á
It is easy to verify that this convolution reinforces the ''plateau'' or ''mesa'' activities in grey level (or white if necessary). Such activities correspond, in medicine, to pathological objects to be considered as targets during the treatment (like tumour in which the same clone of cells gives a homogeneous response in absorbance or resonance) or to physiological objects (like a tissue made of cells having the same function) to be avoided during the treatment. Figure 6 shows the result of a gradient enhancement by the network for a brain tumour. Let us finally remark that we get objects treated at the asymptotes of the network dynamics. We do not need a stop criterion after few steps of processing and the method is easy to parallelize [START_REF] Mattes | Neural networks for parallel contrast enhancement in medical images[END_REF][START_REF] Nabet | Implementation of front-end processor neural networks[END_REF].

Continuous operators

The final aim of these methods is to offer a set of continuous operators adapted to segmentation of grey singularities or grey peaks (0-dimensional objects like micro-calcifications), grey anticlines (1-dimensional objects like vessels) or grey ''mesas'' (2dimensional objects like tumours or functional regions). The problem of segmentation of more complicated objects (fractal objects like diffused tumours affecting, for example, the conjunctive tissue) is open and demands that other variables like texture based one's (e.g. the local fractal dimension or the wavelets coefficients) need be taken into account instead of or along with the grey level.

Let us consider now a compact state set E included in R 2 and a temporal set T included in R + or N, depending on the continuous or discrete version of time used. Let K(E) denotes the set of all compacts of E. If we provide K(E) with the Hausdorff topology (defined by the Hausdorff distance d between subsets), we can define a compact set valued (csv) flow Ø as a continuous application of K(E).T to K(E), which is a semi-group:

V K,t ð Þ [ K E ð Þ:T,V s [ T, K,t ð Þ,s ð Þ ~w K,tzs ð Þ
Because K(E) is a metric space, which is compact if E is compact, we can apply the operators limit and basin as defined in [36,37] to the set valued flow Ø, and hence define the notions of attractor and of stability basin. We will give some examples of csv flows, whose attractors are objects to be contoured in image processing, or final shapes to be obtained at the end of any morphological development, these targets being often the same.

Potential flows. In snake contouring [START_REF] Mumford | Boundary detection by minimizing functional[END_REF][START_REF] Kass | Snakes: active contour models[END_REF][START_REF] Leitner | Dynamic Segmentation: Finding the Edge with Snake-Splines[END_REF][START_REF] Cohen | On active contour models and balloons[END_REF], the aim is to obtain the boundaries of an object of interest by progressively deforming the boundaries of an initial well-known set K(0) (e.g., a sphere) placed outside (respectively inside) the object, and whose deformation K(t) causes the decrease (respectively increase) of a potential function P [START_REF] Kichenassamy | Gradient flows and geometric active contour models[END_REF] such as:

P K t ð Þ ð Þ~aS K t ð Þ ð ÞzbV K t ð Þ ð Þzc ð LK t ð Þ 1= +g x ð Þ k k ½ dx
in which S(K(t)), V(K(t)), hK(t), and g(x) denote respectively the external area, the inner volume, the boundary, and the grey level at the point x of the compact K for iteration t. The gradient iterations of P correspond to a discrete potential flow. For obtaining the continuous version, it suffices to use a potential ''mutational'' equation [START_REF] Doyen | Mutational equations for shapes and vision-based control[END_REF][START_REF] Demongeot | Compact set valued flows I: applications in medical imaging[END_REF][START_REF] Lorenz | Compact set valued operator[END_REF]. We can also add splines-like terms, e.g. d# hK(t) C(x)dx, where C(x) = (h 2 g/hx 1 2 )(h 2 g/hx 2 2 is the mean Gaussian curvature at x (in order to minimize the total variation of the local curvature like for the splines functions), plus a mean square criterion forcing hK(t) to pass in the vicinity of points known a priori with fixed curvatures (in particular singular parabolic or saddle points, if their localization is known a priori).

Mixed potential Hamiltonian Segmentation. The continuous modelling allows stable evolution of differential operators such as gradient or Laplacian. Our segmentation consists in building a differential equation system whose stable manifold is the surface of the object we are looking for. Finding this manifold turns out to be a particular case of the surface intersection problem and provides an immediate analytical representation of the surface. The other major advantages of this method are to perform segmentation and surface tracking simultaneously, to describe complex structures in which branching problems can occur if the segmentation is purely local, and to provide accurate and reliable results.

Let us first consider the 2D problem. The central idea of the method is based on the Thom-Sebastiani conjecture [35] concerning the differential system:

x 0 t ð Þ~F x,y ð Þ, y 0 t ð Þ~G x,y ð Þ
In the neighbourhood of a stable singularity or of a limit cycle of the corresponding velocity vector field supposed to be continuous, let us suppose that we can decompose the system into two parts, a potential and a Hamiltonian one, such as:

x 0 t ð Þ,y 0 t ð Þ ð Þ ~{grad P x,y ð Þzham H x,y ð ÞzR x,y ð Þ,
where the residue R(x,y) tends to 0 when (x,y) tends to the stable singularity or to the limit cycle. Such decomposition has been proven for a large class of Lie ´nard systems [41][42][43][44]. The Thom-Sebastiani conjecture assumes that this result still holds by considering sufficiently regular systems. We will exploit systematically in the following, this possibility to consider a contour as the limit-cycle of a mixed potential Hamiltonian system. In fact, we consider now the boundary surrounding a 2D object with an approximately homogeneous grey level g, thus verifying:

g x,y ð Þ~k,
where k is a constant

The corresponding curve is represented with parametric coordinates by:

x~x t ð Þ; y~y t ð Þ
The continuous modelling implies the existence of the first derivatives of g; so a solution should verify the following equation obtained by differentiation of g(x,y) = k:

x 0 t ð ÞLg=Lxzy 0 t ð ÞLg=Ly~0
A particular solution of this equation is: x'(t) = hg/hy, y'(t) = -hg/hx, but this system does not provide a stable solution; a perturbation (due to noise) moving the curve away from the initial contour line could not be corrected. That is why we add a component which brings the curve back to the contour line defined by g(x,y) = k, according to the steepest slope line of the function (g-k) 2 . We thus obtain:

x 0 t ð Þ~Lg=Ly{bLg=Lx= +g x ð Þ k k, y 0 t ð Þ~{Lg=Lx{bLg=Ly= +g x ð Þ k k
This system consists in two parts: the first one corresponds to an ''edge tracking'' component and the second one is a kind of ''elastic force'' which allows noisy image processing. The b parameter allows to balance these two terms. The system may be solved by numerical analysis methods with initial conditions, like the Runge-Kutta-Gear method. The parametric representation of the curve is then directly obtained. This continuous method can be applied in 3-dimensions to look for particular features of the surface of an object of interest. Let us consider such a surface defined by: f(x,y,g) = constant, parameterized by:

x~x t,h ð Þ, y~y t,h ð Þ, g~h
Our boundary tracking method can be implemented as follows: the algorithm starts with a point on the surface with a grey value h. For each slice of level h, the differential system is solved in order to obtain a closed curve. From some points of this curve, we follow the object surface until the next (k+1) slice by building new 2D differential systems in slice level planes. The algorithm stops when all slices have been processed or when the object surface has been entirely described. This method allows to find automatically all the components of a complex object in which branching problems may occur and to determine how they are linked together. This possibility is one of the major advantages of the method because surface reconstruction from a set of contours is a critical step for complex structures. Classically, interpolation between contours is performed by triangulation techniques or by creating intermediate contours with dynamic elastic interpolation. But these methods need sometimes interaction with the user. In our method the surface modelling is performed in the segmentation step. This algorithm has been tested on MRI images for stereotaxy before stimulation needle introduction or brain tumour puncture [118][119][120][121][122][123][124][125][126][127][128][129].

The remarkable Gaussian line. Homogeneity is not always a stable characteristic of an anatomical structure. So we present now a differential system performing H(g) = 0, where H is an operator similar to the Laplacian or Marr-Hildreth detectors. Let us define the remarkable Gaussian line of a peak as the set of points where the mean Gaussian curvature of the peak vanishes (Figure 10). Its equation writes [41]:

H x,y ð Þ~L 2 g Lx 2 À Á L 2 g Ly 2 À Á { L 2 g LxLy À Á 2 ~0
If H' = |H|, let consider the mixed potential Hamiltonian system [42][43][44] obtained as follows:

dx=dt~{aLH' . Lx H x,y ð Þ . grad g ð Þ k k 2 h i zbLH=Ly dy=dt~{aLH' . Ly H x,y ð Þ . grad g ð Þ k k 2 h i {bLH=Lx
We consider in Figure 10 bottom right the new grey function H(x,y) instead of the function g(x,y) at each pixel (x,y) and we display bottom left the mixed potential Hamiltonian differential system above of which the characteristic line is a limit cycle, called the Hamiltonian contour. Its first term is of steepest descent dissipative nature and along the flow, the trajectories converge to the zeros of H'(x,y). On the set of the zeros of H'(x,y), the second Hamiltonian term of the differential system which is of conservative type, becomes preponderant. Parameters a and b can be used to tune the speed of convergence of the differential system to the limit cycle. The usual Runge-Kutta-Gear discretization scheme yields ultimately for the differential system an algorithm which is quite easy to implement. On each pixel (boundary effects are neglected), the function H(i,j) reads:

H i,j ð Þ~g iz2,j ð Þ{2g iz1,j ð Þzg i,j ð Þ ½ g i,jz2 ð Þ{2g i,jz1 ð Þzg i,j ð Þ ½ { g iz1,jz1 ð Þ {g i,jz1 ð Þ{g iz1,j ð Þzg i,j ð Þ ½ 2
An important property of the remarkable Gaussian line is that in the case of a Gaussian peak, it contours the projection of a volume equal to 2/3 of the total volume of the peak. This property remains available with a good approximation in case of moderate kurtosis and skewness of the peak. An advantage of this technique is that we do not perform a direct segmentation of the grey level. Thus the segmentation is much finer than the corresponding one performed by the watershed lines method or by its variant with markers [103]. We only segment the upper part of the peak and then we multiply by 3/2 the activity integrated inside the remarkable line. This approach is interesting because the lower part of the peak is often noisy. The method seems particularly efficient when the peaks are well separated. If they are close (see Figure 10 bottom right), then we need to tune the parameters a and b and to start the trajectories inside the peaks. For finding a contour line inside, we can: 1) calculate the total variation V(h) = # C(h) ||=g(x)||dx of the gradient norm||2=g|| along a contour line C(h) of level h 2) both decrease and increase h towards two limits h 1 ,h and h 2 . h in order to find an intermediary value of V(h) greater than the two values V(h 1 ) and V(h 2 ) calculated at the extremities h 1 and h 2 . Then C(h 1 ) and C(h 2 ) constitute an annulus whose intersection with the remarkable line is not empty 3) choose the initial condition on C(h 2 ) for starting the simulation of the differential system.

Eventually, we can notice that the remarkable Gaussian lines can serve for matching images or objects of interest, for example, in the case of comparing images to a reference coming from an atlas. They constitute a feature in general more robust than parabolic or saddle singularities sensitive to perturbations causing local skewness of the grey peaks. Watershed contouring. The watershed line is a concept firstly defined by geographers in order to characterize the main features of a landscape: a drop of rain that reaches the ground will flow down to a sea or an ocean. In the case of France, the watershed line splits the country in two parts, the Atlantic zone and the Mediterranean zone. Those zones are called 'catchment basins', and the oceans are the minima of them, i.e. the attraction basins of the gradient operator which corresponds to the gravitational dynamics of the drop on the steepest gradient lines of the relief surface. They define a partition of this relief, and the boundaries of catchment basins define on the pixels plane the watershed lines [105-109]. These lines are confounded in regular cases with the crest lines surrounding the catchment basin. It is easy to understand the interest of this concept in image processing: grey level images can be considered as relief structures, and the watershed lines are a good way to separate light (low grey level) zones from dark (high grey level) ones. It is particularly interesting to determine the watershed lines of the symmetrical reverse landscape obtained by considering the new grey level 1-g, where g is the initial normalized grey level obtained after the contrasting step and after fixing the maximum of g as a normalized value equal to 1. The watershed lines verify variational principles: i) when progressively fulfilling with water a catchment basin, its inner area passes through a series of inflexion points corresponding to the successive saddle points reached by the water. Each inflexion point corresponds to a local maximum of the second derivative of the inner area; ii) for a given inner area, the watershed lines are those containing the maximum of water. The watershed line is computed on a discrete image, by immersion simulation, locating it on the meeting points of several catchment basins (Figure 10). First discrete algorithms of watershed lines computed by immersion simulation were proposed in [105-109] with a discrete operator. In [103,110], the watershed line is computed on the reverse image, in order to have one and only one local maximum of the original image into each catchment basin of the reverse image. The resulting labelling (still not a partition) is done on the original image. We used the Vincent-Soille algorithm [105] on discrete images with a linear complexity (about 7,25 n, where n denotes the number of pixels in the image). It can be used also in 3 dimensions.

Reaction-diffusion contrasting. Several methods of image contrasting by using differential linear or non-linear operators have been proposed [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF][START_REF] Cottet | Mode `les de re ´action-diffusion pour des re ´seaux de neurones stochastiques et de ´terministes[END_REF][START_REF] Catte | Image selective smoothing and edge detection by nonlinear diffusion[END_REF][START_REF] Cottet | Image processing through reaction combined with nonlinear diffusion[END_REF][START_REF] Cottet | Image segmentation using snake-splines and anisotropic diffusion operators[END_REF][START_REF] Weickert | Coherence-Enhancing Diffusion Filtering[END_REF][START_REF] Tschumperle | Restauration d'images vectorielles par EDP[END_REF][START_REF] Tschumperle | Diffusion PDE's on Vector-Valued Images[END_REF]. These methods can be parallelized as for the neural networks and we will show in the following that there exists a deep relationship between the discrete neural network approach and the continuous differential operator approach.

1) The Catte ´-Lions-Morel-Coll non-linear diffusion operator. It is well known that the solution of the heat differential operator:

Lu=Lt~k:Du~k:div grad u ð Þ ð Þ
is the Gaussian kernel, with variance equal to s 2 = 2 kt, by choosing as initial conditions u(.,0) the grey level. This property has suggested [START_REF] Catte | Image selective smoothing and edge detection by nonlinear diffusion[END_REF] the use of another differential non-linear diffusion operator:

Lu=Lt~div g grad G Ã u ð Þk k ð Þ :grad u ð Þ ð Þ ,
where G is a Gaussian kernel and g is a non-negative nonincreasing function on R + verifying g(0) = 1 and g tends to 0 at infinity; in practice, we can choose for g a set function, whose value is 1 on the interval [0,S] and 0 on]S, +'[: there is diffusion if and only if ||grad(G * u)||#S and, after a certain transient, it remains a gradient only on the boundary of sufficiently discriminable objects. For example, Figure 11 presents images after some hundreds of iterations, showing the gradient on the boundary of brain structures. The end of the procedure as for the heat operator (Figure 11 left (a))

shows that diffusion wins, giving a constant grey level at the asymptotic state. In order to improve the method of getting the contrasted image at the asymptotic state of the simulation, we must add a reaction term in order to obtain the final expected image as the attractor of a differential reactiondiffusion operator, like for the iterative discrete neural network as in Figure 11 left (c).

2) The non-isotropic reaction-diffusion operator. By searching a continuous operator having as discrete finite elements scheme a deterministic neural network system similar to that presented in Section 2, it has been proposed [START_REF] Cottet | Mode `les de re ´action-diffusion pour des re ´seaux de neurones stochastiques et de ´terministes[END_REF][START_REF] Cottet | Image processing through reaction combined with nonlinear diffusion[END_REF][START_REF] Cottet | Image segmentation using snake-splines and anisotropic diffusion operators[END_REF] with direct reference to the discrete neural network approach [48,49,52] a new reaction-diffusion operator. Let us recall the deterministic neural network with threshold 0 defined by:

x i tz1 ð Þ~1, if H i t ð Þ~X w ij x j t ð Þ w 0, j[V i ð Þ ~0, if H i t ð Þ~X w ij x j t ð Þ ƒ 0, j[V i ð Þ
where V(i) is a neighbourhood of i. If we suppose the neural network to be 2D and infinite, lets us denote by (i 1 , i 2 ) the position of the neuron i, where i M Z 2 ; if w ij are symmetrical and translation invariant with finite range R, where R is the radius of the neighbourhood V(0) of 0, there exists T defined on [21,1] 2 and valued in [21,1] such as:

w ij ~T i 1 {j 1 ð Þh=R, i 2 {j 2 ð Þh=R ð Þ ,
where h is a strictly positive real number, T has as mean value m, Let us denote now by f a continuous regularized version of the Heaviside function (like the arc-tan) and let us take

F~f {1 , a u ð Þ~lR 4 S h 2 F' u ð Þ À Á , b u ð Þ~{F u ð ÞzlR 2 h 2 mu À Á F' u ð Þ
then the reaction-diffusion operator defined by:

Lu=Lt~a u ð ÞDuzb u ð Þ
has a natural discretization corresponding to the neural network above, by identifying x i (t) and u(ih,t) and by remarking that the neural network system has the same asymptotic behaviour as the differential system:

dx i t ð Þ=dt~l X w ij x j t ð Þ{H x i t ð Þ ð Þ . H 0 H x i t ð Þ ð Þ ð Þ ,
when l is sufficiently large. In [START_REF] Cottet | Mode `les de re ´action-diffusion pour des re ´seaux de neurones stochastiques et de ´terministes[END_REF], it is shown that, for adapted values of R, homogeneous in grey, 1D objects can be enhanced in a heterogeneous environment, in the same way as for a neural network system. In [START_REF] Cottet | Image segmentation using snake-splines and anisotropic diffusion operators[END_REF] and in Figure 11 (right), the same proof is given for 2D objects like the internal cavities of the heart, where a snakes splines procedure is used after contrasting. 3) Proposal for a new image reaction-diffusion-chemotaxis operator. In order to have, like for the previous operator, the final treated image as asymptotic of a differential operator, we propose to consider the grey level u as a chemotactic substrate concentration consumed by animals whose concentration will be denoted by v [START_REF] Demongeot | Multiscale Modeling and Imaging: The Challenges of Biocomplexity[END_REF][START_REF] Demongeot | Dynamical Systems and Bio-medical Image Processing[END_REF]144]. The principle of this method consists in locating initially a uniform concentration v(0) of animals on the initial grey level image u(0) or on its boundary: the substrate u can diffuse with a term eDu and is consumed with a saturation rate equal to: 2Kuv/(u+k); the animal concentration v can diffuse attracted by the substrate with the term DDv, is submitted to a drift in the direction of substrate peaks with the chemotactic term -xdiv(vgradu) and increases (because of the reproduction) with the term K'uv(u+k'). Let us remark that the two first terms ruling the animal motion can be replaced, if we do not want to introduce a drift term, by an attractiondiffusion term like:

D L 2 v Lx 2 : Lu=LxzL 2 v Ly 2 : Lu=Ly À Á
The corresponding differential partial derivative operator is then given by:

Lu=Lt~eDu{Kuv= uzk ð Þ, Lv=Lt~DDv{xdiv vgradu ð Þ zK 0 uv uzk 0 ð Þ
or by the following PDE:

Lu=Lt~eDu{Kuv= uzk ð Þ Lv Lt~D L 2 v Lx 2 : Lu=LxzL 2 v Ly 2 : Lu=Ly À Á zK 0 uv= uzk 0 ð Þ
In the two cases above, the asymptote of u is 0 and the asymptotes of v give the ''treated image''. The corresponding image processing leads to a contrast enhancement before segmentation: in Figure 7, we can see the initial image on the bottom left and the contrasted one on the bottom right. The contours have been then easily obtained by applying a snakes splines procedure [START_REF] Leitner | Dynamic Segmentation: Finding the Edge with Snake-Splines[END_REF][START_REF] Cohen | On active contour models and balloons[END_REF]. If we are adding to the second equation of the differential system a Dupin term like Kv/Du, we will encourage animals to follow Dupin lines, i.e. inflexion curves, which is very suitable for a grey anticlines segmentation (for example in vessels segmentation).

Conclusion

The neuro-mimetic lateral inhibition mechanism and the set valued snakes-like flows allow the generation of various image processing methods (essentially contrast enhancement and contouring). We have given numerous applications of this methodological approach in image processing essentially dedicated to medical imaging and surgical robotics. Further both theoretical and numerical studies have to be completed, in order to show the utility of these new tools in morphogenesis modelling, allowing to generate artificial objects of biological and/or medical interest (like cells, tissues, organs) by using the same operators as for recognizing them in a real image. We conjecture that the spatial information about anatomic organs obtained from the biomimetic image-processing methods, has to do with the morphogens localization, which results from the morphogenetic processes creating these organs combining robust genetic regulatory networks [145,146] ruling their metabolic reactions and cell proliferation, with classical diffusion [147] of morphogens inside their tissues. In particular, the main patterns observed during the embryonic formation can be found in the biomimetic processing of the images by the final adult organ. INTRODUCTION Recent advances in the knowledge about the regulatory genetic networks, in particular in the field of the control by micro-RNAs of central functions like the progression in the cell cycle have shown that the architecture of genetic interactions graphs is often made of a central core which is a circuit between genes, surrounded by an upper tree with multiple roots whose leaves are the genes of the circuit, and of a pendant tree whose roots are these genes and the leaves are the final control genes ruling the progression of the cell in its cycle.

II. DISTANCES BETWEEN GRAPHS

Mathematical distances between trees

Let G=(V,E), where V is a set of vertices -or nodes, like genes i,j,... -and E a set of weighted edges -or interactions like w ij , where w ij >0 (resp. w ij >0) means that the protein expressed by the gene j induces (resp. represses) the activity of the gene i. G is a finite undirected graph, supposed to be without selfcircuits, but containing circuits (or loops, i.e. paths between genes i and itself passing through at least one intermediate gene) and multirooted trees and it is possible to calculate, if G is random, the expectation of the circuits and trees numbers [2]. Let us suppose that there is only a circuit (genes G 1 to G 8 on Figure 1) with an upstream and a downstream tree, called respectively in the following uptree (genes G 9 , G 10 ) and downtree (genes G 11 , G 12 ). Figure 1. Circuit (genes G1 to G8), uptree (G9, G10) and downtree (G11, G12)

Within the framework of this study, we will now introduce distances between trees and circuits: A. Nodal distances between trees [3] Two multirooted phylogenetic trees T, T' on the same set S of genes are isomorphic with a nodal distance equal to 0, if, and only if, for every pair of genes (i,j), the distances between i and j (i.e. the length of the minimal paths between i and j, denoted L(i,j) and L'(i,j) in the undirected trees corresponding to T and T') are the same in T and T'. Any distance d between the matrices L and L' is called nodal distance between T and T'. If T' is defined on a subset B of the genes set of T, the distance between T' and T is equal to those between T' and the restriction of T to B, the same paths length as in T being kept.

G 1 G 2 G 8 G 3 G 7 G 6 G 5 G 4 G 9 G 10

B. p-uplet distance between trees [4]

The quartet distance is the number of p-uplets -sub-trees induced by p leaves -that differs between the trees. The other distances between trees are more complicated; we can use the pq-gram distance [5], the general edit distance between trees [6] or the edge shift distance [7] between trees.

Mathematical distances between circuits

The nodal distances between circuits have the same definition as above in A.

III. EVOLUTION OF THE REGULATORY NETWORKS: THE EXAMPLE OF THE CELL CYCLE CONTROL

The cell cycle control network

The cell cycle is controlled in all realms where prokaryotic or eukaryotic cells are present, i.e. practically everywhere in life systems. In mammals, the first spheric symmetry breakdown of the ovule is due to an irreversible reaction-diffusion process organizing spatial waves of calcium [8], caused by the shock of the spermatozoid head on the pellucid membrane of the ovule: it leads after equiploid cell divisions (the others leading in general to abortive mitoses [9]) occurring on the spatial nodes of calcium on which diffuse many cytoskeleton components as tubuline [10,15], giving to the embryo a cylindric cephalo-caudal structure (Figure 2). The progressive apparition during the evolution of many upstream controllers until the mammals microRNAs or p27 [16] and p53 [17], by keeping nodal distance at the value 0 between Drosophila melanogaster and mammals, thus giving a robust [18][19][20][21] Rb-E2F control (Figure 3).

IV. GASTRULATION

The initial embryonic growth of an organ defined inside its biological boundary is exponential between blastula and gastrula stages [22]. The cell divisions are ruled by the following constraints limiting the cells and embryo growth duration:

1) Thom's laws [23] -A cell divides when its ratio S/V (Surface divided by the Volume inside the cell boundary) reaches a critical value -An organ stops its growth when its ratio S/V reaches a critical value

2) Hayflick's laws [24,25] -A cell lineage in a cell culture stops before a fixed finite number of cell divisions; this number is called the Hayflick's limit k and goes from 7 in C. elegans, to 15 in mice and 110 in Galapagos tortoise. If the embryonic growth is stopped by the Thom's constraint, then we will denote by h the number of divisions remaining for the adult phase. If N e denotes the number of cells observed at the end of the embryonic growth, then h ≈ k -Log 2 N e (by neglecting the apoptosis during the embryonic phase) -A cell population in a stationary organ has an apoptosis rate which compensates the cell divisions.

On Figure 4, we have in blue the initial pattern and in black the pattern after mechanical perturbations due to the first bottle (or flask) cell differentiation : the most basal cell is supposed to be the first to differentiate and to present an apical constriction due to an overconcentration of myosin at this cell pole. Then each cell evolves following the mechanical laws: 1) the sum of the orthogonal forces exerted on the cell walls is equal to its mass (proportional to its inner area) times its orthogonal (to the wall) aceleration. The external forces are the resultant of the extracellular (Pext) and intracellular (Pint) pressure, the contact forces being exerted by the neighbouring cells (sharing a part or the whole of the wall) exerted by other cells and the internal force is due to the cell pressure exerted via the cyto-skeleton and the plasmic membrane [26]. Each force is equal to a coefficient (e.g. the physical pressure) times the length of the wall on which it is exerted 2) If we suppose that the initial cell configuration is in an equilibrium state, we calculate an admissible set of parameters values respecting this equilibrium (Figure 4)

3) Then, we leave the cell system evolve depending on the energetic balance ruling the cytoskeleton apical polymerization [27][28][29] controlled by a specific genetic regulatory network comprizing essentially concertina (cta), actin, myosin, Rho and RhoGEF genes (Figure 5), by choosing a small time step, by updating sequentially each cell and by calculating their displacements respecting the no-overlapping rule (Figure 6). At each step we update the cell common walls by supposing that cell contacts are close, ensured between cells by cadherins and gap junctions [30], and with the extracellular matrix by integrins and adhesins. Cell motion involves a change in its inner area: we suppose that growth occurs where internal forces Fint are larger than external ones Fext (cell has to be stretched), as for constrained growth in continuous media mechanics. When the external forces dominate, the cell is supposed compressible and can be constricted. The cells divide longitudinally or laterally when their ratio perimeter over area is too small (Thom's first law). The revolution symmetry breaking consists in making two cylinders inside the embryo cylinder, one coming from the ventral furrow and giving the digestive tube (this phase is called gastrulation) and the other coming from the dorsal furrow and giving the neural tube (this phase is called neurulation). The first gastrulation step consist in an apical concentration overfluctuation of one of the cytoskeleton components (myosin, actin, tubulin,...) or one of the enzymes (ADenylate Kinase -ADK -or Nucleoside Diphosphate Kinase -NDK) or one of the energy molecules (ATP, GTP) involved. This apical overconcentration diffuses and reaches the extremities [31,32] of the cylindric portion of the embryo at which gap junctions are less important with the cells of the "hemispheric" terminations [30]. Then two first bottle cells can appear at the two extremities of the diffusion line (Figure 7), then this first contraction can propagate until the center of the cylindric part as noticed in [31,32] during the 6 first hours of the gastrulation in Drosophila melanogaster. After this first phase of the ventral furrow formation, for regularizing the tensegrity forces, based on a synergy between balanced tension and compression components of cells, first divisions occur from extremities to center by applying the first Thom's cell law. We can simulate such a process on a transversal slice of the embryo and reproduce the invagination of the gastrulation from both the contraction and consecutively the proliferation phases (Figure 8). The application of the following automaton reproduces rather well the reality (Figure 9) :

1) each cell i has a state at time t represented by 2 bits (x i (t),y i (t)), where x i equals 1 (resp. 0) when it is contracted (resp. normal) and y i equals 1 (resp. 0) in presence of an overconcentration of at least one actor (called generically gastrulation morphogen) of the bottle cell differentiation 2) the initial state of all the cells is initialize at the value (0,0) and any cell state cannot change if a cell is surrounded by 4 cells in the same state 3) in a cell i 0 from the cylindric part of the embryo (in blue on It has to be implemented on a 3D mesh in order to become more realistic and to be able to perform a final closure of the tube from the invaginated streak. A way to obtain it after the morphogen diffusion phase on the 2D manifold corresponding to the embryonic surface consists in using a diffusion reaction Rho Cta Fog Twi D Actin/Myosin GDI equation with at least a Laplacian and a logistic term [26], and moving boundary conditions depending on the new frontier obtained after cell growth and proliferation.

Figure 7 (from [31]). Formation of the ventral furrow in Drosophila melanogaster with the propagation of the contraction wave from embryo extremities to center. Last picture shows the concentration of myosin (green) Figure 8. Left: in blue, the progressive invagination due to the first bottle cell differentiation followed by its consecutive divisions. Left (in green): gastrulation steps with the increase of the intracellular myosin level (http://www.molbio1.princeton.edu/wieschaus/). Right: gastrulation steps from the apical constriction of bottle cells to the constitution of the tube.

On Figure 10, we see the simulation of a pathologic proliferation (after branch lopping or parastic attack) compensating in regions of accumulation of the morphogen auxin. Because the ordinary convexity of peripheral cells fed in the tree The proliferation is supposed to be more important in the concave parts of the tree (the concavity is defined with respect to the sap canals) and less (resp. more) important in the zones of least diffusion (resp. accumulation) of the plant morphogen (auxin) Branch lopping

Cronartium ribicola

Pinus radiata

Auxin accumulation by the sap canals (equivalent to the external feeding of the cylindric embryo by the amniotic liquid in mammals or by the egg vitellus in insects) gives a normal proliferation, the cambial cells blocked in their growth by a pathologic process (branch lopping or parasitic attack) will have an adverse surface to volume ratio in 3D (resp. perimeter to area ratio in our 2D simulations) and then will divide for increasing this ratio [33]. The essential difference between trees and embryo is due to the fact that the "concave" bottle cells have a permanent shape due to their differentiation and hence their divisions will stop only after internalization of the primitive digestive tube, but the peripheral tree cambial cells are concave only due to the differential proliferation of the zones non involved in growth perturbation like parasitic attack or branch lopping. Then the responses are different: in embryo, the bottle cells proliferation leads to a revolution symmetry breaking by creating a small cylinder inside the primitive one; on the contrary, the cambial cells proliferation leads to a revolution symmetry breaking due to a overproduction of cells in the tree vertical concavity causing an absence of restitutio neither in nor ad integrum (compensatory process called also "apocatagenesis"), but a pathologic S shaped tree (see Figure 10 and [34][35][36][37][38]). The end of the cell and tissue growth stops correspond to the end of the morphogenesis process: this can be observed when the second Thom's tissue law is applicable (the surface to volume ratio of an organ becoming adverse) or when the organ is completely covered by an anatomic boundary like an aponeurosis made of fibrous cells or an autoassemblage of extracellular ingredients. In both these cases, a couple of morphogens acting often simultaneously in opposite (e.g. a couple of activator and inhibitor like BMP-7 and BMP-2 in feather morphogenesis in the chicken [39][40][41] can induce the chemotactic motion of fibroblats or the biosynthesis of the elements constituting the auto-assemblage (like proteins and phospholipids). The fact that for a certain value of their viscosity ratio, the morphogens can coexist for a relatively long time in a precise location can greatly favor the birth of anatomic organ boundaries.

V. CONCLUSION

The deep understanding of complex morphogenetic processes like the gastrulation needs the explicitation of multiple interconnected levels. In this paper, we dealt with some of them:

1) the genetic level comprises regulatory networks which have evolve until robust architectures, weakly sensitive to boundary changes (in the graph sense, that is external transcription factors like proteins coming from other metabolisms or traduction factors like RNE-binding oligopetides and micro-Rnas) or to endogeneous perturbation (over or underexpression of a gene leading to an up-or down-regulation of a protein directly involved in the morphogenesis). Some networks are specific even they have general characteristics, like a great proportion of inhibitions among their genes interactions [42][43][44][45] and some other networks are generic, like those controlling the cell-cycle [46] or the nycthemeral rhythms associated with the morphogenesis [47] 2) the cell level is essentially depending on the entrance or not in the cell cycle. The cells blocked in G0 phase can serve secondarily when the cohorte of divinding cells will finished their divisions authorized by the Hayflick's rules, then will stop their lineages and consequently desinhibit a new cohorte of cells in resting state [22]. The cell growth is also under the dependence on the first Thom's law which obliges the cell to respect a surface to volume ratio threshold

3) the tissue level is organized by morphogens acting as chemo-effectors (attractant or repulsant): when they share privilegiated zones in which they present a small diffusion, they are able to favor the coexistence of cells or bricks of future layers, membranes or aponeuroses needed for building the organ anatomic architecture from the functional core to the protective boundary.

The association of this level has to influence the experimental protocol in which only the crucial (needed for understand the mechanisms) as well as critical (causing deep changes in he shape or in the dynamics observed during the morphogenesis) variables [48]. A model driven acquisition process will also allow an easy extraction of the model architecture like the interaction graph of a genetic network [49], and after the study of the mathematical model corresponding to this architecture (well posed direct and inverse problem) [49][50][51][52].

We given in this paper some examples (centered around the gastrulation process) concerning these different steps from experimental acquisition of pertinent data until the interpretation in a mathematical framework of the dynamic or geometric features of the observed forms and functions. We will in the future increase the spatial and temporal resolution of the data by sampling in 3D with a good precision in space (e.g. by using confocal or biphotonic microscopic information), in time (e.g. by using the cinemicroscopy) and in function (e.g. by using the Raman, vibrational, infrared,..., spectromicroscopy). We will study physiologic morphogenetic processes occurring in the development as well as dysgenetic processes observable during ageing, and also in between regenerative processes (healing or compensation after loss of organ due to an infarctus or to a fracture or an ablation) which can lead to a partial reparation (restitutio ad integrum) until a complete rehabilitation, both anatomic and functional (restitutio in integrum or apocatagenesis). INTRODUCTION During the morphogenesis of the gastrula, the second step of the embryogenesis after the blastula stage, we observe in the majority of the animal development dynamics, the following phenomena: i) the cell motion is partly guided by chemotaxis, in order to supply their substrate demand, and also to respect the epigenetic architecture ruled by morphogens, ii) the cell shape is due to a constriction controlled by cell differentiation and iii) the final gastric tube is obtained from cell proliferation relaxing the forces exerted on the cell plasmic membrane and optimizing the cellular "nutritive Area / inner Volume" (A/V) ratio. Concerning the differentiation process, some cells of the embryo start to take the shape of a bottle (bottle or flask cells), decreasing the surface at the interface with their nourishment fluid. At the end of the gastrulation, these bottle cells start to divide and grow, increasing their A/V ratio. In this paper, we are attempting firstly to follow by using image processing techniques the first phase of the gastrulation made of pure cell Manuscript received 15th November 2010. 1 Laboratory AGIM, FRE CNRS University J. Fourier of Grenoble, Faculty of Medicine, 38700 La Tronche, France 2 Laboratory TIMC-IMAG, UMR 5525 CNRS University J. Fourier of Grenoble, Faculty of Medicine, 38700 La Tronche, France *corresponding author: Jacques.Demongeot@imag.fr, e-mails: tayyab.elite@gmail.com, Athanasios.Lontos@imag.fr, Emmanuel.Promayon@imag.fr motion causing the invagination without any proliferation, and secondly to show that the cell constriction results in a streak starting on the two extremities of the embryo and propagating secondly to its central part. In order to simulate the initiation of this phenomenon, we use exclusively laws of physics and mechanics. For this purpose, we have created a threedimensional biomechanical model consisting of a group of cells, forming a structure with two areas: a cylindrical rigid area representing the main embryo body and two curved hemispheric areas constituting the embryo extremities. Each cell is modelled by a set of Newtonian contour particles defining the surface mesh and an inner particle, located at the geometric centre of the cell, mimicking the centrosome. Elastic forces are used to model the rigidifying effect of the tubulin and actin cytoskeleton, and contractile forces to model the action of the myosin fibers. In addition, we have modelled the role of the cadherins by connecting contour particles of adjacent cells. Finally, in order to ensure the symmetry from a partial mesh representing only a part of the embryo, we have applied boundary conditions to the particles located at the lateral extremity edges of the simulated structure. This biomechanical model shows that the inner folding starts at the curved area and then propagates to the rigid one, in accordance with the in vivo process. To compare these results with in vivo experiments, we have monitored the displacement of the centrosome and the cell A/V ratio. Our numerical experiments made during the time lag observed through the microscopic imaging without proliferation suggest that our model adequately simulates thge in-vivo cell behaviour. In a first part, we present image processing techniques and results obtained by applying them on gastrulation microscopic recording of Drosophila melanogaster embryo from [2][3][4][5][6]. In a second part, we describe the biomechanical model of streak formation and the third Section will be devoted to the presentation of numerical simulations confronted to real images of the first invagination stage of the gastrulation. care, if they occur, our regions of interest. Image enhancement techniques need to be applied on the images to get welldefined cell contours. The available images depict actinmyosin networks in colour. Here, we are focussing more on cell boundaries in a particular region, in order to have an idea of cell density; and to observe cells whether they divide in those regions. After evaluating the colour histogram of the image, we filter out the noise by selecting an appropriate colour threshold value. In absence of noise, cell boundaries are easily visualized and hence cells could be tracked. Cell boundaries are further enhanced by applying contrast algorithms on the colour-filtered image again by using simple threshold values after studying the gray-histogram. Cell proliferation results in increase in volume, thereby causing cells to recurrently contract then stabilize. Older cells are more vulnerable to shrink themselves as a result of forces from the neighbouring cells, due to proliferation at the other end. Cell density, as expected, would increase enormously in such regions with respect to other regions (Figure 1), pushing them due to proliferation. Cell density increasing enormously in such regions with respect to other regions which could lead to invagination, close to the point of maximum cell density.

II. IMAGE PROCESSING

Results obtained from cell contouring and counting

Available series of images were processed to have an idea of change of cell density in the region where a constriction appeared at a later stage. After following the image processing steps explained in 2.1., we were able to have an exact visual idea of how the cells move and how changes the cell density in the region where invagination occurs experimentally, and also to observe if there is any division of cell in that region before the development of constriction. This included removal of noise from the images and then applying contrast and contour algorithms, to have an exact visual idea of the number of cells in that region. All cells could be individually tracked, and change of their forms are clearly observed. However, no division of cells was observed in the region where in the later stage, constriction appeared. Instead, cells get squeezed and cell density increased in that particular region before we could see two points in-line, across which this phenomenon occurred, leading finally to the development of a constriction.

I. THE BIOMECHANICAL MODEL

The embryo architecture

The embryo cylindric cephalo-caudal structure (Figure 2) has two zones of cylindricity breaking, the two hemispheric extremities and the cells at the boundary between the two geometries presents a curvature due to a relaxation of the internal rigid properties of their cytoskeleton. Tha appearnce of the gastrulation streak occurs at their level and after propagates to the central cylindric part of the embryo. Taking into account the differential contractility of these boundary cells constitutes the model core and will be implemented over a 3D mesh representing the external embryo surface. 1) the sum of the orthogonal forces exerted on the cell walls is equal to its mass (proportional to its inner area) times its orthogonal (to the wall) acceleration. The external forces are the resultant of the extracellular (Pext) and intracellular (Pint) pressure, the contact forces being exerted by the neighbouring cells (sharing a part or the whole of the wall) exerted by other cells and the internal force is due to the cell pressure exerted via the cyto-skeleton and the plasmic membrane [7]. Each force is equal to a coefficient (e.g. the physical pressure) times the length of the wall on which it is exerted 2) If we suppose that the initial cell configuration is in an equilibrium state, we calculate an admissible set of parameters values respecting this equilibrium 3) Then, we leave the cell system evolve depending on the energetic balance ruling the cytoskeleton apical polymerization [4][5][6] controlled by a specific genetic regulatory network comprizing essentially concertina (cta), actin, myosin, Rho and RhoGEF genes [7][8][9][10], by choosing a small time step, by updating sequentially each cell and by calculating their displacements respecting the no-overlapping rule. At each step we update the cell common walls by supposing that cell contacts are close, ensured between cells by cadherins and gap junctions [11,12], and with the extracellular matrix by integrins and adhesins. Cell motion involves a change in its inner area: we suppose that growth occurs where internal forces Fint are larger than external ones Fext (cell has to be stretched), as for constrained growth in continuous media mechanics. When the external forces dominate, the cell is supposed compressible and can be constricted. The cells divide longitudinally or laterally when their ratio perimeter over area is too small [13][14][15][16][17]. The revolution symmetry breaking consists in making two cylinders inside the embryo cylinder, one coming from the ventral furrow and giving the digestive tube (this phase is called gastrulation) and the other coming from the dorsal furrow and giving the neural tube (this phase is called neurulation). The first gastrulation step consist in an apical concentration overfluctuation of one of the cytoskeleton components (myosin, actin, tubulin,...) or one of the enzymes (ADenylate Kinase -ADK -or Nucleoside Diphosphate Kinase -NDK) or one of the energy molecules (ATP, GTP) involved [18][19][20][21][22][23]. This apical overconcentration diffuses and reaches the extremities [11,12] of the cylindric portion of the embryo at which gap junctions are less important with the cells of the "hemispheric" terminations. Then two first bottle cells can appear at the two extremities of the diffusion line, then this first contraction can propagate until the center of the cylindric part as noticed in [11,12] during the 6 first hours of the gastrulation in Drosophila melanogaster. After this first phase of the ventral furrow formation, for regularizing the tensegrity forces, based on a synergy between balanced tension and compression components of cells, first divisions occur from extremities to center by applying the first Thom's cell law. We can simulate such a process on a transversal slice of the embryo and reproduce the invagination of the gastrulation from both the contraction and consecutively the proliferation phases.

The end of the cell and tissue growth stops correspond to the end of the morphogenesis process: this can be observed when the second Thom's tissue law is applicable (the surface to volume ratio of an organ becoming adverse) or when the organ is completely covered by an anatomic boundary like an aponeurosis made of fibrous cells or an autoassemblage of extracellular ingredients.

In both these cases, a couple of morphogens acting often simultaneously in opposite (e.g. a couple of activator and inhibitor like BMP-7 and BMP-2 in feather morphogenesis in the chicken [24][25][26][27] can induce the chemotactic motion of fibroblats or the biosynthesis of the elements constituting the auto-assemblage (like proteins and phospholipids). The fact that for a certain value of their viscosity ratio, the morphogens can coexist for a relatively long time in a precise location can greatly favor the birth of anatomic organ boundaries.

II. NUMERICAL SIMULATIONS

The myosin is supposed to diffuse in all directions from a cell in which an excess of myosin is synthesized caused by random fluctuations over-expressing its gene, notably along a directrix of the cylinder constituting the body of the embryo of Drosophila melanogaster. This directrix represents the shortest path until the hemispheric extremities of the embryo on the boundary of which cells change of curvature (yellow on Figure 3). The whole model mixes a Multi-Agents Model (MAM) responsible for cell growth, migration and proliferation with a Reaction-Diffusion Partial Differential Equation (RD-PDE) for myosin dispatching, mainly in charge of cell contraction. The results of simulation of the hybrid model are given on Figures 3 and4 showing the same phenomena as those observed in experiments: the progressive invagination starts at the extremities of the embryo and after propagates to the central cylindric part of the embryo. The final step of the gastrulation needs proliferation in order to close the internal of the cell differentiation in bottle (or flask) cells contracted at their apical extremity, provoking during the first minutes of gastrulation a reorganization of the superficial cell layer of the embryo without division, leading to the formation of a streak. Resulting invagination starts experimentally at the 2 extremities of the embryo and propagates to its central part.cylinder which will give birth to the intestinal tube of the adult animal. By following the progressive migration in embryo depth of the centrosomes it appears that the run is faster and deeper for curved cells at the extremities of the embryo than for central cells (Figure 5). The saturation curves representing this displacement behave like the curves representing the evolution of the cell diameter under progressively increasing forces applied externally to the cell, which is a way to induce cell contractions (possibly periodic [30]) similar to those due to myosin (Figure 6) [31,32]. Some divisions of the observed bottle (or flask) cells whose differentiation is due to myosin gene over-expression, suffice to end the gastrulation process if they occur at critical locations as streak lips (red star on Figure 3), located at the boundary of the zero-diffusion domain both observed and simulated for myosin (in yellow on Figure 3). 6: Simulated virtual spherical red blood cell (RBC) suspended in an hypotonic solution. Optical tweezers double trap is simulated by exerting locally a force Fs on two opposite nodes of the cell object contour (upper insert). The variation with load of the cell object diameter D(F) (in µm) in a plane perpendicular to the loading direction is simulated and compared to experimental data published by [31] (black squares). With appropriate scaling of the force, it is possible to adjust the elasticity modulus such that experimental mechanical response of RBC is nicely fitted in the linear elastic regime (red). Increasing the elasticity modulus induces a stiffer response which qualitatively reproduces the departure from the linear regime at larger traction forces (violet).

V. CONCLUSION

We given in this paper some examples (centred around the gastrulation process) concerning different steps of morphogenesis modelling from experimental acquisition of pertinent data until the interpretation in a mathematical framework of the dynamic or the geometric features of observed forms and functions.

We will in the future increase the spatial and temporal resolution of the data by sampling in 3D with a good precision in space (e.g. by using confocal or biphotonic microscopic information), in time (e.g. by using the cine-microscopy) and in function (e.g. by using the Raman vibrational, infrared,..., spectro-microscopy), in order to detect precisely the frontiers between the successive phases of the morphogenetic process of revolution symmetry breaking in the embryo, i.e., i) the morphogen diffusion, ii) the cell migration, iii) the bottle cell differentiation, iv) the streak contraction and v) cell the proliferation ensuring vi) the tube closure.

  détecter et d'observer cette mitose, tout en proposant une intervention interactive humaine en cas d'échec, afin d'éliminer toute fausse inférence. a branch of science which deals with the computer-based analysis of large biological data sets. It is a developing field helping in collecting, linking, and manipulating different types of biological information to discover new biological insight.

Figure 1 - 1 :

 11 Figure 1-1: (a) Schematic representation of the retina and cellular populations (b) SEM of cones and rods ([5]).

Figure 1 - 2 :

 12 Figure 1-2: Cone Rod Density ([9]).

Figure 1 - 3 :

 13 Figure 1-3: Distribution of rods and cones in fovea (left) and in the periphery (right) ([10]).

Figure 1 - 4 :

 14 Figure 1-4: 3D shadow projections of confocal images of the fluorescent cones. Scale bar: 10µm ([11]).

Figure 1 - 5 :

 15 Figure 1-5: Rod and cone cell populations during the degeneration of rd1 mouse (processing done by T. Léveillard, U592 INSERM).

Figure 1 - 6 :

 16 Figure 1-6: Cone mosaic as revealed by anti-cone opsin immunostaining (cones are absent in the circles) ([18]).

5 Figure 1 - 7 :

 517 Figure 1-7: Cone-rod distributions in the two opposite quadrants NW and SE.

  We will not present here definitive conclusions about estimations of w and homogeneity tests, because this part of the research program has been stopped 2 years ago by a dramatic event in the U592 INSERM (http://www.20minutes.fr/paris/296891-Paris-Unblesse-dans-l-incendie-a-l-hopital-des-Quinze-Vingts-a-Paris.php) interrupting suddenly the retina image production.

  restoration and maintenance of retinal function and potential treating of RP patients. So, we would like to take advantage of digitalized images from expression cloning experiments to look for factors directing cell differentiation. In the experiments, cone-enriched retinal cultures from chicken embryos are treated with pools of 100 clones. Each pool of 100 clones from the retinal expression library was tested for its ability to maintain cone viability (live/dead kit from Invitrogen) in 8 independent wells.

Figure 1 - 7 :

 17 Figure 1-7: Cell morphology in cone-enriched cultures (image obtained by T. Léveillard, U592 INSERM).

  been injected in the eye of the rd1 mouse, allowing to follow the differentiation and the function (electro-retinogram) of the injected mice. These experiments conducted in conjunction with Rd-CVF tests could lead to a biotherapy of RP, involving RdCVF and the novel factor in order to maintain respectively cell viability and cell function. In the future an application could be possible on the human retina, if a non invasive procedure (like the fluorescence of non-toxic membrane potential dependent optical dyes) could be performed on the still functional retinal cells, which would indicate the localization of the still living and functional cells during degeneration.

Figure 2 - 1 :

 21 Figure 2-1: Confocal microscopy images from the whole mouse pathologic retina with 236 imagettes ([13]).

Figure 2 - 2 :

 22 Figure 2-2: (a) Confocal microscopy images from the whole mouse pathologic retina with 236 imagettes [13], (b) particular imagette, and (c) its skeletonization. Results after (d) Laplacian filtering, (e) Sobel Filtering, and (f) 5 x 5 Grey Histogram Thresholding.

Figure 2 - 2 :

 22 Figure 2-2: Detected transitions represented by crosses ([1]).

Figure 2 - 3 :

 23 Figure 2-3: Left image shows light grains separated by darker boundaries. Right image shows a mix of light and dark phases ([1]).

1 . 3 .

 13 Circumscription of a circular area of interest within the image to prevent any anisotropy bias induced by the rectangular shape of the image frame, 2. Detection of the transitions inside the mask for a specified direction as explained in Paragraph 2.2.1 and measurement of intercept lengths, Repetition of step 2 by rotating the image and re-sampling the pixels by means of a linear transform.

Figure 2 - 4 :

 24 Figure 2-4: Local histograms and dispatching between internal (black), external (blue) and frontier (red) types.

Figure 3 - 1 :

 31 Figure 3-1: Retinal Explants ([13]).

Figure 3 - 2 :

 32 Figure 3-2: A sample Imagette from Retinal Explants in Figure 3-1.

3 . 2 SEGMENTATION

 32 Segmentation is a technique to reduce the image only to useful information. Segmentation divides image into regions or objects of interest. The idea is to isolate these objects of interest which are cones in our case. The algorithms of segmenting images are based on discontinuity or similarity in pixel values. In discontinuity-based segmentation, an image is partitioned on the basis of abrupt changes in intensity. Thresholding, region-growing are examples of similarity-based segmentations.

Figure 3 - 3 :

 33 Figure 3-3: Sobel with 3 x 3 and 5 x 5 Masks.

Figure 3 - 4 :

 34 Figure 3-4: Results after application of Sobel Filter (Left) using 3 x3 Mask, and (Right) using 5 x 5 Mask.

  employed for edge detection. Before applying the filter to the image, the Gaussian imagesmoothing is done to reduce the sensitivity to noise. Typical 3x3 and 5x5 masks being utilized in this filtering technique are shown in Figure3-5The results of the two masks are shown in Figure3-6.

Figure 3 - 5 :

 35 Figure 3-5: Laplacian Filter with 3 x 3 and 5 x 5 Masks.

Figure 3 - 6 :

 36 Figure 3-6: Results after application of Laplacian Filter (Left) using 3x3 Mask, and (Right) using 5 x 5 Mask.

Figure 3 - 7 :

 37 Figure 3-7: Histogram of the scaled image.

Figure 3 -

 3 Figure 3-8 depicts some histograms of adjacent masks. From these histograms, it is very easy to decide whether the pixel belongs to inside or outside of the object and also if it belongs to the frontier of the object.

Figure 3 - 8 :

 38 Figure 3-8: Local histograms of three different cases in Histogram-Thresholding.

Figure 3 - 9 :

 39 Figure 3-9: Segmentation of cones using Histogram-Thresholding.

Figure 3 - 10 :

 310 Figure 3-10: Conversion of Segmented Cones to Solids.

Figure 3 - 11 :

 311 Figure 3-11: Segmentation of cones for complete retinal explants.

Figure 3 - 12 :

 312 Figure 3-12: Representation of 98 objects from the connected components of the internal pixels set.

Figure 3 - 13 :

 313 Figure 3-13: Horizontal intercept histograms in the 4 quadrants LS (dark blue), RI (light blue), RS (violet) and LI (yellow) (processing done by Y. Usson, TIMC-IMAG).

Figure 3 - 14 :

 314 Figure 3-14: Horizontal intercept distance Gaussian distribution in the 4 quadrants LS (dark blue), RI (light blue), RS (violet) and LI (yellow) ([26]).

Figure 4 - 1 :

 41 Figure 4-1: Analysis of a Healthy Region.

Figure 4 . 2 :

 42 Figure 4.2: Analysis of a Region with Retinal Degeneracy.

Figure 4 - 3 :

 43 Figure 4-3: Analysis of a Region affected partially by Retinal Degeneracy.

4 . 2 PERSPECTIVESFigure 4 - 4 :

 4244 Figure 4-4: Idea of diagnosis of sick portion following the image processing.

3. 1 . 66 CHAPTER

 166 THE EMBRYO ARCHITECTURE.......................................... .................................................62 3.2. THE MODEL (A. LONTOS)............................................................................. ....................62CHAPTER NUMERICAL SIMULATIONS ……………………………………………………………….Morphogenesis is the biological process that causes an organism to develop its shape. It is one of three fundamental aspects of developmental biology along with the control of cell growth and cellular differentiation. The process controls the organized spatial distribution of cells during the embryonic development of an organism.Gastrulation is a phase early in the development of animal embryos, during which the morphology of the embryo is dramatically restructured by cell migration. Gastrulation varies in different phyla. Gastrulation is followed by organogenesis, when individual organs develop within the newly formed germ layers.

Figure 1 - 1 :

 11 Figure 1-1: Gastrulation of a diploblast: The formation of germ layers from a (1) blastula to a (2)gastrula. Some of ectoderm cells (orange) move inward, forming the endoderm (red)[33].

Figure 1 - 2 :

 12 Figure 1-2: Surface view of rabbit embryo ([33]).

  use exclusively laws of physics and mechanics. For this purpose, the authors have created a three dimensional biomechanical model consisting of a group of cells, forming a structure with two areas: a cylindrical rigid area representing the main embryo body and two curved hemispheric areas constituting the embryo extremities. Each cell is modelled by a set of Newtonian contour particles defining the surface mesh and an inner particle, located at the geometric centre of the cell, mimicking the centrosome. Elastic forces are used to model the rigidifying effect of the tubulin and actin cytoskeleton, and contractile forces to model the action of the myosin fibres. In addition, the authors have modelled the role of the cadherins by connecting contour particles of adjacent cells. Finally, in order to ensure the symmetry from a partial mesh representing only a part of the embryo, the authors have applied boundary conditions to the particles located at the lateral extremity edges of the simulated structure. This biomechanical model shows that the inner folding starts at the curved area and then propagates to the rigid one, in accordance with the in vivo process. To compare these results with in vivo experiments, the authors have monitored the displacement of the centrosome and the cell A/V ratio. Numerical experiments made during the time lag observed through the microscopic imaging without proliferation suggest that the model adequately simulates the in-vivo cell behaviour. In a first part, the authors present image processing techniques and results obtained by applying them on gastrulation microscopic recording of Drosophila melanogaster embryo from[2][3][4][5][6]. In a second part, we describe the biomechanical model of streak formation and the third section will be devoted to the presentation of numerical simulations confronted to real images of the first invagination stage of the gastrulation.

Figure 1 - 3 :

 13 Figure 1-3: Transition from ectoderm convex cells to streak concave bottle cells in chicken (after http://www.gastrulation.org/).

Figure 1 - 4 :

 14 Figure 1-4: Transition from an ectoderm convex cell to a concave bottle cell in Drosophila (after http://www.gastrulation.org/).

  -2), applying a filter to remove the green pixels from the image, which do not bear any useful information for us. The resultant image (Figure 2-3 (Top)) is converted to 256 grey level image, to observe the intensities in the image (Figure 2-3 (Middle)) Applying simple contrast algorithms ([35]), helps us to have an exact visual idea of the number of cells in that region.

Figure 2 - 2 :

 22 Figure 2-2: Original image with its colour histogram. As a first step, simple thresholding technique was applied to get the binary image. The value of threshold can be easily chosen from the grey histogram of the image in Figure 2-3 (Middle). Due to noise in the image (green pixels in the original image in Figure 2-2), colour filtering might have removed some pixels on the boundary, hence the binary image might not have a continuous outline for each cell. An easy way to reproduce these pixels is through dilation process. A 3-by-3 mask is chosen to dilate the boundaries. The dilated image is shown in Figure 2-3 (Bottom).

Figure 2 - 3 :

 23 Figure 2-3: Segmentation process: colour filtering (top), gray scale image (middle), dilation and closure (bottom).

Figure 2 - 4 :

 24 Figure 2-4: Boundaries of cells with some intercellular objects (top). Intracellular objects removed (bottom). After dilation, cell boundaries could be drawn as shown in Figure 2-4 (Top). Objects not qualifying to be regarded as cell are removed by using number of pixels contained in each segment. Cell counting was done on the final image (Figure 2-4 (Bottom)). Similar process was employed for the rest of the images in the series. Some results are shown in Figure 2-5.

Figure 2 - 5 :

 25 Figure 2-5: Invagination process between initial and final stage from available data ([2-6]) just before invagination, rectangle showing the area of interest.

Figure 3 - 1 :

 31 Figure 3-1: Top: Cylindrical structure of Drosophila melanogaster embryo of length L and hemispheric left extremities. Middle: differentiation of the first bottle (or flask) cell. Bottom: equilibrium between external and internal forces.The end of the cell and tissue growth stops correspond to the end of the morphogenesis process: this can be observed when the second Thom's tissue law is applicable (the surface to volume ratio of an organ becoming adverse) or when the organ is completely covered by an anatomic boundary like an aponeurosis made of fibrous cells or an autoassemblage of extracellular ingredients. In both these cases, a couple of morphogens acting often simultaneously in opposite (e.g. a couple of activator and inhibitor like BMP-7 and BMP-2 in feather morphogenesis in the chicken[24][25][26][27] can induce the chemotactic motion of

Figure 4 - 1 :

 41 Figure 4-1: Left: cell contraction due to myosin (experiments in[36]). Centre: progressive invagination and streak formation. Right: myosin diffusion along a directrix of the cylindrical part of the embryo (green arrow) provoking the invagination before the tube closure (red star). The zero-diffusion domain for myosin is indicated in yellow.

Figure 4 - 2 :

 42 Figure 4-2: Top: part of the embryo made of cells considered as polyhedra with 8 vertices common with neighbouring cells and a centre (the centrosome)related to the vertices by actin filaments and micro-tubules (constituting the elastic and rigid cytoskeleton) on which myosin using ATP exerts contraction forces. Dark blue points are fixed points of the whole structure. Middle: profile of the structure under the contractile action of the myosin showing the start of the streak at the boundary of the hemispheric extremity of the embryo (red arrow). Bottom: experimental observation of the activity of myosin (fluorescent) with distal invagination The results of simulation of the hybrid model are given on Figures 4-1 and 4-2 showing the same phenomena as those observed in experiments: the progressive invagination starts at the extremities of the embryo and after propagates to the central cylindrical part of the embryo. The final step of the gastrulation needs proliferation in order to close the internal of the cell differentiation in bottle (or flask) cells contracted at their apical extremity, provoking during the first minutes of gastrulation a reorganization of the superficial cell layer of the embryo without division, leading to the formation of a streak. Resulting invagination starts experimentally at the two extremities of the embryo and propagates to

Figure 4 - 3 :

 43 Figure 4-3: Numerical Simulations showing the start of the streak at the boundary of the hemispheric extremity of the embryo its central part cylinder which will give birth to the intestinal tube of the adult animal. By following the progressive migration in embryo depth of the centrosomes it appears that the run is faster and deeper for curved cells at the extremities of the embryo than for central cells (Figure 4-4). The saturation curves representing this displacement behave like the curves representing the evolution of the cell diameter under progressively increasing forces applied externally to the cell, which is a way to induce cell contractions (possibly

Figure 4 - 4 :

 44 Figure 4-4: Displacement of centre particle (virtual centrosome) of a cell located on the cylindrical part of the embryo (blue), substantially smaller and slower than displacement of centrosome of cell located at the curved area (red).

Figure 4 - 5 : 5 CONCLUSION

 455 Figure 4-5: Simulated virtual spherical red blood cell (RBC) suspended in an hypotonic solution.Optical tweezers double trap is simulated by exerting locally a force Fs on two opposite nodes of the cell object contour (upper insert). The variation with load of the cell object diameter D(F) (in μm) in a plane perpendicular to the loading direction is simulated and compared to experimental data published by[31] (black squares). With appropriate scaling of the force, it is possible to adjust the elasticity modulus such that experimental mechanical response of RBC is nicely fitted in the linear elastic regime (red). Increasing the elasticity modulus induces a stiffer response which qualitatively reproduces the departure from the linear regime at larger traction forces (violet).

  of cones [Section II] and cells [Section III], our interest was to see whether we are able to track cells. Detecting and tracking cells constitute an important research area for biologists and scientists to analyze cell-growth or death in response to different diseases. Knowing the lineage of a sick cell could help us to

Figure 1 -

 1 2 depicts several stages of mitosis. Interphase is the "holding" stage or the stage between two successive cell divisions. In prophase, changes occur in both the cytoplasm and nucleus of the dividing cell. In Metaphase, chromosomes align at the metaphase plate of the dividing cell. In Anaphase, the paired chromosomes separate and move to opposite ends of the dividing cell. In Telophase, the nucleus of one cell is divided equally into two nuclei. During cytokinesis, a cleavage furrow is formed perpendicular to the mitotic spindle and gradually splits the cytoplasm and its contents into two daughter cells[7,8].

Figure 1 - 1 :

 11 Figure 1-1: Observing Mitosis With Fluorescence Microscopy; Mitosis In Rat Kangaroo Epithelial Kidney Cells [7]

Figure 2 - 1 :

 21 Figure 2-1: Block Diagram of Method

  separately and then superimposed back on the original image A typical tracking algorithm is shown by the block diagram given in Figure2-2.

Figure 2 - 2 :Figure 2 - 3

 2223 Figure 2-2: Cell Tracking Methodology

Figure 2 -

 2 2 is applied to first segment in the treated image in Figure 2-3. Red circle indicates the position and size of the segmented cell during first iteration. Green circle indicates the updated position and size of the same cell.

Figure 2 -

 2 Figure 2-4 shows a series of results of image treatment depicting cell tracking, also indicating a case where the cell divides in Figure 2-4(f).Figure 2-5 shows positions of

Figure 2 -Figure 2 - 4 Figure 2 - 5

 22425 Figure 2-4 Mitosis confirmed through image treatment[10]

2. 3 DISCUSSIONFrom Figure 2 - 5 ,

 325 we observe that the centroid of the segmented cell changes abruptly when a division occurs. The interdistance '𝑑' between the two consecutive centroids has been used as the indicator which is calculated as: 𝑑 = (𝑥 𝑖 -𝑥 𝑖-1 ) 2 + (𝑦 𝑖 -𝑦 𝑖-1 )2

Figure 2 -

 2 6 depicts a typical case. Most of these errors can be removed by employing a different segmentation technique, which might increase the preprocessing but could provide more consistent and reliable results. Another way to improve the algorithm is to track back the position of the previous segment.

Figure 2 - 6 : 2 . 3 . 4 . 5 .

 262345 Figure 2-6: False Tracking

Fig. 1 -

 1 Fig. 1 -The retina inside the eye structure

Fig. 2 -

 2 Fig. 2 -Contrast illusions: Hermann illusions (top); Mach bands illusion (bottom-left) and lateral inhibition with activation at short range (nearest neighbours neurons) and lateral inhibition at medium range (bottom-right)

Fig. 3 -

 3 Fig. 3 -Physiological and pathological retina. Top-left: Lateral inhibition due to horizontal-cell synapses. Top-right: Confocal microscopy slice of mouse retina having retinitis pigmentosa [34]. Bottom-left: Segmentation of cones and rods showing an important cell deficit in the Left Superior quadrant (LS). Bottom-right: Histogram of the intercept distances showing a significant increase of the inter-cell distance in the LS quadrant (dark blue).

Fig. 4 -

 4 Fig.4-Local organisation of cones and rods in macula (A) and retinal periphery (B)[38]

Fig. 5 -

 5 Fig. 5 -Translation invariant isotropic 4ݔ4 neural network with boundary neurons (violet), external field ݑ , synaptic weight ݑ ଵ ൌ ݓ , and non-linear effect coefficient ݑ ଶ ൌ ݒ

Fig

  Fig. -Top: Confocal slices of the retina: blue square indicates the slice of Fig. (Top right). Bottom: Segmentation of the rods and cones of this slice

Fig. 9 -

 9 Fig. 9 -Contrast enhancement. Bottom-left: Lateral inhibition causes the local contrast enhancement of the yellow square (with medium level in false colors) into a contrasted bright orange square. Top-right: Evolution in time of the DOG function representing an activation near the central neuron, ݅ (green links) and an inhibition (red links) farther away. Bottom: Same processing in grey level with initial image in the middle and contrasted on the right.

Fig. 10 -

 10 Fig. 10 -Contrasting and contouring medical images. Top-left: initial NMR image of a brain tumour. Top-middle: contrast enhancement with apparition of a central activity (blue arrows). Top-right: boundary of the compressed tissue (external snake-spline). Bottom-left: tumour segmentation. Bottom-right: tumour boundary (internal snake spline)

Figure 1 .

 1 Figure 1. Physiological and pathological retina. Top left: lateral inhibition due to horizontal cell synapses [after 3]. Top right: confocal microscopy slice of mouse retina with retinitis pigmentosa coming from T. Le ´veillard & J.A. Sahel [19]. Bottom left: segmentation of cones and rods with a cell deficit in the quadrant Left Superior (LS) [34]. Bottom right: histogram of the intercept distances showing an augmentation of the inter-cell distance in the quadrant Left Superior with respect to others Left Inferior (LI), Right Superior and Inferior (RS & RI) [34]. doi:10.1371/journal.pone.0006010.g001

Figure 2 .

 2 Figure 2. Contrast illusions. Top left: Hermann illusion with bright points at the intersection of grey stripes. Top right: Hermann illusion with grey squares at the intersection of white stripes. Bottom left: Mach bands illusion with enhancement of the vertical lines separating the different grey zones. Bottom right: lateral inhibition with activation at short range (nearest neighbour neurons) and inhibition at medium range (second Manhattan sphere) doi:10.1371/journal.pone.0006010.g002

Figure 3 .

 3 Figure 3. Contrast illusions. Top left: illusion of bright reinforcement at extremities and middle of white stripes in tangential vision. Top right: progressive change of a vertical bright stripe into bright spots (in false colours) during a morphogenetic process with lateral inhibition of morphogens [91]. Bottom left: bright and grey activities, respectively, near the center (vertical black line) and the extremities of the white horizontal diamonds. Bottom right: lateral inhibition simulated by a simple threshold neural network with h = w ii = 2, w ii21 = w ii+1 = 20.5 and a sequential updating from the left to the right hand side doi:10.1371/journal.pone.0006010.g003

Figure 4 .Figure 5 .

 45 Figure 4. Contrast illusions. Kanizsa pyramid: the lateral inhibition causes the sensation of seeing a 3-dimensional pyramid. doi:10.1371/journal.pone.0006010.g004

Figure 6 .Figure 7 .

 67 Figure 6. Contrasting and contouring medical images. Top left: initial NMR image of a brain tumour [48,49]. Top middle: contrast enhancement with occurrence of a central activity (blue arrows). Top right: boundary of the compressed tissue (using a snake spline). Bottom left: tumour segmentation. Bottom right: tumour boundary (internal snake spline) doi:10.1371/journal.pone.0006010.g006

Figure 8 .Figure 9 .

 89 Figure 8. Image attention processing. Top left: hippocampus-like neural network with lateral mixed action. Top-right: from A to F, progressive attentional focus by locally synchronizing the periodic signal associated to each pixel [56]. Bottom-right: desynchronization process between periodic activities of the neurons X i (i from 1 to n) doi:10.1371/journal.pone.0006010.g008

Figure 10 .

 10 Figure 10. Discrete segmentation and continuous contouring. Left: confinement tree [111]. Top-middle: level sets of the confinement tree in a brain tumour NMR slice [111]. Top-right: watershed tree [103]. Middle-left: level sets of the watershed tree [110]. Middle-right: watershed tree and landscape with different water levels. Bottom-left: on the left (respectively right) successful (respectively failed) contour of the remarkable Gaussian line in case of one (respectively two) isolated (respectively close) grey level peak(s) [41]. Bottom-right: 3D image of two close peaks [41] doi:10.1371/journal.pone.0006010.g010

Figure 11 .

 11 Figure 11. Continuous contrasting operators. Left: comparison between heat diffusion (a), Catte ´-Lions (b) [94] and non-isotropic (c) [93,95,96] contrasting operators. Right: non-isotropic contrasting operator and snakes contouring of the cardiac ventricular cavities with initial image (a), Canny-Deriche treatment [146] (b), non-isotropic processing (c), snakes splines contouring (d-f) doi:10.1371/journal.pone.0006010.g011

2

 2 ~M{m 2 , wherein M~ð T 2 y 1 ,y 2 ð Þdy 1 dy 2 :

  International Conference on Advanced Information Networking and Applications Workshops 978-0-7695-4019-1/10 $26.00 © 2010 IEEE DOI 10.1109/WAINA.2010.82

Figure 2 .Figure 3 .

 23 Figure 2. Top: the first shock of a spermatozoid head on the pellucid membrane of the ovule gives calcium waves (in green). Middle: tubulin waves in an artificial cell. Bottom: cylindric structure of an embryo

Figure 4 .Figure 5 .Figure 6 .

 456 Figure 4. Top: differentiation of the first bottle (or flask) cell. Middle: cells equilibrium between external and internal forces. Bottom: successive divisions of the first bottle cell leading to a progressive invagination of the furrow slice

Figure 9 )

 9 occurs for the first time a fluctuation in morphogen concentration passing over a certain threshold, then x i (0)=1 4) after a diffusion time equal to d, the first transition cell k (located at the interface between the cylindric and hemispheric parts of the embryo) receives an overconcentration of the morphogen, then x k (d)=0 and y k (d)=1. After a diffusion time d' (equal to d if the cell i is central), the second interface is reached by the diffusing morphogen at a cell k', then x k' (d')=0 and y k' (d')=1 5) k and k' have one intercell wall with less gap junctions (in red on Figure9), hence they can be transversally contracted after a time c>sup(d,d'), then: x k (c)=0 and y k (c)=1, x k' (c)=0 and y k' (c)=1 6) at time c+1, the left wall of the cells k and k' is contracted then the cell at the right (resp.left) hand size of k (resp. k') then they can contract and take the same state as k (resp. k') 7) at time p, the cells k a,d k' differentiated as bottle cells start to divide into 2 sister cells having the same state as their mother (for respecting the first Thom's cell law) 8) the mitosis wave is coming from the extremities of the streak toward the center and causes the invagination (see the slices of the Figure8). The automaton above is very caricatural but it renders into account the essential features of the ventral furrow formation.

Figure 9 .

 9 Figure 9. Results of an automaton simulation mimicking the progression from embryo extremities to center of the streak contraction and of the proliferation inside the ventral furrow

Figure 10 .

 10 Figure 10. Plant dysgenesis showing an absence of restitio ad (and a fortiori in) integrum during the correction of a gysgenesis due to a parasite (Top left: Cronartium ribicola) or to a branch lopping (Bottom right: Pinus radiata).The proliferation is supposed to be more important in the concave parts of the tree (the concavity is defined with respect to the sap canals) and less (resp. more) important in the zones of least diffusion (resp. accumulation) of the plant morphogen (auxin)

Figure 11 .

 11 Figure 11. Coexistence of two agonist (in blue) and antagonist (in red) morphogens with indication of the zones of least diffusion. When these zones coincide (top) for an precise viscosity ratio, the mean co-presence time is large, which can favor autoassemblages of functional structures (like transmembrane proteins inside a phospholipidic matrix) whose elements are induced by the morphogens

2. 1 .

 1 Image processing techniques: cell countouring and counting Many cell-based research studies require the counting of cells prior to beginning an experiment. Estimation of cell density in various regions of embryo is thus an integral part of such studies. Profile counts or stereological techniques could be used to have an estimate of the cell density in a particular region. The regions where cell density increases enormously 2011 Workshops of International Conference on Advanced Information Networking and Applications 978-0-7695-4338-3/11 $26.00 © 2011 IEEE DOI 10.1109/WAINA.2011.152

Figure1.Figure 2 . 2

 22 Figure1. Invagination process between initial and final stage from available data, [2-6] just before invagination, rectangle showing area of interest

Figure 3 :

 3 Figure 3: Left: cell contraction due to myosin (experiments [36]). Centre: progressive invagination and streak formation. Right: myosin diffusion along a directrix of the cylindric part of the embryo (green arrow) provoking the invagination before the tube closure (red star). The zero-diffusion domain for myosin is indicated in yellow.

Figure 4 :

 4 Figure 4: Top: part of the embryo made of cells considered as polyhedra with 8 vertices common with neighbouring cells and a center (the centrosome) related to the vertices by actin filaments and micro-tubules (constituting the elastic and rigid cytoskeleton) on which myosin using ATP exerts contraction forces. Dark blue points are fixed points of the whole structure. Middle: profile of the structure under the contractile action of the myosin showing the start of the streak at the boundary of the hemispheric extremity of the embryo (red arrow). Bottom: experimental observation of the activity of myosin (fluorescent) with distal invagination

Figure 5 :

 5 Figure 5: Displacement of centre particle (virtual centrosome) of a cell located on the cylindrical part of the embryo (blue), substantially smaller and slower than displacement of centrosome of cell located at the curved area (red).

Figure

  Figure6: Simulated virtual spherical red blood cell (RBC) suspended in an hypotonic solution. Optical tweezers double trap is simulated by exerting locally a force Fs on two opposite nodes of the cell object contour (upper insert). The variation with load of the cell object diameter D(F) (in µm) in a plane perpendicular to the loading direction is simulated and compared to experimental data published by[31] (black squares). With appropriate scaling of the force, it is possible to adjust the elasticity modulus such that experimental mechanical response of RBC is nicely fitted in the linear elastic regime (red). Increasing the elasticity modulus induces a stiffer response which qualitatively reproduces the departure from the linear regime at larger traction forces (violet).

LI LS RS RI E(X)

		45.4	47.4	36.95	32.15
	STDV(X)	26.6	28.8	25.00	18.47

Table 3 .

 3

	1: Characteristics corresponding to the horizontal intercept distances in the 4 quadrants LS,
	RI, RS and LI of Figure 3.14.

  STATE OF ART ……………………………………………………………………..........................52

		Table of Contents
	CHAPTER INTRODUCTION
	1.1.	SOME DEFINITIONS …………………………………………………………...………………….....50
	1.2.
	CHAPTER IMAGE PROCESSING

2.1. IMAGE PROCESSING TECHNIQUES: CELL CONTOURING AND COUNTING………………...…55 2.2. RESULTS OBTAINED FROM CELL CONTOURING AND IMAGE COUNTING...........................55 CHAPTER THE BIOMECHANICAL MODEL

PLoS ONE | www.plosone.org

June 2009 | Volume 4 | Issue 6 | e6010

ACKNOWLEDGEMENTS

I wish to express my gratitude to my supervisor, Mr. Jacques Demongeot, who has supported me throughout the course of my work. I thank him for introducing me to this

Acknowledgments

We are indebted to J. Mattes for helpful discussions and comments and to E. Greene for many suggestions and improvements.

ACKNOWLEDGMENTS

The data concerning the Pinus radiata tree deformation were obtained by F. Padilla (CONAF and CMM) and concerning the gastrulation by E. Wieschaus (cf. for example the sites http:www.ibioseminars.org and http:www.molbio1.princeton. edu/wieschaus/).

This work has been supported by the EC project Alfa IPECA and by the EC NoE VPH. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

SECTION II IMAGING AND MODELLING OF

Keywords: genetic regulatory networks, cell cycle control, microRNAs, graph frontier, circuits, gastrulation

The living organisms are very complex -part digital and part analogy mechanisms. J. von Neumann [1].

Modeling and Image