Cathy Danesin

Marina Ferreira

Pierre Degond

Eric Theveneau
email: eric.theveneau@univ-tlse3.fr

Anteroposterior elongation of the chicken anterior trunk neural tube is hindered by interaction with its surrounding tissues

Keywords: neural tube, chicken embryo, growth, tissue interaction, proliferation

Anteroposterior elongation of the chicken anterior trunk neural tube is hindered by interaction with its surrounding tissues

Cathy Danesin, Marina Ferreira, Pierre Degond, Eric

Introduction

In vertebrates, from gastrulation, the growth of the embryo is biased along the anteroposterior axis giving rise to an extensive elongation. The neural tube, which is the precursor of the central nervous system, sits at the center of the developing embryo and is no exception to this biased growth. In the chick embryo, between 2 days and 4 days of development, the posterior region of the neural tube (from somite 15 to the tip of the tail) grows by a factor 3.56 along the anteroposterior (AP) axis [START_REF] Ferreira | Interkinetic nuclear movements promote apical expansion in pseudostratified epithelia at the expense of apicobasal elongation[END_REF].

By contrast, the neural tube only grows by a factor 2.45 and 2.3 in the dorsoventral (DV) and apicobasal (AB) axes, respectively [START_REF] Ferreira | Interkinetic nuclear movements promote apical expansion in pseudostratified epithelia at the expense of apicobasal elongation[END_REF]. Several phenomena could contribute to this favored anteroposterior growth. Caudally, there is the tail bud, a progenitor zone that adds cells to the growing neural tube and mesoderm in a rostrocaudal manner [START_REF] Sambasivan | Neuromesodermal Progenitors: A Basis for Robust Axial Patterning in Development and Evolution[END_REF]. In addition, it has been proposed that the nascent mesoderm and the underlying notochord would exert a local pressure on the forming neural tube favoring its elongation along the anteroposterior axis [START_REF] Xiong | Mechanical Coupling Coordinates the Co-elongation of Axial and Paraxial Tissues in Avian Embryos[END_REF]. Further, analyses of cell division at early stages of chicken neural tube development (up to 2 days, HH12) suggested that neural cells preferentially divide along the anteroposterior axis [START_REF] Sausedo | Role of nonrandomly oriented cell division in shaping and bending of the neural plate[END_REF]. Furthermore, neural progenitors display a differential proliferation rate along the anteroposterior axis. Indeed, as the neural tube develops, proliferation slows down and differentiation kicks in such that anterior portions of the neural tube have proportionally fewer progenitor cells with an increased mean duration of their cell cycle [START_REF] Molina | Playing with the cell cycle to build the spinal cord[END_REF].

Altogether, these observations suggest that the anteroposterior bias of the neural tube growth is primarily due to fast proliferation with oriented divisions and local pressure from mesoderm in the posterior region which would concur to promote anteroposterior elongation. While this model of growth is widely accepted, numerous questions remain open. Is AP-biased growth a feature of the neural tissue and thus also observed outside of the most posterior region? As mesoderm develops into somites, is the neural tube still under lateromedial pressure favoring AP elongation? Are AP oriented divisions restricted to early stages or a feature of the neural plate/tube regardless of the stage? What is the actual amplitude of the lengthening of the mean cell cycle duration in a given region over time and could it bias growth in a given axis or account for changes in growth rates?

To address some of these points, we used a combination of classical embryology techniques coupled with mathematical modelling of the proliferating epithelium. First, we found that growth of the neural tube is still anisotropic outside of the tail bud region. However, by contrast with the posterior region, AP elongation is slower than growth in the other axes between 2 and 4 days of incubation. Second, extracting a portion of the anterior trunk from the embryo or detaching the neural tube from its surrounding tissues leads to an AP expansion of the neural tube. This indicates that AP growth of the anterior neural tube is hindered by neighboring tissues and its adhesion to the matrix. Third, we found no evidence of a directional bias of cell division along AP axis in the developing neural tube between 2 and 4 days of development. Fourth, we calculated the expected mean cell cycle duration to account for the growth of the neural tube facing the forelimb region from 2 to 4 days of development. We used a mathematical model to simulate different regimes of proliferation using cell cycle durations extracted from our experimental data or from the literature. We found that the growth of the neural tube at forelimb level between 2 and 4 days of development corresponds to a shift of mean cell cycle length from 8h50 to 13h30min in neural progenitor cells together with a proportionally decreasing progenitor pool.

This shifts is likely due to a global change of mean cell cycle length rather than to the emergence of lengthy cycles that would skew the mean towards higher values. Interestingly, the scale of this change is of the same order of magnitude as the observed change in tissue growth suggesting that the lengthening of the cell cycle is one of the main driver of the slowdown of neural tube elongation between 2 and 4 days.

Altogether, our data indicate that neural tube growth goes from a favored AP elongation in caudal-most regions, most likely due to addition of new cells from the activity of the tail bud and lateral pressure from the mesoderm, to a hindered AP growth due interaction with the matrix. Concomitantly, overall growth slows down due to changes in cell cycle duration.

Results

1. Anterior trunk neural tube growth is slower along the anteroposterior axis than in the other axes Does neural tube growth switch from being AP-biased in posterior regions to isotropic in anterior regions? To assess this, we monitored the growth of the neural tube anterior to the tail bud. To measure AP growth, we injected the fluorescent tracer DiI at different locations along the neural tube (Fig. 1A) in embryos at 48h (Fig. 1B) or 72h of incubation (Fig. 1C). From these, we measured the distance between the different DiI points in ovo at t0 and 24h later and calculated the mean growth rate along the AP axis in each time window (Fig. 1D-E). In parallel, we assessed neural tube growth in the DV and AB directions by measuring these axes on transversal sections of embryos at identical region and stages (Fig. 1F-J).

Interestingly, the mean AP growth rates of 1.44 between 48 to 72h and of 1.27 between 72 and 96h are significantly lower than the observed growth rates for the DV and AB axes at both time points (Fig. 1K).

These results indicate that neural tube growth is not isotropic outside of the tail bud region. Instead, the anterior neural tube displays a reversed anisotropic growth in which AP growth is slower than growth in the other axes.

Anterior trunk neural tube anteroposterior elongation is hindered by interaction with surrounding tissues

In posterior regions, local pressure from surrounding tissues, including the nascent mesoderm and the underlying notochord, has been proposed to favor AP growth [START_REF] Xiong | Mechanical Coupling Coordinates the Co-elongation of Axial and Paraxial Tissues in Avian Embryos[END_REF]. Thus, we asked whether interaction with the surrounding tissues might have a similar role in the neural tube facing somites or on the contrary might contribute to the observed slower AP growth. First, we measured the length of the neural tube of individual embryos along the 5 most posterior somites in ovo and after dissection and fixation (PFA4%) of that same portion at 2 days of development (Fig. 2A). The mean AP length of the neural tube was significantly higher ex ovo compared to in ovo with a 14% length increase upon extraction (Fig. 2B). Noticeably, the mean somite length also significantly increased after dissection but only by 6.9% in the same explants (Fig. 2C), which was significantly lower than the neural tube expansion (see figure legends for statistics). Given that aldehyde fixatives lead to tissue shrinkage, these values are likely to be slightly underestimated.

Next, we assessed the role of the interaction of the neural tube with the surrounding tissues. For this, we again harvested portions of the trunk corresponding to the 5 most posterior somites at 2 days of development that we either fixed immediately after dissection or after 2h in culture medium (DMEM) with or without Dispase, which digests collagens and Fibronectin and thus cuts links between the neural tube and underlying basal lamina (Fig. 2D-E). After 2h in culture, in explants in DMEM, the neural tube kept its overall shape and length compared to explants fixed at t0 (Fig. 2D-E). By contrast, the Dispase treatment led to drastic changes of the explants shape, visible to the naked eye, with the superficial ectoderm detaching and rolling up. In these Dispase-treated samples, the neural tube buckled and its AP length significantly increased by 19% compared to control DMEM explants (Fig. 2D-E). Therefore, adhesion of the anterior neural tube to the extracellular matrix plays a role in restricting its AP growth.

Interfering with the cell-matrix interaction also prevents normal interaction with the surrounding tissues. To assess whether the elongation observed in Dispase-treated samples was due to the loss of physical interaction with the surrounding tissues and/or to the loss of adhesion to the matrix per se, we performed a similar assay but this time, the surrounding tissues were removed by microdissection.

Importantly, the proximal parts of the somites were kept intact to ensure that the basal lamina of the neural tube was not damaged by the procedure (Fig. 2F). In this case, the dissected neural tube elongated by only 5% compared to explants with surrounding tissues, a difference that was not found to be statistically significant, indicating that the neural tube AP extension obtained after Dispase treatment was mostly due to the loss of interaction with the basal lamina (Fig. 2F-G).

Collectively, these data indicate two levels of physical constraints that hinder AP elongation in the anterior neural tube. One associated with being part of the whole embryo in ovo that is released upon extraction of a whole piece of trunk and a second related to the anchorage to the matrix that is released upon enzymatic degradation of collagens and fibronectin.

3. Anchorage of the anterior trunk neural tube to the extracellular matrix favors apicobasal elongation Next, in order to assess the consequences of the loss of cell-matrix interactions on the other axes, we monitored the size of the neural tube in DV and AB directions and its pseudostratification on transversal sections in DMEM and Dispase-treated samples. For that, we stained DNA using DAPI and Actin with Phalloidin (Fig. 3A-B). The DV height, AB width and number of nuclear pseudolayers were similar in controls (t0) and DMEM incubated samples (Fig. 3C-E). By contrast, in Dispase-treated samples the DV height increased by 13.6% (Fig. 3C) while the AB width (Fig. 3D) and pseudostratification (Fig. 3E) decreased by 44 and 38% respectively. This shows that anchorage to the matrix is essential for normal AB elongation of the cells and the associated pseudostratification. In addition, matrix degradation led simultaneously to AP lengthening, DV expansion and AB shortening of the neural tube. Importantly, these changes take place in a short time window (2h), during which one can reasonably assume that cell volumes remain constant. Therefore, this result suggests that AB elongation and pseudostratification occur at the expense of AP and DV growth of the neural tube, in line with our previous observations [START_REF] Ferreira | Interkinetic nuclear movements promote apical expansion in pseudostratified epithelia at the expense of apicobasal elongation[END_REF].

4. Mitoses orientations in the neural tube are not biased along a specific direction between 2 and 4 days of development Next, we wondered whether changes in cell proliferation could account for tissue growth dynamics either due to a change of orientation of cell divisions or to a change of the pace of proliferation.

First, we assessed whether cell divisions displayed a preferred orientation in the neural tube as its AP growth gets hindered. Previous analyses on fixed samples were performed using DNA labelling. It was proposed that in the trunk of embryos up to stage HH12 (2 days of development), nearly 55% of all divisions in the neural tube were oriented along the AP axis [START_REF] Sausedo | Role of nonrandomly oriented cell division in shaping and bending of the neural plate[END_REF]. To analyze the orientation of cell division, we performed an immunostaining for γ-tubulin together with DNA labelling using DAPI, mounted the samples in open book configuration (Fig. 4A) and measured the angles of divisions according to the AP axis of the neural tube (Fig. 4B). As it is known that metaphasis plates rotate until the last minute before division, we restricted our analysis to situations with a clear separation of the genetic material to two sides of the cell (anaphase) to avoid the possible misleading interpretation of early mitotic figures. We used the axis between the virtual line passing through the two γ-tubulin dots and the AP axis of the neural tube to measure the angle of division . With this strategy, we analyzed the orientation of anaphases in the neural tube in 48h embryo at the level of the presomitic mesoderm (psm) or the forelimb. In these two regions, and contrary to what had been previously proposed, the divisions were evenly distributed along the AP, DV and intermediate/diagonal orientation (Fig. 4E-F; 48h and 48h psm). Next, we analyzed the orientation of cell divisions at the level of the forelimb along further development, at 72 and 96h, and found again evenly distributed orientations of cell division along AP and DV axes at both stages (Fig. 4E-F; 72h and 96h). Therefore, these results indicate that the orientation of cell divisions is unlikely to contribute to the switch from AP-biased to AP-hindered growth of the neural tube from posterior to anterior regions or to changes in growth in the anterior neural tube over time.

Dynamics of the anterior trunk neural tube growth over two days at forelimb level corresponds to a moderate and homogenous lengthening of the mean cell cycle duration

In absence of oriented cell divisions, we wondered whether a change of the pace of proliferation might account for the observed slowdown of growth between 2 to 4 days of development. Indeed, our measurements indicated an important change of anterior neural tube 3D growth between 72h and 96h as compared to 48-72h growth rate (Fig. 1K). Between 2 to 3 days the mean 3D growth rate is 4.2 (1.44 x 1.7

x 1.72) while between 3 to 4 days it is only 2.74 (1.27 x 1.52 x 1.42) which correspond to a 35% drop. This raises the question of whether cell proliferation dynamics in anterior trunk neural tube could contribute to tissue growth temporal variations. Mean cell cycle length has been reported to vary widely in time and space during chicken embryo development and in particular in the neural tube (Fig. 5A). There is a general consensus that the duration of the cell cycle increases as the neural tube develops but in the chicken embryo, there has been no global study in a given region over time to confirm this or to estimate the amplitude of the change [START_REF] Molina | Playing with the cell cycle to build the spinal cord[END_REF] or whether such changes match observed changes in growth.

To estimate cell cycle durations in the anterior chick neural tube between 48h and 96h of development, we quantified the progenitor and total cell populations. We used Sox2 immunostaining to label progenitor cells and DAPI to label all cells at 48h, 72h and 96h of development at the forelimb level, on transversal sections (Fig. 5B-C). Since we found no oriented cell divisions at these stages of development, between two time points, observation of cell numbers on series of transversal sections is a proxy of cell numbers in a 3D volume. Indeed, without a bias of cell division orientation, the proportion of daughter cells that are excluded from a given transversal section is the same as the proportion of newly incorporated daughter cells coming from divisions that took place in adjacent transversal sections. Thus, the evolution of cell number in sections over time reflects the evolution of the whole population.

Transversal sections at all time points were generated with the same thickness (14 µm). However, because of tissue growth, a 14 µm-thick section of the neural tube at 48h of development will no longer be 14 µm-thick at 72 or 96h. According to the AP elongation rates that we calculated from our DiI injections (Fig. 1B-E), a 14 µm-thick section will develop into a 20.16 µm section (14 x 1.44) at 72h and into a 25.6 µm-thick section at 96h (14 x 1.44 x 1.27). Therefore, the number of cells observed on 14µm sections at 72 and 96h needs to be corrected by the observed AP growth rates, respectively x1.44 at 72h and x1.83 (1.44 x 1.27) at 96h. It is also important to note that cell translocation can occur independently of cell division, as previously reported in the cephalic neural tube in mice [START_REF] Mcshane | Cellular basis of neuroepithelial bending during mouse spinal neural tube closure[END_REF].

Therefore, extrapolating the dynamics of the whole population from transversal sections will only be correct if such translocations events are very rare or evenly distributed between all axes.

We first calculated the number of progenitor cells and total cell number dynamics from histological sections of the neural tube (Fig 5C). From these data, we calculated the expected doubling time required for the progenitor region at a given time point to generate the whole population at the next time point (see methods for the calculations) with the two aforementioned assumptions that there are no oriented divisions (our data) and that there no significant neural cell translocations (the simplest hypothesis). From this, we calculated a mean cell cycle duration of 8h50 for neural progenitors between 48 and 72h of development and of 13h30 between 72 and 96h of development (Fig. 5A, green bars).

To compare the relevance of these cell cycle values with the previously published values, we ran a computational model that we formerly validated [START_REF] Ferreira | Interkinetic nuclear movements promote apical expansion in pseudostratified epithelia at the expense of apicobasal elongation[END_REF] to simulate the growth of the neural tube cell population over two days (Fig. 5D-E). Importantly, and unless otherwise specified, all mean cycle values used for simulations are allowed a narrow range of 1h above and below the mean. We first assumed that all cells are proliferative and simulated growth regimes with the fastest (5h, Fig. 5D, open circles) and the slowest cell cycle published values (28h, Fig. 5D, black dots) [START_REF] Fujita | Kinetics of cellular proliferation[END_REF][START_REF] Wilcock | Mitotic spindle orientation distinguishes stem cell and terminal modes of neuron production in the early spinal cord[END_REF]) and an intermediate regime in which neural progenitor cells can proliferate with a cell cycle of any duration between 5 and 28h (Fig. 5D, grey dots). All data were normalized to t0 such that each population virtually starts with one cell, allowing direct comparison of growth dynamics (Fig. 5D). This shows that none of the proliferation regimes tested accounts for the actual growth of the population.

Cycles of 5h give rise to an excessive growth whereas cycles of 28h, previously observed in vivo using live imaging [START_REF] Wilcock | Mitotic spindle orientation distinguishes stem cell and terminal modes of neuron production in the early spinal cord[END_REF], or an intermediate 5-28h cell cycle length do not allow production of enough cells.

We then use our calculated values of cell cycle length to perform new simulations of neural cell growth between 48 and 96h of development, with all cells proliferating, normalized all datasets to t0, and compared them to the biological data (Fig. 5E). We modeled 3 conditions, constant cell cycle length of 8h50 (Fig. 5E, open circles), constant cell cycle length of 13h30 (Fig. 5E, black dots) or a mixed regime with a cell cycle of 8h50 for the first 24h and then switching to 13h30 for the next 24h (Fig. 5E, grey dots). The mixed regime was the closest to the biological data (Fig. 5E, green dots). However, the mixed regime still produced too many cells compared to the observed cell number over 2 days (+21%), with a divergence between the predicted growth rate and the observed cell numbers that is mostly observed in the last 24h (Fig. 5E). This discrepancy could be due to an underestimation of the actual cell cycle length between 3 and 4 days of development or to the fact that we modelled a population with 100 % of progenitors while the biological data indicate that the proportion of proliferative cells (Sox2+) progressively shrinks over time (Fig. 5B-C).

Thus, we plotted the percentages of differentiated cells at 48, 72 and 96h of development and used a LOWESS fit to extrapolate the proportion of differentiated cells at intermediate stages of development (Fig. 6A). From this, we calculated the proportion of available progenitor cells over time (Fig. 6B). We then used this information to repeat simulations of the mixed 8h50-13h30 proliferation regime with different proportions of proliferative cells (Fig. 6C): 100% throughout (Fig. 6C, light grey dots), 100% for 24h and then 85% (Fig. 6C, dark grey dots), 100% for 24h, 85% for 12h and 70% for another 12h (Fig. 6C, black dots). The first two conditions give rise to an excess of 21% and 8.5% of cells compared to the biological data after two days, respectively. The last condition is the closest to the biological data with a difference in total cell number inferior to 3% with the biological data. These results show that a slight lengthening of the cell cycle duration (+50%) coupled with a progressively declining proportion of progenitor cells is sufficient to reproduce the dynamics of the neural tube cell population at forelimb level between 2 and 4 days of development.

The long cell cycle value that we calculated between 3 and 4 days of development is significantly shorter than some values that have been previously proposed in the literature (up to 28h). The mean increase of cell cycle duration between 2 and 4 days of development could either occur because all cell divisions take slightly longer or by an increase of the range of cell cycle durations. In the last case, some cells would still cycle fast while others would take a very long time skewing the mean towards higher values while accounting for observed lengthy cycles. We tested these two scenarios by repeating simulations with the declining progenitor populations in which the slowdown of cell cycle takes place with the aforementioned narrow range (+/-1h) around the mean (Fig. 6D, black dots) or a wide range (+/-5h) around the mean (Fig. 6D, open circles). In this last scenario, cells can proliferate with any cell cycle duration between 8h30 and 18h30. Despite the mean duration being the same at 13h30, the wide range simulation generates 16% extra cells compared to the biological data. This indicates that the shift of cell cycle duration is most likely driven by an overall lengthening of cell cycle in all proliferative progenitors rather than a widening of the range of cell cycle values towards lengthy cycles.

Can these moderate changes in cell cycle pace account for the observed slowdown of the neural tube growth over two days? We compared the change in growth rates in 3D, calculated from the growth rate in each axis, with the change in growth rates in terms of cell number (which is driven by the lengthening of the cell cycle) during the same period. We have a mean 3D growth rate of the neural tube of 4.21 between 2 and 3 days of development and of 2.74 between 3 and 4 days (Fig. 1K). This gives a ratio of 0.65 for 3D growth between these two periods, corresponding to a 35% drop. By comparison the cell cycle goes from 8h50 (8.8h) between 2 and 3 days to 13h30min (13.5h) between 3 to 4 days, which also gives a ratio of 0.65 of cell cycle variation. This indicates that the slowdown of the cell cycle and that of the 3D growth are of similar amplitude. This strong correlation between cell cycle variation and tissue growth suggests that the lengthening of the cell cycle is likely to be one of the main driver of the declining rate of 3D growth in the anterior neural tube. However, one should keep in mind that other changes that are taking place concomitantly could contribute to the slowdown of tissue growth such as the declining number progenitors and the delamination and subsequent differentiation of neurons.

In conclusion, our data indicate that the switch from AP-favored to AP-hindered growth of the neural tube primarily depends on interaction with the matrix whereas the subsequent slowdown of growth observed over two days is likely to be driven by a moderate and homogenous lengthening of the duration of the cell cycle.

Discussion

Anteroposterior growth is notoriously biased along the AP axis and this is mainly driven by the progenitor zone which adds cells to the growing tissues. Surprisingly, we found that rather than switching to isotropic growth outside of the tail bud region, the neural tube displays a reversed anisotropic growth with a slower AP elongation relative to AB and DV axes. More precisely, AP growth is hindered by anchorage to the matrix. Surprisingly, we found that simply extracting a portion of the trunk was sufficient to promote an immediate AP elongation of the neural tube, suggesting that it is under slight AP compression in ovo. In addition, impairing interactions with the matrix also led to an additional increase in AP length and DV height. This occurred concomitantly to a reduction of the AB width of the cells. This echoes our previous observation that neuroepithelial cells in the forelimb region see their volume shrink between 2 and 4 days of development while their size along the AB axis more than doubles in the same time frame [START_REF] Ferreira | Interkinetic nuclear movements promote apical expansion in pseudostratified epithelia at the expense of apicobasal elongation[END_REF]. This dramatic change of cell shape while cells are getting smaller means that most of the cell material is used along the AB axis. Thus, any mechanism that favors or impairs AB elongation would impact growth in the other axes. One example is interkinetic nuclear movements (INM) (see below).

In posterior regions, it has been proposed that the nascent mesoderm and the underlying notochord contribute to a caudally oriented force on the growing neural tube favoring AP elongation [START_REF] Xiong | Mechanical Coupling Coordinates the Co-elongation of Axial and Paraxial Tissues in Avian Embryos[END_REF]. Our dissection and Dispase experiments suggest that this relationship no longer stands in anterior regions where the neural tube faces somites. One possibility could be that somitogenesis itself contributes to this hindered AP elongation. Somites aggregate by epithelialization from a loose mesenchyme called the presomitic mesoderm. One could think that the epithelialization and the local reorganization of the matrix that occurs during the process [START_REF] Martins | Dynamic 3D cell rearrangements guided by a fibronectin matrix underlie somitogenesis[END_REF] might lead to a local compressive force along the AP axis when somite precursors cluster to epithelialize. Indeed, the region of the anterior presomitic mesoderm that contains the next somite to be formed is approximately 90 µm long whereas a newly epithelialized somite is 50-60 µm long [START_REF] Martins | Dynamic 3D cell rearrangements guided by a fibronectin matrix underlie somitogenesis[END_REF]. This local "shrinkage" or "condensation" of the paraxial mesoderm during somitogenesis might contribute to a local compressive force aligned with the AP axis that could participate to the hindered AP elongation of the neural tube at this level. Interestingly, we also found that somites elongate in AP direction after dissection, but to a significantly lower extent than the neural tube. A possible explanation for this is that the neural tube is a long epithelial structure with bilateral anchoring to the matrix, while paraxial mesoderm consists of a succession of epithelial somites separated by matrix boundaries that may buffer compression forces.

Regarding the cell cycle, we found, as previously reported, that the mean cell cycle duration increases over time. However, our calculations and simulations indicate that the actual amplitude of the proliferation slowdown at the level of the forelimb is about 35% over 2 days. Interestingly, these results are in better agreement with analyzes in mice where cell cycle length was monitored over time in one region [START_REF] Kicheva | Coordination of progenitor specification and growth in mouse and chick spinal cord[END_REF] than with studies in chicken that compared different regions along the AP axis [START_REF] Olivera-Martinez | Major transcriptome re-organisation and abrupt changes in signalling, cell cycle and chromatin regulation at neural differentiation in vivo[END_REF]. This suggests that there may be important regional differences in cell cycle length. Of note, our estimated cell cycle values are line with what was measured using none invasive methods such as EdU/BrdU incorporations but differ significantly from cell cycle values that were measured by live imaging. While some interpreted this difference as a new level of complexity discovered by live imaging, one could also argue that these long cell cycle values are actually induced by the protocol of imaging itself. For live imaging, neural tubes have to be cooled down, sliced, cultured under a microscope and intensively exposed to laser along z stacks over long periods of time. All these treatments might impact on cell fitness and the pace of the cell cycle (e.g. changes in temperature, phototoxicity). Some cells may not be affected but some may actually be struggling to divide normally leading to artefactual lengthy cell cycle durations. Strikingly, we found that the reduction of cell proliferation between 2 and 4 days of development, is of the same order of magnitude as the observed slowdown of 3D growth during the same time period indicating that the slowdown of proliferation is likely to be the main driver of this change. Nevertheless, the relationship between cell proliferation and tissue growth is complex. As we found no preferred orientation of cell division, the simplest assumption would be that proliferation impacts growth of the tissue equally in all directions. However, while this might be true for a ball of cells, it is less likely for a pseudostratified epithelium. Indeed, pseudostratified tissues are monolayers and all new cells are added within the plane of the epithelium. Thus, proliferation should solely impact apical expansion, which in the case of the neuroepithelium corresponds to growth along the AP and DV axes. In addition, proliferation in neuroepithelial cells involves interkinetic nuclear movements (INM) with cells sending their nuclei apically prior to mitosis. We previously showed that INM favor expansion of the apical domain [START_REF] Ferreira | Interkinetic nuclear movements promote apical expansion in pseudostratified epithelia at the expense of apicobasal elongation[END_REF] which should further regulate AP and DV expansion.

However, this is mitigated by cell shape changes over time, neuroepithelial cells growing extensively along their apicobasal axis while experiencing a reduction of their volume [START_REF] Ferreira | Interkinetic nuclear movements promote apical expansion in pseudostratified epithelia at the expense of apicobasal elongation[END_REF]. Therefore, while new cells are added in the plane of the epithelium a significant part of the cell material is redistributed along the AB axis, thus competing with INM-driven apical expansion.

We previously discussed that changes in the size of the apical domain, either by contractility or by INM, could impact pseudostratification and AB elongation [START_REF] Ferreira | Interkinetic nuclear movements promote apical expansion in pseudostratified epithelia at the expense of apicobasal elongation[END_REF]. The logic being that enlarging the apical domain would allow cells to redistribute their nuclei within the 2D plane of the epithelium while a reduction of the apical domain would force pilling of the nuclei along the AB axis.

However, this relationship is most likely to apply only in cells with a short AB axis as very elongated cells will be less influenced by small variations of their apical domain. Therefore, in the case of the anterior neural tube between 2 and 4 days of development, the increase of the mean cell cycle duration would primarily contribute to the pace of AP and DV expansion rather than controlling the balance between apical expansion and AB elongation. Changes in cell cycle dynamics have been linked to fate decisions in the neuroepithelium [START_REF] Molina | Playing with the cell cycle to build the spinal cord[END_REF]. Given the link between INM and cell and tissue shape, it is also possible that changes in the dynamics of the cell cycle directly contribute to the morphogenesis of the neural tube itself. This way, proliferation would coordinate identity and shape of neuroepithelial cells ensuring the formation of a functional organ.

Materials and Methods

DiI injections

After 48h of incubation, eggs were windowed using curved scissors. A solution of CM-DiI (Molecular Probes C7000) diluted 1 : 10 from a 0.5% stock solution in Ethanol 100%, in 0.3M sucrose solution as previously described [START_REF] Kulesa | In ovo time-lapse analysis of chick hindbrain neural crest cell migration shows cell interactions during migration to the branchial arches[END_REF] was then injected using a pulled glass needle.

Two to three injections were performed along the AP axis. A picture was taken before closing the eggs and putting them back in the incubator for 24h. The center of each DiI spot was used for measurements of distance.

Whole-Mount Immunostaining

After fixation in 4% paraformaldehyde, embryos are washed in PBS1X / 0.1% Tween 20 and incubated 60 minutes in PBS1X/1% Triton/5% Bovine Serum Albumin (BSA). Embryos are then incubated with primary antibody in PBS1X/1%BSA/0.1% triton 24-48 hours at 4°C (anti-Sox2, Merck Millipore AB 5603; anti-γ-tubulin, Sigma GTU-88), washed in PBS1X/0.1% triton X100 and incubated with secondary antibodies. All secondary antibodies were coupled with Alexa fluorophores and purchased from Invitrogen Molecular probes (Life Technology) and incubated with DAPI or Phalloidin in PBS1X/0.1% triton X100 overnight at 4°C. Embryos are washed in PBS1X/0.1% triton X100 before mounting in mowiol and observation or embedding for histology.

Cryosections and immunostaining

After fixation in 4% paraformaldehyde (PFA 4%) or immunostaining, embryos are washed in Phosphate Buffer (PB) and embedded in PB/15% sucrose/7.5% gelatin as previously described [START_REF] Andrieu | MMP14 is required for delamination of chick neural crest cells independently of its catalytic activity[END_REF]. Individual blocks are cut under a dissecting microscope to position the embryo in the desired orientation for sectioning. Blocks are frozen in isopentane (Sigma, 615838) at -70°C and stored at -70°C until sectioning. Cryosections were performed using a cryostat Leica CM1950 as previously described [START_REF] Theveneau | Ets-1 confers cranial features on neural crest delamination[END_REF]. Sections of 14µm were incubated in PBS1X at 42°C until gelatin was totally removed, treated with PBS1X/1% triton X100/5% BSA for 1 hour for permeabilization and blocking.

Primary antibodies were diluted as mentioned above and disposed on sections overnight at 4°C under a coverslip. Secondary antibodies were diluted as mentioned previously and applied on sections for 2 hours at room temperature. All washes were done in PBS1X.

Measurement of neural tube and somite length

Neural tube and somite lengths were measured using the line tools in FIJI. For all explants, images were taken dorsally except for DMEM explants (Fig. 2D) which tend to slightly curl. Thus, dorsal views lead to an underestimation of the actual anteroposterior length of the tissue. To circumvent this problem, neural tubes from DMEM explants were separated and imaged laterally.

Quantification of cells

Mean area occupied by a DAPI stained nuclei was estimated in FIJI by manually segmenting at least 20 nuclei either in Sox2+ or in Sox2-areas, on multiple cryosections at 48, 72 and 96h of development. At each stage, Sox2+ and Sox2-areas were measured manually on each cryosection and estimation of the number of nuclei was obtained by dividing this surface by the corresponding average nucleus surface.

Estimation of mean cell cycle durations

At 48h, a 14 µm-thick transversal section of the neural tube contains 318.6 cells on average whereas at 72h there are 1208 cells. Since we found no preferred orientation of cell division, daughter cells that may be produced outside of the transversal section would be replaced by equivalent daughter coming from divisions in adjacent sections. In addition, the relationship between the evolution of the population in transversal sections and the whole tissue only stands if there are no spatially biased postmitotic intercalation of cells. We postulated that such intercalation are negligible. Therefore, using these two assumptions, we can extrapolate the actual number of cells at 72h by multiplying the number of cells on a 14 µm-thick section by the AP elongation factor between 2 and 3 days of development (1.44). This gives a mean cell number at 72h of 1739. From these cell numbers, we can calculate the growth rate by dividing the total number of cells at 72h (1739) by the total number of cells at 48h (318.6), in this case 5.46. Then, from this, the expected doubling time is : 24h (the interval between 48 and 72h of development), times 2, divided by the actual growth rate between the two considered time points, here (24h x 2) / 5.46 = 8.8 hours of doubling time (8h50). We can apply the same calculation from the next interval between 72 and 96h. Although this time, the total number of cells generated between 72 and 96h corresponds to the total number of cells at 96h minus the number of differentiated cells (Sox2-negative cells) already present at 72h corrected by the AP elongation factor between 3 and 4 days of development (1.27). This gives a total of 3645.7 cells at 96h that were generated by the mean 1026 Sox2+ neural progenitors at 72h. From this, we get a growth rate of 3.55 which gives a doubling time of (24 x 2) / 3.55 = 13.5 hours (13h30).

Mathematical Modelling

The complete description of the mathematical model is available in [START_REF] Ferreira | Interkinetic nuclear movements promote apical expansion in pseudostratified epithelia at the expense of apicobasal elongation[END_REF] In the present study, we only used the cell cycle clock to model cell population growth. Each simulation was done 20 times and mean values are plotted. Each initial condition was modelled for 12h. The mean growth rate was computed and applied to calculate the evolution of the cell population over the next 48h.

Image acquisition and analysis

Images were taken on a stereomicroscope Leica MZF10F equipped with a camera Leica DFC450C and the LAS software, a Nikon Eclipse 80i equipped with a camera DXM1200C and the NIS-elements software, a Zeiss Axio Imager 2 equipped with a Hamamatsu ORCA Flash 4 camera and the Zen2 software, B, transversal sections of the chicken neural tube at the level of the forelimb region at 48, 72 and 96h of development stained with DAPI (DNA, green) and Sox2 (neural progenitor cells, magenta). C, mean number of total neural tube cells (black line) and Sox2-positive progenitors (grey line) per time point, per transversal section (n 48h = 3 embryos-28 sections, n 72h = 3 embryos-28 sections, n 96h = 4 embryos-40 sections), values for 72 and 96h were corrected by the corresponding AP growth rates (see methods). D, simulation of neural tube population growth versus biological data (green crosses) for 48h normalized to t0 using fastest (5h, open circles) and slowest (28h, black dots) published durations of cell cycle in the trunk neural tube and an intermediate regime (grey circles). E, simulation of neural tube population growth versus biological data (green crosses) for 48h normalized to t0 using mean cell cycle durations estimated from our biological data, i.e. constant 8h50 (open circles), constant 13h30 (black dots) and a mixed regime of 8h50 for the first day and 13h30 the second day (grey dots), considering that all cells proliferate throughout the simulation. Scale bars in B, at 48h, 50 µm ; at 72 and 96h, 80 µm. References in panel A: [START_REF] Fujita | Kinetics of cellular proliferation[END_REF][START_REF] Kicheva | Coordination of progenitor specification and growth in mouse and chick spinal cord[END_REF][START_REF] Langman | Behavior of neuro-epithelial cells during closure of the neural tube[END_REF][START_REF] Dreau | The strength of SMAD1/5 activity determines the mode of stem cell division in the developing spinal cord[END_REF][START_REF] Olivera-Martinez | Major transcriptome re-organisation and abrupt changes in signalling, cell cycle and chromatin regulation at neural differentiation in vivo[END_REF][START_REF] Peco | The CDC25B phosphatase shortens the G2 phase of neural progenitors and promotes efficient neuron production[END_REF][START_REF] Saade | Sonic hedgehog signaling switches the mode of division in the developing nervous system[END_REF][START_REF] Smith | Cell cycle and neuroepithelial cell shape during bending of the chick neural plate[END_REF][START_REF] Wilcock | Mitotic spindle orientation distinguishes stem cell and terminal modes of neuron production in the early spinal cord[END_REF]. Figure 6. Dynamics of the anterior neural tube growth fits with a homogenous lengthening of the cell cycle duration together with a proportionally declining progenitor population A-B, LOWESS curve fitting from our biological data to predict the evolution of the progenitor cell proportion over time. C, Simulations with a mixed proliferation regime (cell cycle of 8h50 for the first day and 13h30 the second day) with declining proportions of proliferating cells over time. Black curve and dots: 100% for the first 24h of simulation, 85% between 24 and 36h and 70% between 36 and 48h. Dark grey dots and curve: 100% for the first 24h of simulation, 85% between 24 and 48h. Light grey curve and dots are the same as the mixed regime with 100% proliferating cells displayed in Figure 5 panel E, provided for comparison. D, Simulations with a mixed proliferation regime and the declining proportions

  a confocal Zeiss 710 with Zen software and Leica SP8 with LasX software. Images were then processed and analyzed with FIJI. support by the Engineering and Physical Sciences Research Council (EPSRC) under grants no. EP/M006883/1 and EP/N014529/1, by the Royal Society International Exchanges under grant no. IE160750, by the Royal Society and the Wolfson Foundation through a Royal Society Wolfson Research Merit Award no. WM130048 and by the National Science Foundation (NSF) under grant no. RNMS11-07444 (KI-Net). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Figure 1 .

 1 Figure 1. Anterior trunk neural tube growth is slower along the anteroposterior axis than along the

Figure 2 .

 2 Figure 2. Anterior trunk neural tube anteroposterior length is restricted by interaction with the surrounding

Figure 3 .

 3 Figure 3. Adhesion of neural tube cells to the extracellular matrix is required for apicobasal elongation and

Figure 4 .

 4 Figure 4. Mitoses are randomly oriented in the neural tube between 2 and 4 days of development.

Figure 5 .

 5 Figure 5. Mean cell cycle length passes from 8h50 to 13h30 between 48h and 96h of development in the

Acknowledgements

We thank Dr. Bruno Glise for critical reading of the manuscript. ET acknowledges support from the Fondation pour la Recherche Médicale (FRM AJE201224), the Region Midi-Pyrénées (13053025), Toulouse Cancer Santé (DynaMeca), the CNRS and Université Paul Sabatier (UMR5547). PD acknowledges

of proliferating cells over time (100% for the first 24h of simulation, 85% between 24 and 36h and 70% between 36 and 48h) with either a narrow range around the long cell cycle values (+/-1h, black dots and curve) or a wide range around long cell cycle values (+/-5h, open circle and curve). All simulations were done 20 times, mean is plotted.


Alexandre Bonnefond

Olivier Simonin

Isabelle Guérin-Lassous

Extension of Flocking Models to Environments with Obstacles and Degraded Communications

In this paper, we study existing flocking models and propose extensions to improve their abilities to deal with environments having obstacles impacting the communication quality. Often depicted as robust systems, there is yet a lack of understanding how flocking models compare and how they are impacted by the communication quality when they exchange control data. We extend two standard models to improve their ability to stay connected while evolving in environments with different obstacles distributions. By taking into account the radio propagation, we model the obstacles impact on communications in a simulator that we use to optimize flocking parameters. The simulation results show the efficiency of the proposed models and how they adapt to different environmental constraints.

I. INTRODUCTION

Since its formulation by Craig Reynolds [START_REF] Reynolds | Flocks, herds and schools: A distributed behavioral model[END_REF] in 1987, the flocking theory has been of increasing interest. Based on three simple interactions between agents that are, the separation, the alignment and the cohesion (see Fig. 1), the author designed a decentralized system that can evolve in different types of environment. If this control strategy was originally designed for perception-based systems, the development of UAVs (Unmanned Aerial Vehicles) combined with improved communication abilities has made possible the adaptation of flocking to these types of vehicles. The agents are now able to exchange, via a high rate communication channel, crucial information such as their positions and velocities to their nearby "flockmates". As a consequence, most of the recent papers in the literature [START_REF] Olfati-Saber | Flocking for multi-agent dynamic systems: algorithms and theory[END_REF][3][4] [START_REF] Wen | A connectivitypreserving flocking algorithm for multi-agent dynamical systems with bounded potential function[END_REF] consider communicationbased flocking but without taking into account the real features of communications and how they are impacted by the presence of obstacles. The use of simple assumptions on communications appears to be one of the main issue when it comes to real experiments. The work of [START_REF] Fraser | Simulating the effect of degraded wireless communications on emergent behavior[END_REF] gives a good idea of how realistic communications can affect an emergent behavior. However, the controller involved in this paper is quite simple and hence doesn't reflect how standard flocking models are affected. The objective of this article is to understand how the flocking models compare and how they are impacted by the quality of the communication, in terms of received power, when exchanging control data, this quality varying according to the distribution of obstacles. We also propose to extend two standard models (Olfati-Saber [START_REF] Olfati-Saber | Flocking for multi-agent dynamic systems: algorithms and theory[END_REF] and Vásárhelyi/Vicsek [START_REF] Vásárhelyi | Optimized flocking of autonomous drones in confined environments[END_REF]) to improve their robustness *This work is supported by the Inria-DGA DynaFlock project. 1 INSA Lyon, CITI Lab., Inria Chroma team, 6 avenue des Arts, 69621 Villeurbanne, France firstname.lastname@inria.fr 2 Lyon 1 University, LIP Lab., Inria Dante team, 46 allée d'Italie, 69364 Lyon Cedex 07, France firstname.lastname@ens-lyon.fr Fig. 1: Flocking : The 3 interactions defined by Reynolds (illustration from [START_REF] Bouraqadi | Flocking-based multi-robot exploration[END_REF]) in the presence of obstacles while designing decentralized flocking. We aim to maintain properties such as cohesion (avoiding fragmentation), non collision, alignment, and quality of communication. For this purpose, we use the stochastic framework of [START_REF] Virágh | Flocking algorithm for autonomous flying robots[END_REF] to simulate and evaluate flocking models. We integrate in this simulator both realistic propagation models and various constrained environments including obstacles. Combined with an appropriate optimization process of the parameters, we show that the proposed flocking models remain robust even in highly constrained environments.

In Section II, we go through a non exhaustive list of existing flocking models with a focus on two specific models that we use in this article. In Section III, we describe the extension of the flocking models we propose. Then in Section IV, we present the models we consider for communications before redefining the notion of neighbourhood. We then compare these new models in Section V within different environments under varying constraints, specifically in terms of obstacles. We finally conclude with a discussion on how this framework could be extended.

II. EXISTING FLOCKING MODELS

Flocking models are usually applied to multi-agent systems in which a collective motion is aimed. Those models rely on three fundamental interactions introduced by Reynolds [1]:

• Collision Avoidance (or separation): avoid collision with nearby flockmates • Velocity Matching (or alignment): attempt to match velocity (heading and amplitude) with nearby flockmates • Flock Centering (or cohesion): attempt to stay close to nearby flockmates where the flockmates are any agent in the system. Models are usually composed of N agents. Each agent i is referred to by its position q i and velocity v i (bold notation is used for vectors). It also has its set of neighbours with whom it interacts, denoted by:

N i = {j ∈ [[1, .., N ]]; j = i : | q j -q i ≤ r} ( 1
)
where r is the interaction range. Most of the flocking models use these rules and definitions as a starting point, but when it comes to dynamics and control, different models are proposed. Among all the flocking models found in the literature, multiple strategies have emerged such as the leader-follower method [START_REF] Gamage | Leader follower based formation control strategies for nonholonomic mobile robots: Design, implementation and experimental validation[END_REF][5], the potential field method [START_REF] Ge | Queues and artificial potential trenches for multirobot formations[END_REF], the behaviorbased method [START_REF] Balch | Behavior-based formation control for multirobot teams[END_REF], the virtual structure method [START_REF] Li | Formation uav flight control using virtual structure and motion synchronization[END_REF] or the reinforcement learning method [13][14].

Because the flocking is now used to tackle problems with different goals, Reynolds' laws are embedded into complex multi-agent dynamics such as in the paper of Olfati-Saber [START_REF] Olfati-Saber | Flocking for multi-agent dynamic systems: algorithms and theory[END_REF] (denoted O-S hereafter). In his work, the author designs an appropriate mathematical framework to ensure the stability of his model. This aspect is important in the control theory and ensures a non erratic behavior which is highly desired for UAV-based applications for instance. O-S also introduces features like obstacle avoidance and a socalled "navigational feedback" giving a global direction to the agents, key parameter to avoid the fragmentation of the flock (i.e. the separation into multiple clusters). For these reasons, this model has been repeatedly used in the literature [START_REF] La | Flocking control of a mobile sensor network to track and observe a moving target[END_REF][14] [START_REF] Wen | A connectivity-preserving flocking algorithm for nonlinear multi-agent systems with bounded potential function[END_REF] with slight modifications and different purposes. The dynamics used in O-S is based on a double integrator with a velocity consensus protocol:

qi = v i vi = u i , (2)
and the control input u i , applied to each agent i, relies on three terms:

u i = f g i + f d i + f γ i , (3)
where f g i is a gradient-based term merging both separation and cohesion, f d i is the velocity alignment term which acts as a damping force and f γ i is the so-called navigational feedback. This latter term acts as if there was a virtual leader whose parameters (location and velocity) are known by all the agents. As shown in Section III, we use part of the O-S's model as a first layer to build a new model that satisfies our objectives. Fig. 2 shows the regular lattice agents achieve when applying O-S's model.

Another flocking mechanism that we use in this paper is the one introduced by Vicsek [START_REF] Vicsek | Novel type of phase transition in a system of self-driven particules[END_REF] in 1995, which considers that the particles or agents are "self-propelled". This characteristic, derived from ferromagnetic material observations, will be very important in our consideration as it leads to a pure flocking (i.e. purely decentralized). Based on this idea of self-propelling particles, Vásárhelyi et al. developed a framework [START_REF] Vásárhelyi | Optimized flocking of autonomous drones in confined environments[END_REF] including a simulator which goes deep into the physical representation of the environment and the related constraints such as communication delays, inertia, refresh rate of sensors, locality of the communication, inaccuracy of on-board sensors and outer noises. The model is defined by an equation determining the acceleration a i of each agent i (see [START_REF] Virágh | Flocking algorithm for autonomous flying robots[END_REF]):

a i = η i + v d i -v i -v s i |v d i -v i -v s i | • min v d i -v i -v s i τ CT RL , a max , (4)
where η i models the outer noise (like the wind effect), v s i represents the inner noise on the velocity measurements (due to GPS approximations), τ CT RL is the relaxation time of the velocity controller, a max is the limit in acceleration and v d i is the desired velocity:

v d i = ṽd i |ṽ d i | min{|ṽ d i |, v max }, (5)
where v max is the maximal allowed speed of each agent and ṽd i is defined as:

ṽd i = v i |v i | v f lock +v rep i +v f rict i + s v wall is + s v obstacle is ( 6
)
where v i is the real velocity of the ith agent, v f lock is the preferred speed of the agents, v rep i is the repulsion term based on the half spring model, v f rict i is the alignment term that synchronizes motion within the flock, v wall is and v obstacle is are the terms used for obstacles and walls avoidance. The combination of this control law as well as the integration of real outdoor experiment constraints makes this framework very attractive. This is why we use it in this paper to propose extended versions.

Although those different models reach good performances (like, for instance, stability, robustness, smooth behaviour) in simulation, they have issues when scaling up to real outdoor experiments. Indeed, as pointed out in [START_REF] Vásárhelyi | Optimized flocking of autonomous drones in confined environments[END_REF], in practice, multiple communication outages happened and it is a major drawback when working with UAVs.

III. FLOCKING MODELS EXTENSION

In order to avoid the pitfall experienced by the models presented in Section II when used in a realistic environment, we propose new flocking models. Those models exploit the fundamental interactions of the two models described previously and include additional features expected to be robust to degraded communications due to the presence of obstacles.

A. VAT Model: Vásárhelyi + Attraction

Backed by a robust behaviour, the velocity controller defined by Vásárhelyi (see Eq. 6) was originally designed without any attraction term. We observed that, within environments containing many obstacles, the repetitive obstacles avoidance maneuvers combined with degraded communications eventually lead to fragmentation. We hence decided to introduce an attraction term based on the half-spring model in order to enhance the cohesion whenever it is possible:

v att ij =    p att (r att 0 -q ij ) q i -q j q ij , if q ij > r att 0 0, otherwise (7)
where p att is a linear gain, q ij = |q i -q j | is the distance between agents i and j and r att 0 is the attraction range under which agents stop attracting each other so as to prevent additional collisions. Also we have v att i = j v att ij . Adding this term to the input control must be done carefully regarding the repulsion term in order to not have oscillations. Hence, we set r att 0 large enough for the agents to have a slack distance between attraction and repulsion. Our first proposition is hence a velocity-based model denoted VAT:

ṽd i = v i |v i | v f lock + v rep i + v att i + v f rict i + s is (8)
The obstacle and wall interactions are concatenated in the last term for convenience and Eq. 5 remains valid. This resulting desired velocity encompasses 12 parameters, acting on different terms of the law and producing different behaviours. We will see in Section V-B how those parameters are optimized to reach an optimal behaviour (see the attached video).

B. VOS Model: Vásárhelyi + Olfati-Saber

The second model we extend, denoted VOS, is based on a combination of O-S's model, exploiting the stability of the behaviour derived from using its control input (3), and some features of Vásárhelyi's model used to add self-propelling and obstacle avoidance abilities. To this end, we come up with a nested control structure merging both controllers:

     u i = αf g i + βf d i ( 9a
)
ṽd i = v i |v i | v f lock + s v wall,obst is + Ts u i dt (9b)
where α and β are some linear coefficients which compensate the delay introduced by the integration and T s = 1/f s with f s = 100Hz the working frequency of the controller. The other terms have been defined in Section II. The objective here is to avoid the fragmentation by adding the selfpropelling term (v f lock .v i /|v i |). We recall that in its original paper [START_REF] Olfati-Saber | Flocking for multi-agent dynamic systems: algorithms and theory[END_REF], O-S introduces a navigational feedback term to deal with the fragmentation (f γ i in Eq. 3). It is something we want to avoid here as it goes against our will of developing a fully decentralized system performing pure flocking. Even if we had kept the f γ i term, the stability of O-S's model would not be guaranteed because it relies on free-flocking [START_REF] Olfati-Saber | Flocking for multi-agent dynamic systems: algorithms and theory[END_REF], i.e. flocking in environments without obstacles. However we will see in Section V-C that VOS presents a good stability even with many obstacles. As for the VAT model, the 10 parameters of the VOS model will be optimized before evaluation (see Sec. V-B).

IV. COMMUNICATION AND NEIGHBOURHOOD IN PRESENCE OF OBSTACLES

As mentioned earlier, most of the flocking simulations consider ideal communications and do not take into account the impact of the environment on the communication quality which usually results in issues when scaling up to real experiments. A majority of the simulations uses a distance threshold under which the communication is possible and above which it is not (parameter r of Eq.1) and the communication is either perfect [START_REF] Olfati-Saber | Flocking for multi-agent dynamic systems: algorithms and theory[END_REF] [START_REF] La | Flocking control of a mobile sensor network to track and observe a moving target[END_REF], or more realistic (with added random delays and packet losses) [START_REF] Vásárhelyi | Optimized flocking of autonomous drones in confined environments[END_REF] [START_REF] Virágh | Flocking algorithm for autonomous flying robots[END_REF]. In [START_REF] Li | Flocking control for multi-agent systems with communication optimization[END_REF], the authors use the BER (Bit Error Rate) combined with the SNR (Signal to Noise Ratio) to characterize the communication quality and integrate it into a performance index expressing a trade-off between inter-agent distance and quality of the communication. In [START_REF] Frew | Hybrid rf propagation model using itm and gaussian processes for communication-aware planning[END_REF], the authors come up with an hybrid system combining a prediction step and a learning step of the path loss field. They use the SPLAT!1 tool to predict how the path loss field is affected by the topology of an area. Although this method gives realistic insights for the RF measurements, it does not tell how the multi-agent system is impacted by these fluctuations. On this subject, in [START_REF] Fraser | Simulating the effect of degraded wireless communications on emergent behavior[END_REF], the authors explain how an emergent behavior, such as flocking, is impacted when we use a relatively accurate propagation model and tune the inherent losses due to shadowing effects, SNR, jamming and so on. However, these studies consider only simplistic flocking models, far from realistic fleet of UAVs.

We propose to use a communication model based on the path loss estimation, including the presence of obstacles. This allows us to redefine the neighbourhood definition and consequently to better estimate with whom each agent interacts.

In the following, we introduce a prediction step using the log distance path loss model (denoted LDPL) for both free space propagation but also inside obstacles and we see how heterogeneous environments affect the quality of the communications and consequently the flocking strategy.

A. Log Distance Path Loss model with environment losses

Instead of using distance-based communication model, we aim at modelling how radio waves dynamically evolve with the environment and affect the neighbourhood. To do so, we consider the LDPL model that predicts the path loss a signal encounters within different types of environments [START_REF] Rappaport | Wireless Communications: Principles and Practice[END_REF]. It is formulated as follows: where, P L is the path loss measured in decibels (dB), P T x dBm is the transmitted power in dBm at the emitter, P Rx dBm is the received power in dBm at the receiver, P L 0 is the path loss at the reference distance d 0 calculated using the Friis free-space path loss model, d is the distance between the emitter and the receiver, γ is the path loss exponent that depends on the environment (we use γ = 2 in free space) and X g reflects the variations of the path loss caused by shadowing effects and/or multiple paths. This latter is a zero mean Gaussian random variable with a standard deviation of σ = 2dB. For the sake of clarity, we will refer to this communication model as regular. The green curve in Figure 3 represents the evolution of P Rx dBm on a path where there is no obstacle and P T x dBm = 20dBm. We can see that it takes approximately 300 meters to reach the power threshold P min = -70dBm which is often considered as the minimum level required to enable a reliable communication between agents in Wi-Fi for instance.

P L = P T x dBm -P Rx dBm , (10)
P L =    P L 0 + 10γlog 10 d d 0 + X g , if d ≥ d 0 P L 0 , otherwise (11)
In the previous formula, the environment is modeled via, among others, the parameter γ and X g . In order to better take into account the impact of the different obstacles on the communication quality, we add a dissipation model inside the obstacles leading to high loss of the signal strength when the signal goes through an obstacle. The resulting path loss is expressed as follows:

P L =              P L 0 + 10γlog 10 ( d -d obst d 0 ) + X g + P L obst (d obst ), if d -d obst ≥ d 0 P L 0 + P L obst (d obst ), otherwise (12)
where, d obst is the length of the path inside the obstacle and P L obst is the path loss in dB due to obstacle attenuation. By using a simplified model of the path loss, this latter is defined as:

P L obst (d obst ) = 10γ obst log 10 (d obst ) + K (13)
where γ obst is the path loss exponent for the obstacle and K is a constant. We assume that d obst ≥ 1m. This communication model will be called degraded. The orange and the blue curves in Figure 3 represent the evolution of P Rx dBm on a path where respectively d obst = 1.5m and d obst = 3m. In this configuration we have γ obst = 4, P T x dBm = 20dBm and K = 0. We can see that this model has a strong impact on the communication's quality as it takes less than 30 meters to first lose the communication among agents with a 3 meters obstacle between them.

B. Neighbour Filtering

Thanks to Eq. 12, we are now able to redefine the neighbourhood N i seen by agent i:

N i = {j ∈ [[1, .., N ]]; j = i : P L ij ≤ P threshold } ( 14
)
where P L ij is the path loss (Eq. 12) measured from agents j to i and P threshold ( = P min ) is the power threshold above which the communication is no longer possible. Also we introduce another threshold N max , adjustable, which acts on the cardinality of N i :

N i = N i , if |N i |≤ N max N f iltered i , otherwise (15)
where N f iltered i corresponds to a subset of N i where the worst neighbors in terms of path loss are recursively removed until |N i |= N max . This feature is very important as it limits the computational costs of each agent while focusing the interactions on the most reliable flockmates. To add more realism, this neighbourhood is updated at the same frequency as the GPS since position information are needed to compute inter-distances. With this new definition of the neighbourhood, based on a more realistic communication model, we can study the evolution of a flock in highly constrained environments with multiple obstacles.

V. COMPARISON OF MODELS

A. Simulation framework

The simulator we use is an extension of the tool2 developped by Viragh et al. [START_REF] Virágh | Flocking algorithm for autonomous flying robots[END_REF], illustrated in Fig. 4. It is a multiagent based simulator which introduces random variables for modelling noises on different variables, adding a complexity for both the resolution of the stochastic differential equations (SDE), using Euler-Marayuma method, but also the parameter optimization (see Section V-B). We decided to use this simulator instead of off the shelf simulators such as ARGoS or GAZEBO for two main reasons. The first one is that the original Vásárhelyi's model was already implemented in this simulator and hence would give us a fair comparison with our models. The second one, and probably the most important, is that we mostly use this simulator on servers to perform parameters optimization which requires hundreds of runs and this tool is well designed for that. The visualisation is only used to check the behavior of the agents for a given set of parameters. The hyper-parameters of an experiment such as its duration T or its calculation step for the Euler method are adjustable as well as the obstacle distribution defined by, among others, its density D which we measure as follows:

D = 100 × s A s A arena , (16)
where A s is the area of the sth obstacle and A arena is the area of the arena. Figure 4(b) shows a 3D view of a so-called "Forest" environment as it has many small obstacles.

Several features have been added to this simulator to comply with our needs like the clusters detection using Depth First Search algorithm (DFS) over the graph of connected agents, the ability to switch on the fly the flocking model or the communication model (regular/degraded). See the attached video illustrating models simulation.

B. Parameter optimization

As we mentioned in Section III, the models we propose and study have several parameters that must be optimized to reach their best possible flocking behaviour. We took some inspiration from [START_REF] Vásárhelyi | Optimized flocking of autonomous drones in confined environments[END_REF] to define different metrics so as to rate a flocking model:

• Velocity: Ψ vel = 1 T 1 N T 0 N i=1 |v i (t)|dt,
which characterizes the average velocity of the flock.

• Collision (between agents):

Ψ col = 1 T 1 N (N -1) N 0 N i=1 j =i Θ(q ij (t) -r coll )dt,
evaluates the ratio of collision, with Θ the Heaviside step function and r coll the minimum allowed distance between two agents.

• Connected:

Ψ clust = 1 T T 0 |O(t)|dt, O(t) = {i ∈ [[1, .., N ]] : |N i |= 0}, O(t)
is the set of connected agents at t, ie. the complement of isolated agents O(t). • Wall (collision):

Ψ wall = T 0 N i=1 Θ(q is (t))q is (t)dt T 0 N i=1 Θ(q is (t)dt
, where qis takes positive values outside the arena and inside an obstacle, negative ones otherwise. Ψ wall characterizes the average distance spent inside the obstacles.

• Cluster Correlation:

Ψ corr = 1 T 1 N T 0 N i=1 1 |Ni|-1 j∈Ni vi•vj |vi||vj | dt,
which gives a measure on the average alignment per cluster.

• Cluster Received Power:

Ψ pow = 1 T 1 N T 0 N i=1 1 |Ni|-1 j∈Ni P Rxij (t)dt,
which gives a measure on the average received power per cluster.

Based on this, the following conditions are required to obtain an optimal flocking behaviour: Ψ vel → v f lock , Ψ col → 0, Ψ clust → 1, Ψ wall → 0, Ψ corr → 1, Ψ pow ≥ P min . Those metrics are then passed into non-linear functions to map them between 0 and 1 (the sigmoid function for example), 0 being the worst case and 1 the best, to fit into a multi-objective optimization framework.

In this paper, we used NSGA-III [START_REF] Deb | An evolutionary many-objective optimization algorithm using reference-point-based nondominated sorting approach, part i: Solving problems with box constraints[END_REF][21] which is a fast non-dominated sorting genetic algorithm designed for multiobjective problems. This algorithm is however not suited for stochastic optimization but we modified it by evaluating the same solution several times and using the computed variances of the metrics as a constraint of the optimization. This leads to more robust solutions, less impacted by the randomness of an experiment.

Our strategy is to optimize both VAT and VOS models on a "generic environment", i.e. with an average density of obstacles, and with the degraded communication model. The number of agents used for the optimization step and the simulations whose results are given in Section V-C is the same: N = 10. Then we use those optimized versions on specific obstacle distributions to compare their performances. The generic environment has N obst = 90 obstacles uniformly distributed on a 1250 m 2 arena. Each obstacle is represented by a square whose edge length L follows a normal (Gaussian) distribution L ∼ N (µ, σ), with mean µ = 24m and a standard deviation σ = 12m. The average density D = 3.3%. The optimization process is represented in Fig. 5 which tracks the evolution of the average of the metrics of the Pareto set at each generation for both models (VOS in blue and VAT in red). The color-bar gives an index on how advanced we are in the optimization, based on a maximum number of generations. Although VOS seems to be outperformed by VAT, it is only because its Pareto front is bigger with broader elements than VAT one. Table I gathers the average metric values after 100 evaluations in the generic environment with degraded communication using the "best" solutions derived from the optimization for both models. One can observe that both models have close performances except for the velocity metric which is lower in VOS model. These solutions are carefully compared in the following subsection.

C. Simulation results

We aim here at comparing the VAT and VOS models inside different environments and using the different communication models we proposed (regular and degraded). The following comparison will be based on two different environments (see Fig. 4):

1) Forest: N obst = 500; L ∼ N (4, 2); D = 0.7%.

2) City: N obst = 15; L ∼ N (200, 20); D = 36.5%. We will also make a comparison to Vásárhelyi's model, denoted VA, that we simulate using its optimal set of parameters given in [START_REF] Vásárhelyi | Optimized flocking of autonomous drones in confined environments[END_REF]. Tables III and IV list the results from 100 stochastic evaluations with the optimized set of parameters for each model. The average values of the mapped metrics (1 is best, 0 is worst) as well their relative variances are displayed. The last line corresponds to the average inter-agent distance in meters. To remain consistent, we use the same time of experiment T = 600s and the same accuracy of the Euler method DeltaT = 0.01s. In addition, we consider a delay t d = 0.2s in the communications (ie. access delay to the radio medium, delay of data packet transmission and propagation delay) and reckon that agents exchange information at the GPS frequency f GP S = 5Hz.

The first important observation we can draw here is that the VA model has much lower performances than the original results of [START_REF] Vásárhelyi | Optimized flocking of autonomous drones in confined environments[END_REF] reported in Table II. This can be explained by the fact that the author relies on small arenas in his simulations, leading to an indirect effect of cohesion, the agents having no choice but to head back toward the center of the arena. In our considerations, and to avoid this advantageous behaviour, we set a large enough value (1250m) for our arena length. Also the distribution of the obstacles is not explicitly mentioned in [START_REF] Vásárhelyi | Optimized flocking of autonomous drones in confined environments[END_REF] which obviously impacts most of the metrics. Finally, taking into account a radio propagation model also adds a new constraint to the model that was not considered in the initial work. The drawbacks of not adding an attraction term are even more visible on the forest environment where the cluster parameter is relatively low, meaning there has been some fragmentation (Clust=0.637 in Table III). One may also wonder why adding an attraction term for the VAT model does not create more collisions. This can be explained by the fact that this attraction term is only effective above a precise inter-agent distance with a bounded resulting inertia preventing oscillations and hence collisions. The fact that the VAT model is also designed and optimized for degraded communications, and therefore applies the neighbour filtering, is one more reason to explain VA and VAT differences. 1) Maintaining the Connection: One of our goal in this article was to implement models capable of preventing the undesired effect of perpetual fragmentation. We can see in both tables that this objective is mostly satisfied. Indeed, for both VAT and VOS, the average cluster metric is high, 0.868 ≤ Clust ≤ 0.985. This means that even if fragmentation occurred, the resulting clusters managed to repair broken connections and got back to a single cluster. Partial fragmentation was to be expected especially in the forest environment as obstacles create frequent communication outages and interaction interruptions 2) Received Power and Collision: The inter-agent distance is a resulting parameter of all interactions and is mostly optimized by the received power metric but also the collision one. Agents must lie within a restricted range of distances between each other in order to be far enough to avoid collisions but close enough to remain connected. Please note that in these tables, a low value of the received power doesn't mean that the connection is lost. For instance in Table IV, we have Pow = 0.573 for the VOS model with degraded communications, but this is equivalent to a value of P R = -59.4dBm, meaning that the communication was possible on average. As a comparison, in Table III, for VAT model with degraded communications, Pow = 0.045 -71.7dBm (-70dBm being the lowest acceptable value). This confirms the fact that, even if the communication is mostly broken, whenever they can, agents manage to get back to a nominal flocking behaviour.

3) Environment Overfitting: As mentioned earlier, the models we simulate are not optimised on the environments we use for the comparison. Therefore, there is a risk of overfitting a particular obstacle distribution. We can observe it for the VOS model in Table III where the average velocity metric is low, 0.032 ≤ Vel ≤ 0.055, equivalent to 60% of the desired velocity v f lock . Although, in Table I, Vel= 0.612, the number of obstacle is around 5 times smaller. Within the forest environment, there are many obstacles close to each other leading to a "cage effect" and hindering the flock movements. It also creates many oscillations 0.531 ≤ Corr ≤ 0.54 and potential collisions 0.793 ≤ Col ≤ 0.799 (see video). The VOS model is hence not robust to the variation of obstacle distribution.

4) Regular/Degraded: The resulting effects of the degraded communication model can be observed in both tables. Moreover, the last line which concerns the interagent distance shows that this mode will increase the average distance between agents and even more when in the forest environment. However, the VAT model seems to be very resilient and robust to losses introduced by the degraded mode on its received power which makes it a valuable model.

VI. CONCLUSION

In this article, we proposed two flocking models based on the state of the art (VAT Vásárhelyi + Attraction and VOS Vásárhelyi + Olfati-Saber) to deal with degraded communications in presence of obstacles. We compared models via simulations integrating a Log Distance Path Loss model that takes into account the propagation inside the obstacles. With this evaluation environment, it was possible to analyse how flocks are affected or robust to degraded communications in complex environments. From the results, we can claim that the VAT model shows more encouraging performances than VOS and the original Vásárhelyi's model.

Even if the optimization has been conducted here on a single "generic" environment, it could be interesting to systematically optimize the models in each used environment. We then plan to design adaptive flocking strategies capable of optimally evolving in various environments with different types of communication. We also plan to carry out experiments with quadrirotor UAVs, in order to measure model efficiency in real environments.

Fig. 2 :

 2 Fig. 2: Illustration of the Olfati-Saber model in the simulation framework, with 20 agents and communication links in white. The top green bar indicates the spatial scale.

Fig. 3 :

 3 Fig. 3: Evolution of the received power as a function of interagent distance. The green line represents a path without obstacle, the orange line has a 1.5 meters long obstacle on its path and the blue one has a 3m long obstacle on its path

  (a) 3D View of a "City" environment (b) 3D View of a "Forest" environment

Fig. 4 :

 4 Fig. 4: Simulator views

Fig. 5 :

 5 Fig. 5: Evolution of the metrics during optimization using degraded communication in the genetic env. (VAT in red, VOS in blue)

TABLE I :

 I VAT and VOS optimized parameters in GENERIC env.

		Vel	Col	Clust	Wall	Corr	Pow
	VAT	0.999	0.99	0.94	1	0.86	0.202
	VOS	0.612	0.88	0.92	1	0.784	0.26

TABLE II :

 II VA performances from[START_REF] Vásárhelyi | Optimized flocking of autonomous drones in confined environments[END_REF]

		Vel	Col	Clust	Wall	Corr	Pow
	VA	0.938	0.945	1	0.997	0.916	N/A

TABLE III :

 III FOREST

			Regular			Degraded
		VA	VAT	VOS	VA	VAT	VOS
	Vel	0.975	0.946	0.032	0.975	0.955	0.055
	σ	0.010	0.03	0.052	0.008	0.023	0.063
	Col	0.301	1	0.793	0.286	0.961	0.799
	σ	0.429	0.099	0.382	0.429	0.172	0.380
	Clust	0.637	0.98	0.923	0.508	0.938	0.868
	σ	0.148	0.085	0.13	0.118	0.121	0.163
	Wall	1	1	1	1	1	1
	σ	0	0	0	0	0	0
	Corr	0.133	0.712	0.54	0.128	0.698	0.531
	σ	0.068	0.106	0.12	0.06	0.12	0.112
	Pow	0.033	0.53	0.517	0	0.045	0.031
	σ	0.049	0.11	0.074	0	0.079	0.066
	Dist	132.3	37.9	49.9	132.5	42.11	54.37

TABLE IV :

 IV CITY

			Regular			Degraded
		VA	VAT	VOS	VA	VAT	VOS
	Vel	0.999	1	0.484	0.999	1	0.501
	σ	1.1e-5	0	0.095	4.8e-5	0	0.101
	Col	0.71	0.999	0.918	0.654	1	0.86
	σ	0.442	0.09	0.264	0.47	0	0.33
	Clust	0.671	0.985	0.939	0.54	0.915	0.925
	σ	0.18	0.065	0.142	0.197	0.17	0.165
	Wall	1	1	1	1	1	1
	σ	0	0	0	0	0	0
	Corr	0.605	0,912	0.886	0.73	0.945	0.894
	σ	0.165	0.067	0.082	0.112	0.04	0.073
	Pow	0.512	0.797	0.684	0.335	0.777	0.573
	σ	0.205	0.092	0.073	0.286	0.16	0.204
	Dist	82.19	29.4	44.3	86.76	33.8	45.55

Information about SPLAT!: https://www.qsl.net/kd2bd/splat.html

Access to the original project here: https://github.com/csviragh/robotsim