Etienne Loiseau

Simon Gsell

Aude Nommick

Charline Jomard

Delphine Gras

Pascal Chanez

Umberto D'ortona

Laurent Kodjabachian

Julien Favier

Annie Viallat
email: annie.viallat@univ-amu.fr

& A Viallat

Active mucus-cilia hydrodynamic coupling drives self-organization of human bronchial epithelium

Introduction

The airways are protected from inhaled pollutants and pathogens by a layer of mucus, a complex fluid that is continuously transported along the bronchial walls before being cleared out and swallowed. The mucus is propelled via the continuous beating of active cilia exposed by the epithelial ciliated cells. This process, called mucociliary clearance, is strongly impaired in chronic respiratory diseases, such as primary dyskinesia, severe asthma [START_REF] Chanez | Severe asthma is an epithelial disease[END_REF] , chronic obstructive pulmonary disease 2 and cystic fibrosis [START_REF] Boucher | Airway surface dehydration in cystic fibrosis: pathogenesis and therapy[END_REF] , which affect hundreds of millions of people [START_REF] Burney | The global burden of chronic respiratory disease in adults[END_REF] .

Effective directional transport of mucus requires both a high surface density of active [START_REF] Matsui | Coordinated clearance of periciliary liquid and mucus from airway surfaces[END_REF][START_REF] Khelloufi | Spatiotemporal organization of cilia drives multiscale mucus swirls in model human bronchial epithelium[END_REF] epithelial cilia to overcome friction forces on the epithelial wall and a directional coordination of ciliary beats. Recently, centimetric mucus swirls and a circular orientational order of the underlying ciliary beats [START_REF] Khelloufi | Spatiotemporal organization of cilia drives multiscale mucus swirls in model human bronchial epithelium[END_REF] were observed in bronchial epithelia reconstructed in vitro only when the ciliary density was high [START_REF] Matsui | Coordinated clearance of periciliary liquid and mucus from airway surfaces[END_REF][START_REF] Khelloufi | Spatiotemporal organization of cilia drives multiscale mucus swirls in model human bronchial epithelium[END_REF] . Yet, the relationship between cilia density and the spatial scale on which mucus is transported is unknown. Furthermore, the biological and physical mechanisms involved in the emergence of the orientational order of ciliary beats during ciliogenesis, its stability and long-term repair are much debated. The ciliary beat direction is set by the polarity of the basal body of the cilium and is considered to be fixed during development by the planar cell polarity (PCP) proteins that regulate the polarization of epithelial tissue [START_REF] Spassky | The development and functions of multiciliated epithelia[END_REF][START_REF] Walentek | Planar Cell Polarity in Ciliated Epithelia[END_REF][START_REF] Boutin | Biology of multiciliated cells[END_REF] . However, daily variations of ciliary-driven cerebrospinal fluid flow patterns, recently reported in epithelia of mouse brain ventricles, showed that ciliary beats can change direction [START_REF] Faubel | Cilia-based flow network in the brain ventricles[END_REF] . Moreover, under an external fluid flow, the ciliary beat directions on the skin of the Xenopus embryo before PCP is in place and on maturing cell cultures of mouse ependyma redirected along the flow streamlines [START_REF] Guirao | Coupling between hydrodynamic forces and planar cell polarity orients mammalian motile cilia[END_REF][START_REF] Mitchell | A positive feedback mechanism governs the polarity and motion of motile cilia[END_REF] . Cilia can then actively respond to external hydrodynamic cues. Interestingly, the bronchial epithelium is naturally continuously subjected to the endogenous mucus flow generated by ciliary beats. This raises crucial questions: could mucus flow affect the long-range ordering of ciliary beats by creating a hydrodynamic coupling among cilia? Would hydrodynamic feedback on cilia lead to an enhancement of orientational ciliary beat order, and in turn result in a better mucus transport?

Here, we address these questions by combining experiments on in-vitro reconstituted human bronchial epithelia and a hydrodynamic model. We observe and quantify the dynamics of formation and growth of swirly patterns of mucus flow during ciliogenesis. We show that, at a physiological cilia density, a macroscopic swirl of mucus develops associated with a circular coordination of ciliary beat directions and a circular pattern of planar cell polarity of the tissue.

By establishing that the macroscopic ciliary-beat order is lost and recovered by respectively removing and adding mucus, we demonstrate that cilia/mucus hydrodynamic interactions drive the collective dynamics of ciliary-beat directions. Based on these observations, we developed a two-dimensional hydrodynamic model which highlights two relevant physical parameters of the epithelium, the density of cilia and mucus viscosity, allowing the establishment of a phase diagram of mucus transport.

Results

The reconstituted bronchial epithelium is cultured in-vitro at the air-liquid interface from human primary cells (Fig. 1a). The pseudo-stratified tissue consists of basal cells, goblet cells and ciliated cells each carrying 200-300 active cilia gathered in bundles on the apical surface [START_REF] Vladar | Airway epithelial homeostasis and planar cell polarity signaling depend on multiciliated cell differentiation[END_REF] .

After cellular differentiation, the tissue exhibits the characteristics of the bronchial epithelium [START_REF] Gras | An ex vivo model of severe asthma using reconstituted human bronchial epithelium[END_REF] such as active cilia which beat at a frequency around 8 Hz and the production and transport of mucus on the epithelium surface. In the culture chamber, the basal side is in contact with the culture medium through a porous membrane and the apical part is at the air interface.

Mucus flows and beats of bundles of cilia are measured by videomicroscopy.

Ciliogenesis dynamics: emergence and growth of mucus swirls and ciliary-beat coordination.

Density of active cilia:

The percentage of the culture area covered by the active cilia, φ, was measured daily from the full confluency of the cell culture and throughout ciliogenesis on 10 cultures from 5 donors (see Fig. 1b and Table S1). After eight days, φ is lower than 5% (Fig. 1b). Then, the ciliary density increases by successive steps of 20 to 25% (Fig. 1c) to reach a final value of about 70-80% after two more weeks (Fig. 1d). Interestingly, the dynamics of ciliogenesis is almost identical for all 10 cultures (Fig. 1b). The temporal evolution of the spatial distribution of the ciliated patches is shown in Extended Data Fig. 1. We observe the appearance of a percolating cluster for a cilia density above 60 %, in coexistence with clusters of finite size. This analysis shows that ciliated cells appear at random positions during ciliogenesis, and not by forming organized structures of large circular or compact clusters.

Spatial range of mucus transport:

From day 8 to day 11, when φ increases from 3% to 15%, the mucus layer above the cilia vibrates but is not transported (see Movie S1). A local mucus transport is observed on day 12 (Fig. 1e and Movie S2), for φ ∼ 25%. Strikingly, from day 13 and 14, while φ remains almost constant, mucus swirls of a few hundred micrometers in diameter (Fig. 1f, g and Movies S3 and S4) are observed, together with small paths along which mucus is transported towards the edges of the culture. When φ gradually increases, the number of swirls decreases but their diameter increases by capturing the mucus transported nearby (Fig. 1 h & i, Movie S5) and reaches the millimeter scale for φ in the 35-50% range.

Generally, we observe simultaneously the presence of 3 swirls of several 100s µm in diameter per sample. The size of the culture chamber, 5 mm in diameter, limits the number of swirls that can coexist at the same time. A typical example consists of 1 swirl of 1.3 mm in diameter and 2 smaller swirls of ~600 µm in diameter in sample #8 on day 15. For φ greater than 50%, a single mucus swirl extends over the entire surface of the culture (Fig. 1i and Movie S6). The characteristic size of mucus swirls during ciliogenesis is indicated in Fig. 1b. One should notice that swirls rotate clockwise or counterclockwise with equal probability (see Movie S7). To understand the emergence of these fluid flow patterns, we examined how microscopic cilia self-organize and coordinate their beating underneath the mucus.

Cilia beat frequency: An analysis of the frequencies on a large surface reveals a distribution with a mean frequency of f= 7.6 +/-1.3 Hz. This results in a spatial cut-off, illustrated on a frequency map (Extended Data Fig. 2), which sets a maximum length scale (order of 8-10 cells) for the existence of a phase correlation. Furthermore, even at short range, phase correlation is very unstable over time (Extended Data Fig. 2). This prevents the propagation of metachronal waves over long distances and in the following we will focus on the dynamics of alignment interactions at long distances.

Orientational order of the ciliary beats: the beat direction of each bundle of cilia is represented by a rod in Fig. 2a and 2b (see also Movie S8). As illustrated in Fig. 2b, the directions of the ciliary beats present a strong circular order under the mucus swirls, with only a few localized defects mainly located in the swirl center. Areas where ciliary beats are aligned are also observed (Fig. 2a). The beat orientational order is quantified in a square surface of length L by 2 parameters classically used for nematic materials [START_REF] Kawaguchi | Topological defects control collective dynamics in neural progenitor cell cultures[END_REF][START_REF] Saw | Topological defects in epithelia govern cell death and extrusion[END_REF] . The first one is a local order parameter S(L) = <cos 2θ>² + <sin 2θ>², where θ is the angle between the beat direction and a reference direction, and the average < > relates to all bundles in the domain of size L. S varies between 0 and 1; 0 and 1 are respectively an isotropic state with random beat directions and a fully ordered state with aligned beat directions. The second parameter is the density of topological defects defined as local discontinuities in the director field of beat directions (Fig. 2a). They are associated with topological charges + ½ (Fig. 2a red circles) or -½ (Fig. 2a blue triangles), corresponding to a clockwise or counterclockwise π-rotation of the director field around the singularity respectively.

The ciliary-beat orientational order under a mucus swirl is computed at different length scales L. S(L) is averaged on a set of square tiles of length L that totally pave the mucus swirl, and is compared to S(L) computed for a full order when all ciliary-beat directions are tangential to the swirl. The comparison of the two curves shown in Fig. 2c for a typical ciliary beat 100 µmvortex reveals a strong tangential orientational order of the ciliary beats.

At the end of ciliogenesis, below the macroscopic mucus swirl, we observed in its central part a strong circular orientational order of ciliary beats originating from a previous local swirl that has grown overtime. Yet, farther from the center, at a distance of 1.5 mm, the ciliary-beat orientational order is low (S~ 0.1) (Fig. 2d). This order measured at different length scales (L= 50, 100 and 200 µm) increases very significantly within four days before reaching a plateau where S ranges between 0.4 to 0.55 (Fig. 2d). Simultaneously, the surface density of topological defects decreases threefold (Fig. 2e). These results clearly show that the ciliarybeat directions continue to coordinate for several days on mature cultures. It should be noted that during this increase in directional order, mucus rotates on the epithelial surface where it exerts a continuous hydrodynamic stress on the top of cilia.

Planar cell polarity of the tissue

We performed immunofluorescent staining against the core PCP protein Vangl1 to evaluate whether the cultures were polarized and to measure the main directions of planar cell polarity.

We observed large domains where Vangl1 was aligned (fig. 3a-c, typical scale 200 µm) or circularly organized (Fig. 3d-f, typical scale 100 µm). These domains likely correspond to areas of linear mucus transport towards the culture edges and to local mucus swirls respectively. At a millimetric scale, a circular order of Vangl1 is also observed under an extended millimetric mucus swirl, as illustrated in Fig. 3g.

In summary, the characteristic size of mucus transport estimated by the size of the spontaneously emerging mucus swirls increases strongly with cilia density to reach the entire culture size at physiological cilia density. These swirls are associated with both a circular order of ciliary-beat directions and a circular planar cell polarity of the tissue below the swirls.

To understand the specific role of the hydrodynamic forces exerted by mucus on the coordination of ciliary beat directions, we removed and added mucus on mature epithelial surfaces.

Mucus drives the loss and recovery of the ciliary orientational order

We removed rotating mucus on the surface of mature cultures by an extensive apical washing.

Such washing empties the goblet cells thus limiting the production of new mucus for a week.

A week after mucus removal, the giant ciliary-beat orientational order spanning the culture chamber was lost and several local circular patterns emerged instead (Fig. 4a). Many topological defects appeared at a high density of about 6.10 2 to 8.10 2 defects/mm². S decreases to 0.1 over a length scale of L= 300 µm compared to S=0.4 before mucus removal (Fig. 4b). At the tissue level, immunofluorescent staining of the core PCP protein Vangl1 clearly shows that no general circular PCP pattern (similar to the one shown in Figure 3g) is detectable in absence of mucus. Instead, we could observe multiple local PCP swirly patterns reminiscent of the local ciliary swirls reported in this condition (Extended Data Fig. 3).

Strikingly, a similar loss of ciliary-beat orientational order was observed when mucus got dehydrated on top of the epithelium and was no more transported although cilia still beat below (Fig. 4b and movie S9). This shows that mucus transport is mandatory to maintain the macroscopic ciliary-beat orientational order.

Next, the mucus stored at the edge of the culture chamber was diluted and redistributed over the entire epithelial surface. After two days, a new millimetric mucus swirl was reconstituted at the epithelial surface (see movie S10), associated with a strong circular orientational order of underlying ciliary beats (Fig. 4c and Supplementary Fig. 1). The mucus flow forces ciliary beats to recover a large-scale orientational order. Over 24h we observed re-orientations of beat directions up to 30-35°. These results show that the direction of ciliary beats is not frozen after ciliogenesis in in-vitro reconstituted bronchial epithelium and supports the idea that cilia actively respond to the hydrodynamic force exerted by the mucus.

To understand the underlying hydrodynamic mechanisms, we developed a simple hydrodynamic 2D-model.

The hydrodynamic coupling between mucus and cilia drives the emergence of a global order

The modeling approach is based on the schematic view shown in Figure 5a. The epithelial surface is described using hexagonal unit elements of size D, which represent the elementary patch of neighbouring ciliated cells that have a common direction of ciliary beat, next to nonciliated cells (club and goblet cells).The fluid in contact with the epithelium is decomposed in two layers [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF] : the low-viscosity periciliary layer (PCL), with a thickness similar to cilia length and where ciliary beating mostly occurs, and the high-viscosity mucus layer, which contacts cilia tips. We focus on long-range flow phenomena developing in directions parallel to the epithelium, over length scales much larger than the typical fluid layer thickness. Therefore, the mucus flow is considered to be two-dimensional. In this configuration, the PCL behaves as a slip layer. This approach, illustrated by the schematized velocity profiles in Figure 5a, is consistent with prior results issued from three-dimensional numerical simulations [START_REF] Ding | Mixing and transport by ciliary carpets: a numerical study[END_REF][START_REF] Chateau | Transport efficiency of metachronal waves in 3D cilium arrays immersed in a two-phase flow[END_REF] .

In this model, ciliary activity is averaged over time and space. We consider only the global driving force that propels the mucus in a specific direction by averaging over time the force resulting from the stroke recovery cycle of beats. At the scale of one hexagonal element, we average the driving forces exerted by all the cilia of the considered element.

The flow is modeled by the incompressible two-dimensional Navier-Stokes equations and simulated using a lattice-Boltzmann method. The effect of the epithelial wall on the flow is modeled through a force density (i.e. momentum source/sink) fn (in N. m -3 ) exerted on the mucus. Two types of momentum sources are considered. The friction force relates to the normal shear stress in the PCL. It applies over the entire epithelial plane and is assumed to linearly depend on the mucus velocity U, fv = -κU, where κ is a PCL friction parameter. From a dimensional analysis, the parameter κ writes as κ~μp/δp², where μp and δp are the viscosity and the thickness of the PCL, respectively.The ciliary force f c , locally driving the mucus in a specific direction as a result of time-averaged ciliary beating, is applied on ciliated elements only. Its magnitude is constant and its orientation is forced to gradually align with the local mucus velocity, as detailed in Materials and methods. This streamwise-alignment rule provides a hydrodynamic interaction mechanism between ciliated elements. The total momentum source in the flow equations is expressed as

f n = f v + f c
The model depends on three non-dimensional physical parameters, φ, the Reynolds number Re = ρD||fc||/κμ, with ρ and μ the density and viscosity of the mucus, and a non-dimensional hydrodynamic interaction length λ = (μ/κ) 1/2 /D, that is also written λ~(μ/μp) ½ (δp/D). The latter parameter represents the typical range of influence of a ciliated element and defines the maximum interaction length between ciliated cells. Here, the flow solution is computed in a square domain of ∼7000 elements by choosing parameter values in physiological ranges. φ ∊ [0.1, 0.7], D ≈2 10 -5 m, δp≈10 -5 m and μp ≈10 -3 Pa.s are estimated from the ciliogenesis experiment and in Button et al. [START_REF] Button | A periciliary brush promotes the lung health by separating the mucus layer from airway epithelia[END_REF] ., respectively. The literature reports widely scattered measurements of mucus viscosity [START_REF] Lai | Micro-and macrorheology of mucus[END_REF] . We have varied it in the range 5. 10 -3 -5 10 -2 Pa.s as recently found for healthy mucus [START_REF] Jory | Mucus Microrheology Measured on Human Bronchial Epithelium Culture[END_REF] . Therefore, λ is in the range [1, 4]. Finally, Re = 0.1 prevents inertial effects. The initial location of ciliated elements and their initial orientation are randomly set. A free-slip boundary condition is used on the limits of the square domain in order to reproduce a confinement similar to that of the experiments. We obtain a phase diagram which captures the observed mucus flow patterns (Fig. 5b). Indeed, mucus flows exhibit swirly patterns, whose typical length scale increases as a function of cilia density as experimentally observed during ciliogenesis. For λ = 2, when φ ≥ 0.5 a global swirl spanning the entire chamber emerges, similar to what we observe for cilia densities ≥ 50 % in presence of mucus. In this long-range organization regime, the chamber walls enforce a global circulation of the mucus flow. For φ ~ 0.3 simulations predict the emergence of a couple of large swirls whose size ranges between 600-1000 µm in diameter, in the same order with what we observe experimentally with swirls of mean diameter of 1.25 +/-0.2 mm. For densities lower than 0.3, simulations predict a myriad of local swirls whose size ranges between 300-600 µm in diameter, in agreement with the observed swirls with mean size of 500 µm. At a fixed cilia density of 70 %, the decrease of the interaction length λ from 2 to 1, equivalent to a decrease of the viscosity of the transported fluid, drives the transition from a global swirl to multiple local swirls. This λ dependence of the loss of the global orientational order is reminiscent of the observed transition which occurs when we remove the mucus and replace it with culture medium (Fig. 5c). The combination of the experiment with the hydrodynamical model provides a strong evidence that long range hydrodynamic interactions are sufficient to establish a global orientational order in the bronchial epithelium.

Discussion

In this work we first described the fascinating spontaneous emergence of mucus swirls on reconstituted bronchial epithelia. We quantitatively measured the formation and growth of these swirls during ciliogenesis as well as the final long-range mucus transport observed for a cilia density higher than 50%. We highlight that the local milling patterns observed during ciliogenesis are not determined by the circular geometry of the culture chamber because they are numerous and randomly spatially distributed. However, the giant swirl observed at the end of ciliogenesis 'feels' the whole culture chamber and its circular shape is determined by the chamber geometry. Associated with these swirls, we showed that the ciliary beats and the planar cell polarity of the tissue have a circular orientational order. As already reported in previous studies [START_REF] Feriani | Assessing the Collective Dynamics of Motile Cilia in Cultures of Human Airway Cells by Multiscale DDM[END_REF][START_REF] Sanderson | Ciliary activity of cultured rabbit tracheal epithelium: beat pattern and metachrony[END_REF] , metachronal waves do not develop on long distances. This is probably due to the beat frequency scattering that prevents long-range synchronization and to the fact that ciliated cells do not form a dense carpet since they are separated from each other by goblet cells. On this point, this system is very different from swimming microorganisms (paramecium, volvox...) for which ciliary beats synchronize and metachronal waves develop [START_REF] Brumley | Flagellar synchronization through direct hydrodynamic interactions[END_REF][START_REF] Wan | Reorganization of complex ciliary flows around regenerating Stentor coeruleus[END_REF][START_REF] Gueron | Cilia internal mechanism and metachronal coordination as the result of hydrodynamical coupling[END_REF] . In these systems, metachronal waves optimize the velocity of the propelled fluid which might be crucial to search for food or escape predators. In bronchial epithelium, the absence of metachronal waves on long distances suggests that the transport velocity of the mucus is not a critical parameter in healthy situation. Yet, what is important in this system is the long range coordination of the ciliary beat directions which relies on the hydrodynamics active coupling between mucus and cilia, constrained by the planar cell polarity of the tissue.

We proposed a model that describes the mechanism of mucus swirl formation and coordination of ciliary-beat directions based on our observations that mucus is essential to maintain and recover large-scale transport. The model relies on the active coupled response of cilia whose beat direction locally tends to align with the mucus flow resulting from all surrounding ciliary beats. It predicts the emergence of mucus swirls through a circular orientation of ciliary beats, their growth with cilia density and with the length of hydrodynamic interactions, directly related to mucus viscosity. Noteworthy, the ciliary beat system belongs to the class of active nematic living materials, recently observed in several living organisms [START_REF] Kawaguchi | Topological defects control collective dynamics in neural progenitor cell cultures[END_REF][START_REF] Saw | Topological defects in epithelia govern cell death and extrusion[END_REF][START_REF] Duclos | Topological defects in confined populations of spindle-shaped cells[END_REF] . Here, the bronchial epithelium is an even more complex material composed of a double nematic system. Indeed, cilia are not only subjected to the mucus-driven hydrodynamic force, but are also constrained at their base by the PCP protein long-range organization. While it cannot be excluded that an intrinsic self-organization of PCP proteins may force the coordination of ciliary beat directions, our observations rather suggest that the organization of PCP proteins is the result of hydrodynamic forces. On a local scale however, one could put forward that before tissue confluence, the intrinsic turbulent dynamics of jamming cells could generate small circular structures [START_REF] Blanch-Mercader | Turbulent Dynamics of Epithelial Cell Cultures[END_REF] that could facilitate the emergence of local swirls. At large scale, in-vivo, internal forces exerted on the epithelium during embryogenesis are involved to set the planar polarization [START_REF] Butler | Planar cell polarity in development and disease[END_REF] . In in-vitro cultures, these forces do not exist and it is the hydrodynamic forces that could lead to the large-scale self-organization of PCP proteins. This would explain the circular organization of PCP and the plasticity of mature cultures allowing the reorganization of the directions of ciliary beats under the effect of hydrodynamic forces.

We conclude with two remarks, First, our model extends the class of Vicsek-like models [START_REF] Vicsek | Novel Type of Phase Transition in a System of Self-Driven Particles[END_REF][START_REF] Costanzo | Spontaneous emergence of milling (vortex state) in a Vicsek-like model[END_REF] based on simple rules between moving individuals and used to describe the collective behavior observed in active biological matter, such as bacterial colonies, cytoskeletal elements and epithelial monolayers [START_REF] Schaller | Polar patterns of driven filaments[END_REF][START_REF] Sokolov | Concentration dependence of the collective dynamics of swimming bacteria[END_REF][START_REF] Sumino | Large-scale vortex lattice emerging from collectively moving microtubules[END_REF] . Second, from an applied perspective, our results could provide important insights for adapting gene therapy strategies, for example to determine the number of ciliated cells that need to be repaired in primary ciliary dyskinesia.

Materials and methods

In-vitro bronchial epithelium -Cultures of human bronchial epithelium reconstituted from primary cells, were either bought from Epithelix (MucilAir) or cultured from human transplant donor lungs deemed unsuitable for transplantation and donated to medical research. The ethics committees of the institutions involved approved this study (CERC-SFCTCV-2018-5-6-9-8-32-DjXa). Primary human bronchial epithelial cells were isolated by protease digestion of human airways, and cells were cultivated under Air-Liquid Interface (ALI) conditions, as previously described [START_REF] Gras | An ex vivo model of severe asthma using reconstituted human bronchial epithelium[END_REF][START_REF] Vachier | Glucocorticoid Receptors in Bronchial Epithelial Cells in Asthma[END_REF] . All experiments were performed in accordance with relevant guidelines and regulations.

The tissue is reconstituted at the Air Liquid Interface in transwells of 6 mm in diameter. The apical side of the tissue is at the air interface and the basal side is in contact via a porous membrane (pore size of 0.4 µm) with the culture medium: Epithelix differentiation culture medium during ciliogenesis, then, Epithelix MucilAir culture medium; Pneumacult ALI medium (Stem cell) for homemade cultures. The culture medium (700 μl) is replaced every two days.

To avoid excessive dehydration of mucus overtime, which greatly impacts the mucus transport velocity (Supplementary Fig. 2), we add 3µl of culture medium on the apical side twice a week.

Removal of mucus

For experiments in absence of mucus, we remove the mucus by performing an apical washing.

We add 200 μl of culture medium on the apical side for 10 minutes, followed by three rinses with the culture medium. After an apical washing, a thin layer of culture medium (∼ 10 -40 μm) remains on top of the cilia. Eventually, if patches of mucus remain stuck on the surface of the epithelium we add 10 µl of culture medium on the apical side that we leave overnight.

Imaging -We imaged the cultures in bright field on an inverted Nikon Eclipse Ti microscope with a x20 objective and a Luminera infinity camera. The temperature is maintained at 37 °C and a humidified air flow with 5% CO2 is applied to maintain a physiological pH. To image the cultures at the same locations on different days as well as to image the samples on large areas (5x5 tiles), we used a motorized stage calibrated in the xy plane and saved the different positions. First, we visualized the streamlines of the mucus flows by following cellular debris trapped in the mucus to characterise the emergence and growth of mucus flow patterns. Then, after slightly diluting the mucus, by adding of 3-4 µl of culture medium on the apical side to get rid of the largest debris, we imaged the underneath beating cilia in the same area. To quantify the ciliogenesis we imaged on the central part of each sample 25 tiles resulting in an area of about 2mm x 1.5mm. For imaging planar polarity of Vangl1, samples were mounted between a glass slide and a coverslip embedded on a drop of ProLong mounting media. Images were taken with the confocal microscope ZEISS LSM 880 reverse standing AxioObserver 7. The objective used were Zeiss 20x Plan Apochromat 20x/0.8. Two-colors confocal z-series images were acquired using sequential laser excitation, and a z-slices interval 1µm, using Tile Scan mode covering areas of ~2x2 mm². Images were converted into single plane projection (MIP) and analyzed or edited using ImageJ software.

Immunofluorescence staining of Vangl1 and ZO1 -

Samples were first fixed in 100% methanol, 10 minutes at -20°C. Then, washed 3x3min in PBS at 4°C and blocked with BSA 3% in PBS 30min at RT. Primary antibodies against Z01 (mouse monoclonal IgG1 Thermo Scientifique 33-9100) and Vangl1 (rabbit Sigma Aldrich HPA025235) were diluted at 1/200 and 1/500 respectively in blocking medium and incubated 45min at RT. After 2x5min washes with PBS-Triton 0.1%, secondary antibodies: Alexa mouse IgG1 568 (ThermoFisher A21124) diluted at 1/800 and Alexa Fluor anti-rabbit 488 (Life Technologies A21206) diluted at 1/500 were incubated in blocking medium 30min at RT. Image processing and quantification -Cilia density 𝜑𝜑, beating directions of cilia and localisation of topological defects were quantified using in-house image processing programs implemented in python and applied on videomicroscopy acquisitions (40 fps). We first performed a standard deviation projection over 40 frames to reveal the trajectory of the tip of the cilia. We binarised this image by applying a threshold to compute 𝜑𝜑. We implemented the algorithm described in Püspöki et al. [START_REF] Püspöki | Transforms and Operators for Directional Bioimage Analysis: A Survey[END_REF] based on the structure tensor computation to determine the direction of each pixel on previously obtained standard deviation projection images. The director field of ciliary beats was computed by averaging the orientation of individual pixels in a box of 15x15 pixels, corresponding to the size of a single ciliated cell. Then we localized the topological defects from the director field by computing the winding number [START_REF] Huterer | Distribution of singularities in the cosmic microwave background polarization[END_REF] , i.e. we compute how much the director field rotates around a considered cell (see Supplementary Fig. 3 for a detailed explanation).

To visualize the streamlines of the mucus flow on a still image, we perform a temporal standard deviation projection from video acquisition on 100 to 300 frames depending on the frame rate.

Hydrodynamic model -In our model, the two-dimensional mucus flow is driven by two force densities, denoted by fc and fv. The force fv results from the shear stress in the PCL; it linearly depends on the mucus velocity U, namely fv = -κ U, and is applied everywhere. The ciliary force fc is only applied on ciliated regions. Its orientation is defined on each ciliated element through the angle 𝜃𝜃c i (t), which is the angle between the ciliary force and a given direction on the i-th ciliated element at time t. The angle difference between the local mucus flow and the forcing direction Δ𝜃𝜃 i =𝜃𝜃f i -𝜃𝜃c i , where 𝜃𝜃f i (t) is the mucus velocity averaged over the i-th ciliated cell.

The orientation of the ciliary forcing is updated following a streamwise-alignment rule, namely S is computed in a sliding window of size L with L varying from Lmin= 6 µm, the size of a ciliated cell to Lmax = 170 µm, the size of the swirl. The window is moved over the entire surface of the image by L/3 steps. The comparison of the nematic order of the cilia, blue curve, with the nematic order computed for a perfect circular orientation, red curve, shows that underneath a swirl of mucus, the cilia exhibit a strong circular order. d, Dynamics of ordering of cilia in a sample exhibiting a single swirl spanning the entire culture. The parameter order S is measured at mid-distance (~1.5 mm) between the center of rotation and the border of the chamber and is averaged on a total area of ~1mm². Day 1 defines the first day the swirl spans the culture. At a length scale of 50 µm we observe a sharp increase of order within 24h, this indicates that cilia get more aligned overtime. At a larger length scale of 200 µm there is a 24h delay time before we observe the same increase of S. After four days, S reaches a plateau at ~0.5. e, Temporal evolution of the density of topological defects in the same experimental conditions than in d. Underneath the rotation of mucus, the number of defect decreases within a couple of days and then remains constant. The 2-dimensional epithelium is paved with ciliated elements and non-ciliated elements of characteristic size D. The PCL is a shear region with the fluid velocity in the xy plane equal to zero on the epithelial surface and is equal to the mucus velocity at the PCL/mucus interface. The mucus is transported in the xy plane. b, Phase diagram of the mucus flow patterns in the λ, φ space. Panels show streamlines and iso-contours of the non-dimensional flow velocity magnitude, normalized in each case by the maximum velocity in the domain and ranging from 0 (off-white) to 1 (dark blue). In all cases, the solution has reached a steady state, and the motile force on motile elements is aligned with the streamlines. At a given λ, increasing φ from 0.1 to 0.7 results in a transition from a local and small swirly patterns to a large global swirl. At a given φ, an increase of the hydrodynamic interaction length results in an increase of the characteristic size of the mucus flow patterns. Red rectangle indicates the part of the computational domain examined in c. Scale bar = 300 µm, and all panels have the same size. c, Comparison between simulation and experiments. At a fixed cilia density of 0.7, typical of a mature culture, increasing λ from 1 to 2 results in a transition from multiple swirls to a single global swirl. In-vitro we observe the same mucus flow patterns on mature cultures in absence or in presence or mucus, respectively. Scale bars = 300 µm.

Extended Data

Finally, samples are

  washed 2x5min. The membrane of the Transwell was cut and mounted in the coverslip. Quantification of Vangl1 orientation -We measured manually the orientation of Vangl1 protein on each individual cell in the regions of interest using imageJ software. Then we plotted the distributions of beat directions according to the measured Vangl1 orientations with beat direction = Vangl1 orientation + 90°.

Fig. 1 |

 1 Fig. 1 | Emergence of mucus flow patterns during ciliogenesis. a, Schematic of the ALI chamber: the basal side of the bronchial epithelium is in contact with the culture medium through a porous membrane and the apical side points towards the air. The epithelium has a pseudo stratified structure made of basal cells, goblet cells, which produce the mucus and ciliated cells. During the forward stroke, cilia tips poke in the mucus layer while their recovery motion takes place in the periciliary layer. b, Dynamics of the ciliogenesis. The curve represents the temporal evolution of the percentage of surface covered by active cilia after initiation of the differentiation. Data are from samples #1 to #10 from 5 patients, the color represents a patient and different symbols of the same colors represent different cultures from the same patient. Each data point is an average over at least an area of 1mm x 1mm and error bars are the standard deviation. Characteristic length scale of mucus transport associated to the cilia density are indicated in green, the values represent the average diameter of the swirls +/-standard deviation. c-d, Images representative of two cilia densities, 𝜑𝜑 = 30% (day 14) in c and 𝜑𝜑 = 80% (day 22) in d. Images are the result of a standard deviation projection over 40 frames acquired at 40 fps. Scale bars = 30 µm. e-g, Formation of a swirl of mucus from day 12 to day 14 in sample #13. The three images are the same field of view. Images are the result of a standard deviation projection over 80 frames at 0.2 fps in e (movie S2), 80 frames at 0.2 fps in f (movie S3) and 100 frames at 1 fps in g (movie S4), respectively. Scale bars = 50 µm. h-i, Growth of a mucus swirl between day 12 and day 15, in sample #8, the diameter of the mucus swirl increases by successive accumulation of mucus transported in the vicinity of the swirl. Images results from the standard deviation projection over 100 frames at 0.2 fps in h and over 150 frames at 40 fps in i. Scale bars = 100 µm.

Fig. 2 |

 2 Fig. 2 | Orientational order of ciliary beats. a, On the background image the white strokes reveals the trajectories of the tip of bundles of cilia. The director field of the beating directions is plotted as an overlay. Red circles and blue triangles are located where a discontinuity arises in the director field and are +½ and -½ nematic topological defects respectively. A schematic of +½ and -½ defect is drawn in inset. The beating directions of cilia are revealed by performing a standard deviation projection over 40 frames at 40 fps. The director field is computed by averaging the orientation in boxes of 15x15 pixels, corresponding to the size of a single ciliated cell. Day of observation #20. Scale bar = 30 µm. b, Beating directions of cilia underneath a local swirl of mucus, sample #11. The colors code for the direction of beating according to the orientation scale in the upper right corner. The director field averaged over 15x15 pixels reveals a robust circular orientational order. Only few topological defects are present. Scale bar = 20 µm. c, Nematic order S, blue curve, associated to the cilia beating orientations of the panel b. S is computed in a sliding window of size L with L varying from Lmin= 6 µm, the size of a ciliated cell to

Fig. 3 |

 3 Fig. 3 | Planar cell polarity organisation. a & b, Long range linear polarization of the PCP Vangl1 protein (blue-green label). Tight junctions, ZO1, are labelled in red. Scale bars = 20 µm. In insets, a zoom on a few cells is shown. Scale bars = 5 µm. c, Quantification of the direction of ciliary beats according to the orientation of Vangl1: the beat direction is orthogonal to the Vangl1 edge. The angle distribution of beat directions corresponds to the field in a & b, the red line indicates the median direction of 146° and the standard deviation is 34°. d & e, Vangl1 exhibits a circular pattern reminiscent of the observed local swirls of ciliary beats. The yellow lines are a guide for the eyes and represent the main orientations orthogonal to vangl1. Scale bars = 20 µm. f, Angle distribution of beat directions, computed from Vangl1 orientation, corresponding to the field in d & e, angles are well distributed over the 0-180° and the red line indicates the median direction of 99°. A perfect circular order would lead to a median of 90°. g, PCP protein Vangl1 over an area of ~1.5 mm x 1 mm in a sample where a millimetric mucus

Fig. 4 |

 4 Fig. 4 | The transport of mucus drives the self-organization of the epithelium. a, Size distribution of the local swirls that form one week after mucus removal. Data are pooled from samples #6, 10 and 11. b, Comparison of the nematic order computed at different length scales. One week after the mucus has been removed or became immobile due to dehydration, we observe a loss of order at length scales greater than ~25 µm. c, Orientation map of the beating direction of cilia 3 days after addition of a diluted mucus on top of a disorganized culture (sample #15). The colors code for the beating direction according to the orientation scale in the inset. The director field is plotted by averaging the directions over a group of 10 cells. The bronchial epithelium exhibits a circular orientation of the ciliary beats over a millimetric distance. Scale bar = 200 µm.

Fig. 5

 5 Fig.5Hydrodynamic model of the bronchial epithelium. a, Schematic of the modelled epithelium. The 2-dimensional epithelium is paved with ciliated elements and non-ciliated elements of characteristic size D. The PCL is a shear region with the fluid velocity in the xy plane equal to zero on the epithelial surface and is equal to the mucus velocity at the PCL/mucus interface. The mucus is transported in the xy plane. b, Phase diagram of the mucus flow patterns in the λ, φ space. Panels show streamlines and iso-contours of the non-dimensional flow velocity magnitude, normalized in each case by the maximum velocity in the domain and ranging from 0 (off-white) to 1 (dark blue). In all cases, the solution has reached a steady state, and the motile force on motile elements is aligned with the streamlines. At a given λ, increasing φ from 0.1 to 0.7 results in a transition from a local and small swirly patterns to a large global swirl. At a given φ, an increase of the hydrodynamic interaction length results in an increase of the characteristic size of the mucus flow patterns. Red rectangle indicates the part of the computational domain examined in c. Scale bar = 300 µm, and all panels have the same size. c, Comparison between simulation and experiments. At a fixed cilia density of 0.7, typical of a mature culture, increasing λ from 1 to 2 results in a transition from multiple swirls to a single global swirl. In-vitro we observe the same mucus flow patterns on mature cultures in absence or in presence or mucus, respectively. Scale bars = 300 µm.

Fig. 1 |

 1 Spatial distribution of ciliated clusters during ciliogenesis. a, Size distributions of ciliated cluster for different cilia densities φ. The size of clusters is expressed as the number of ciliated cells with a mean cell area of 40 μm2. At the beginning of ciliogenesis, most clusters are formed by 1-2 ciliated cells. Upon differentiation of new ciliated cells, more clusters are formed, which grow larger and eventually connect. We observe the appearance of a percolating cluster (for a density above 60 %), in coexistence with clusters of finite size. The size distribution for φ = 68 corresponds to the clusters of finite size, and the percolating cluster has disappeared. For each cilia density, φ and the cluster size distributions are computed on 9 tiles of size 440 μm x 330 μm. φ is the mean density + /-standard deviation. Data are from samples #5 and #6. b-f, Characteristic spatial distribution of ciliated cells at different φ. For φ = 4 %, only sparse clusters of few cells are present (b). While φ increases, the number and the size of clusters increase (c-e). For φ = 70 % we observe a percolating cluster, in coexistence with clusters of finite size (f). The colorbar indicates the number of ciliated cells per cluster and the largest cluster appears in yellow. Scale bars are 50 μm. Extended Data Fig. 2 | Cilia beat frequency (CBF). a, Frequency map on an area where a swirl is present. The colour codes for the CBF and the CBF distribution is shown in b. Black regions correspond to locations where flowing opaque debris in the mucus prevent the visualization of beating cilia underneath. Scale bar = 40 μm. b, CBF distribution, CBF = 7.6 + /-1.3 Hz. A lower cut-off has been defined at 4 Hz to prevent artefacts due to flowing debris. c, Kymograph plotted along the blue path displayed in a, corresponding to a ciliated area with the same CBF. Patterns on the kymograph allows to identify ciliary beats with the same phase. This is highlighted by the green lines at the beginning. Yet, the synchronization is not stable over time and the loss of synchronization is highlighted by the blue lines. Spatial scale bar is 20 μm and temporal scale bar is 0.1 s. Extended Data Fig. 3 | Planar cell polarity patterns in absence of mucus. PCP protein Vangl1 over an area of ~1.3 mm x 1.1 mm in a sample (#21), where a millimetric swirl was initially present and then the mucus had been removed by performing successive washing. The sample has been fixed after one week in absence of mucus. The yellow lines are a guide for the eyes and represent the main orientations orthogonal to Vangl1 (featuring cilia beat directions). The millimetric circular order is lost and PCP exhibits instead numerous swirly patterns. Angle distributions of ciliary beats, computed from Vangl1 orientation, are plotted for 4 quadrants as in Fig. 3 (dashed triangles, arc of a circle with radius of 500 μm and subtending an angle of 35°). The red lines indicate the median of the distributions with α = 90 for top left, α = 96 for top right, α = 96 for bottom right and α = 83 for bottom left. Scale bar = 100 μm.

Acknowledgements

The authors thanks Dr. Pascal Thomas and Dr. Xavier Benoit D'Journo, service de chirurgie thoracique, Hôpital NORD, AP-HM, for providing biological material. The MACBION project leading to this publication has received funding from Excellence Initiative of Aix-Marseille University -A*MIDEX, a French Investissements d'Avenir programme, and from the French National Research Agency (ANR) (SINUMER project ANR-18-CE45-0009-01). The research leading to these results has received funding from the People Programme (Marie Curie Actions) of the European Union's Seventh Framework Programme (FP7/2007-2013) under REA grant agreement n. PCOFUND-GA-2013-609102, through the PRESTIGE programme coordinated by Campus France. A.N. was supported by a fellowship from Fondation pour la Recherche Médicale.

Data availability

All data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.

Code availability

Custom codes that were used to analyze experimental data within this manuscript are available from the corresponding author upon reasonable request.

where Δt is the model time step, 𝛺𝛺 is a fixed angular velocity and θ0 is an angle threshold. The angle θ0, which is kept to very small values in this study, is introduced in order to allow the emergence of steady solutions. More details on the initialization of the model, definition of the parameters and numerical method are present in the supplementary information.

Author contributions E.L and A.V designed the study. S.G, U.D and J.F conceived the computational model. E.L designed and carried out image analysis. E.L and C.J. acquired experimental data. A.N. acquired the data on PCP. E.L, CJ and A.N analyzed experimental data. S.G. implemented and performed simulations. D.G. performed bronchial epithelium cultures. All authors discussed the results and the manuscript. E.L, S.G and A.V wrote the paper.

Competing interests

The authors declare no competing interests


Damien Bigourd

Philippe Morin

Jérôme Dubertrand

Coralie Fourcade-Dutin

Hervé Maillotte

Yves Quiquempois

Géraud Bouwmans

Parametric gain shaping from a structured pump pulse



Abstract-We investigate optical parametric amplification in a fiber pumped by a shaped pulse. The temporal structure of the pump impacts directly the instantaneous spectral gain and its temporal measurement is achieved in the spectral domain thanks to a chirped signal with a broad spectrum.

Index Terms-Four wave mixing, Parametric amplifier, fiber optics I. INTRODUCTION Fiber optical parametric amplifier (FOPA) has been widely investigated in the context of telecommunication [START_REF] Hansryd | Fiber-based optical parametric amplifiers and their applications[END_REF] for its interesting features such as the instantaneous amplification with large bandwidths [START_REF] Marhic | Broadband fiber optical parametric amplifiers[END_REF][START_REF] Ho | 200nm-bandwidth fiber optical amplifier combining parametric and Raman gain[END_REF], high gain values [START_REF] Torounidis | Fiber-optical parametric amplifier with 70-dB gain[END_REF] and low noise figures [START_REF] Voss | Measurement of the photon statistics and the noise figure of a fiber-optic parametric amplifier[END_REF]. In most cases, the amplification is achieved with continuous waves (CWs) involving a pump with a constant peak power, an amplified signal and a generated idler. In order to reach higher peak power, FOPA is preferably pumped in ps or ns pulsed regime and we often assume that the pump power is constant during the pulse duration [START_REF] Bigourd | High-gain, optical-parametric, chirped-pulse amplification of femtosecond pulses at 1 µm[END_REF][START_REF] Fourcade-Dutin | 12 THz flat gain fiber optical parametric amplifiers with dispersion varying fibers[END_REF][START_REF] Bigourd | Ultra-broadband fiber optical parametric amplifier pumped by chirped pulses[END_REF]. However, the temporal shape of the pump pulse plays a crucial role in the amplification process and it can be of great use to target a specific signal profile [START_REF] Finot | Influence of the pump shape on the modulation instability process induced in a dispersion-oscillating fiber[END_REF]. Indeed, the instantaneous parametric gain and its spectral shape depends strongly on the power. Therefore, it also has a temporal distribution since a specific gain shape is expected at each time of the pump profile. Various systems take benefit from this behavior. For example, a train of ps pulses can be obtained at high repetition rate when a sinusoidal pump is coupled in a FOPA with a CW signal [START_REF] Vedadi | Theoretical study of high repetition rate short pulse generation with fiber optical parametric amplification[END_REF]. It is also well known from several investigations achieved in nonlinear crystals, that the features (as efficiency, bandwidth, signal to noise ratio) of the parametric amplifiers can be optimized by adapting the pump waveform [START_REF] Batysta | Spectral pulse shaping of a 5 Hz, multi-joule, broadband optical parametric chirped pulse amplification frontend for 10 PW laser system[END_REF] or the durations of the pump and the signal [START_REF] Moses | Temporal optimization of ultrabroadband high-energy OPCPA[END_REF]. In the manuscript, we investigate the role of the pump shape in a FOPA to tune the instantaneous parametric gain. However, the time distribution measurement of the spectral gain is not straightforward. When the parametric fluorescence (which is representative of the spectral gain) is recorded with a spectrograph, the gain temporal structure cannot be directly inferred since the measurement is a time integrated (averaged) spectrum. Therefore, in our case, the structured gain is measured thanks to a chirped signal where the instantaneous frequency evolves quasi-linearly with time. The temporal gain shape influences directly the signal in the temporal and spectral domains. This scheme should be of prime interest to tune a pulse in any experiment requiring a chirped signal as in lidar [START_REF] Delfyett | Chirped pulse laser sources and applications[END_REF], optical time-stretch imaging [START_REF] Wong | Asymmetric-detection time stretch optical microscopy (ATOM) for ultrafast high-contrast cellular imaging on flow[END_REF] or laser amplifier [START_REF] Mecseki | Flat-top picosecond pulses generated by chirped spectral modulation from a Nd:YLF regenerative amplifier for pumping few-cycle optical parametric amplifiers[END_REF]. The letter is organized as follows. Section II explains the principle of the temporal gain shaping from a structured pump pulse through numerical simulations. Section III presents the experimental set-up and some examples of shaped gains. Finally, Section IV provides concluding remarks.

II. PRINCIPLE

The general principle relies on the instantaneous variation of the gain due to the temporal profile of the pump. The main focus is to generate parametric amplification using one strong pump at a single wavelength with a temporal variation. In order to explain how the spectral gain changes with the pump profile, we performed numerical simulation by integrating the nonlinear Schrödinger equation along the propagation z

2 3 4 2 2 3 4 2 3 4 1 2! 3! 4! E i E E i E EE z                    (1)
describing the evolution of the slowly varying total electric field E in an optical fiber. Equation ( 1) is solved with the standard split-step Fourier method including a pump pulse with various shapes, a signal and a weak random initial condition that mimics the quantum fluctuation. The temporal resolution is 20 fs and the number of point is 2 18 . The experimental parameters, as the pump and fiber characteristics, is included [START_REF] Morin | µJ-level Raman-assisted fiber optical fiber parametric chirped pulse amplification[END_REF] to illustrate qualitatively some experimental behaviors. 20, 30 and 40 are the second, third and fourth order dispersion terms at the pump angular frequency and equal to -6.1x10 -28 s 2 /m,30=4.7 x10 -39 s 3 /m, 40= -7.37x10 -53 s 4 /m, respectively. The nonlinear coefficient is =1. Bouwmans, E. Hugonnot performed for the fiber length used during the experiment (L=13 m) and for a shorter one (L=6 m) to evaluate the role of the saturation. The spectral gains as a function of time is simulated for several pump shapes and in all cases, the energy equals to 1 µJ. With such energy level, some of us showed that the Raman process has a weak contribution compared to the parametric effect [START_REF] Morin | µJ-level Raman-assisted fiber optical fiber parametric chirped pulse amplification[END_REF]. As this observation was also confirmed by our preliminary simulations, the Raman contribution is neglected in the following simulations. For all pump shapes, the gain is measured by injecting a quasi-CW signal with a top-hat profile and an energy of 6.2 pJ. The frequency offset from the pump Ω/(2 is tuned from 0.5 THz to 4.5 THz. At the output of the fiber, the amplified signal is selected with a tunable bandpass filter with a bandwidth equals to 4 GHz (FWHM) and Fourier transformed. A pump pulse is firstly injected with a tophat shape (Fig. 1.a-blue line) with a duration of 1.4 ns at Full Width Half Maximum (FWHM) and a peak power equals to 717 W. In the case of a top-hat shape (Fig. 1.a blue line) and for L=6 m, the gain profile is identical for all the time where the peak power is constant (Fig. shapes. The gain at 3.6 THz is not maximum as for L=6 m due to the pump depletion and the temporal evolution of the phasematching that varies with the power. The gain peaks also shift toward the pump frequency (Ω/(2 =0 THz) due to the saturation [START_REF] Inoue | Signal wavelength dependence of gain saturation in a fiber optical parametric amplifier[END_REF]. As expected, the depletion and the peak shift is time dependent for the Gaussian and asymmetric power profiles (Fig. 1.e,g). The temporal shape of the amplified signal is also affected by the gain temporal profile at their frequencies. For example, when a CW signal is injected at 3.6 THz together with the shaped pump (Fig. 1.b-g-white dashed lines), the profile of the amplified signal is significantly different in the four cases (Fig. 2.a-solid lines). In addition, the temporal variation of the amplified signal is different when a chirped signal is injected together with the shaped pump since the signal has a timefrequency distribution. Fig. 2.a (dashed lines) shows also the examples when the chirped signal has a spectrum centered at 3.6 THz and a chirped rate equals to 1.2 THz/ns (Fig. 1.b-ggreen lines). The gain shape depends strongly on the pump profile. In these cases, the horizontal abscissa represent both the temporal and spectral domains. For a flat-top profile, the temporal gain has a low order super-Gaussian shape (Fig. 2 blue dashed line) due to the unmatched parametric process at Ω/(2≠3.6 THz (Time≠0 ps). In the three other cases, the gain is structured since it depends both on the temporal evolution of the pump and the time-frequency distribution of the signal. For longer fiber (Fig. 2. b), the gain profiles becomes well different from Fig. 2.a due to the saturation effect and the time dependence of the phase matching.

III. EXPERIMENTAL IMPLEMENTATION

A. Experimental set-up

The experimental setup is displayed in Fig. 3. The signal is generated from a mode-locked oscillator delivering chirped pulses with a spectral bandwidth of 8 nm centered at 1053 nm (Fig. 4-grey line). The repetition rate is 38.5 MHz and the average power is 25 mW. The linearly chirped pulse has a duration of 4 ps at the output of the oscillator and then is additionally stretched by an Öffner type stretcher to 960 ps (dispersion~120 ps/nm). The repetition rate is decreased to 10 kHz with an acoustooptic modulator (AOM) to match the one of the pump. The pump source is obtained by a CW laser with a wavelength at 1034.2 nm. Another AOM generates a pulse train at 10 kHz with a pulse duration of 140 ns which is synchronized with the pulse train from the oscillator. A first electrooptic modulator (EOM) modulates the phase with a sine shape at 2 GHz to increase the stimulated Brillouin scattering threshold. In general, the pump chirp affects the signal and idler [START_REF] Bigourd | Ultra-broadband fiber optical parametric amplifier pumped by chirped pulses[END_REF][START_REF] Wiberg | Idler Chirp Optimization in a pulse-pumped parametric amplifier[END_REF][START_REF] Vanvincq | Ultrabroadband fiber optical parametric amplifiers pumped by chirped pulses. Part 1: analytical model[END_REF][START_REF] Fourcade-Dutin | Ultrabroadband fiber optical parametric amplifiers pumped by chirped pulses. Part 2: sub-30 fs pulse amplification at high gain[END_REF]. In our condition and for all pump pulse shapes used in following, this impact cannot be observed since the pump spectral spread (~28 GHz) is negligible compared to the signal bandwidth. Then, the pump pulse is amplified in a first ytterbium doped fiber amplifier (YDFA) and then is shortened to 1.3 ns by another electrooptic modulator driven by an electric pulse that can be shaped by an arbitrary wave generator. This system delivers programmable temporal waveforms with a 12 bit amplitude resolution, 100 ps steps over the pulse envelope and a 200 ps rising time. Finally, the pulse energy is raised up to 2 µJ in two YDFAs. The excess of amplified stimulated emission (ASE) is reduced by filters with a spectral width of 1 nm (FWHM) after each YDFA. A second AOM is also used to temporally filter the ASE between successive pulses. In the initial condition, the temporal pulse shape is optimized to be near a flat-top profile. The shape is shown in the inset of Fig. 5 (black line). The maximum variation is less than 10%. Finally, the pump and signal pulses are injected in the photonic crystal fiber (PCF) thanks to a wavelength division multiplexer (WDM). The PCF is an all solid photonic band gap fiber [START_REF] Bouwmans | Fabrication and characterization of an all-solid 2D photonic bandgap fiber with a low-loss region (<20 dB/km) around 1550 nm[END_REF] which has been specifically made to simultaneously obtain a zero dispersion wavelength at around 1033 nm such as the pump wavelength lies in the anomalous dispersion regime and a core diameter of 10 µm [START_REF] Morin | µJ-level Raman-assisted fiber optical fiber parametric chirped pulse amplification[END_REF][START_REF] Bouwmans | Fabrication and characterization of an all-solid 2D photonic bandgap fiber with a low-loss region (<20 dB/km) around 1550 nm[END_REF] to handle higher peak power than standard PCF with a smaller core [START_REF] Bigourd | High-gain, optical-parametric, chirped-pulse amplification of femtosecond pulses at 1 µm[END_REF][START_REF] Fourcade-Dutin | 12 THz flat gain fiber optical parametric amplifiers with dispersion varying fibers[END_REF][START_REF] Bigourd | Ultra-broadband fiber optical parametric amplifier pumped by chirped pulses[END_REF]. The fiber length is 13 m and the amplifier operates in the saturated regime to get a relatively high conversion efficiency. All components are polarization maintaining until the WDM. In front of the PCF, a polarization controller is added to align the pump and signal to one of its birefringent axes. At the output, the spectrum is recorded by an optical spectrum analyzer (OSA). In the time domain, the pump pulse is recorded with a fast photodiode and a 6 GHz oscilloscope in front of the last YDFAs (Fig. 3). When the pump energy is ~1.07 µJ, the parametric fluorescence spectrum (Fig. 4-black line) has two lobes centered at ~1049 nm and 1020 nm. Therefore, we expect that the small-signal gain is not centered on the signal spectrum. Below 1080 nm, a characteristic Raman gain peak is also observed. The blue line in Fig. 4.a shows the amplified signal and the generated idler when the seed energy is 6.2 pJ in front of the WDM. At the PCF output, the energy of the amplified signal is 224 nJ and is deduced from the spectrum and the total energy measurement. Therefore the integrated (energy) gain is 45 dB and the conversion efficiency is ~21 %. The corresponding spectral gain is shown in Fig. 4.b. The gain is obtained from the ratio between the amplified signal at the PCF output and the input signal. The gain value is in the range of 38-48 dB from 1045 nm to 1055 nm. At higher wavelength part of the signal is not amplified as expected from the fluorescence spectrum.

B. Parametric gain shaping 1) Time to spectrum calibration

In the following, we optimize the gain shape through the modification of the pump profile. As the signal is chirped, each part of the pump envelope amplify a portion of the signal spectrum. Therefore, the synchronization and the delay between the pump and the signal is highly important. In fact, the relationship between the time and the spectrum corresponds to the chirp rate of the signal. In order to measure carefully the temporal distribution of the spectral gain, the chirp of the signal must be exactly calibrated. The chirp can be measured directly with the FOPA by applying a dip in the pump profile. When the signal is simultaneously injected in the FOPA, a local decrease is observed in the spectral domain and therefore the chirp should be inferred. For example, the inset of Fig. 5 (measured with the photodiode in Fig. 3) shows a hole placed at the center of the pump profile (numerically centered at 0 ps). The peak power is decreased by 66 % at the center and the width (FWHM) is ~264 ps. In the spectrum measured at the PCF output, the hole appears at the center of the signal spectrum; i.e p=1053 nm. When the temporal position of the hole is at -200 ps and ~200 ps, the corresponding spectral holes occur at p ~1055 nm and p~1051.5 nm, respectively. Fig. 5.b (black circles) displays the relation between the time and the spectrum when the hole is tuned from -200 to 400 ps. For a longer scan, the spectral location of the hole is difficult to obtain exactly because the gain decrease at these wavelength even if the dip is not applied (Fig. 4). The chirp of the stretcher obtained from the geometric parameters are also displayed in Fig 5 .b and the calibrated value obtained directly with the set-up is in very good agreement.

The creation of the hole in the spectrum can also have other interests. For example, the application of a spectral hole in a seed is very interesting to limit the gain narrowing appearing in a laser amplifier [START_REF] Mecseki | Flat-top picosecond pulses generated by chirped spectral modulation from a Nd:YLF regenerative amplifier for pumping few-cycle optical parametric amplifiers[END_REF].

2) Example of gain shaping

From the previous section, we have confirmed that a local variation of the pump influences directly the spectral gain and the chirp has been measured to calibrate the time to wavelength distribution. Therefore, the spectral gain can be locally tuned by applying a special shape in the pump profile. For example, the gain bandwidth can be enhanced when the instantaneous power From the calibrated time-frequency distribution, we deduced that the instantaneous power should be higher at around -200 ps (Fig. 6.a-blue line). In this case, the total energy is 1.07 µJ and is identical to the one used with the flat-top profile. Consequently, the gain has been extended to higher wavelength (Fig. 6.c-blue line) and the bandwidth is now increased by ~30 % from its initial value. This gain bandwidth improvement is a very interesting step forward to amplify an ultra-short pulse in a FOPA with a chirped signal [START_REF] Bigourd | High-gain, optical-parametric, chirped-pulse amplification of femtosecond pulses at 1 µm[END_REF]. Alternatively, the amplified signal can be injected in a laser chain. The required bandwidth is defined by the laser gain material. Therefore, the signal (from the FOPA) does not need to have a bandwidth considerably larger than the one from the laser amplifier. Thus, the pump can be temporally narrowed in order to get a higher peak power (for the same energy) and gain. For example, the gain bandwidth is reduced to 4 nm (Fig. 6.c-red line) by applying a bell shape in the pump profile with a temporal width of ~400 ps centered at 0 ps (Fig. 6.a-red line). The maximum gain centered at 1053 nm increases by 25 dB although we intentionally decrease the current of the last laser amplifier in order to avoid any damage. The energy is 0.7 µJ. The relevant gain can be measured between ~1050 and ~1056 nm. Outside this range, the pump pulse does not overlap with the signal spectrum, and thus no amplification can occur. The observed spectral signal corresponds to the fluorescence.

IV. CONCLUSION

We have investigated the shaping of parametric gain in a PCF thanks to the modification of the pump temporal profile. The local change of the power influences directly the instantaneous parametric gain. With this technique, we demonstrated the possibility to enhance or narrow the bandwidth by tuning locally the pump peak power at a specific time. This scheme can be implemented in front of a laser system to prepare or tune the spectrum of the seed. In addition, the shaping to a specific gain profile is easy and straightforward when the signal has a known and calibrated chirp.

  This work was supported by the ANR (FiberAmp project ANR-16-CE24-0009, ANR-11-EQPX-0017, ANR-11-LABX-0007), the EUR EIPHI program (ANR-17-EURE-0002), the Franche-Comté and Hauts de France councils and the European Regional Development Fund (ERDF) (CPER Photonics for Society P4S, HEAFISY). D. Bigourd, C. Fourcade-Dutin and H. Maillotte are with the Institut FEMTO-ST, Département d'Optique, UMR 6174 CNRS-Université Bourgogne Copyright (c) 2018 IEEE

  1.b) and the maximum gain is obtained at 3.6 THz. For a Gaussian shape (Fig. 1.a green line) and L=6 m, the gain changes with time and is maximum at the center (Fig. 1. d). For a pump with an asymmetric hole in its flat-top profile (Fig. 1.a red line), the gain is decreased locally at a given time (Fig. 1. f). For a longer fiber (L=13 m), the temporal distributions of the gain (Fig. c, e, g) have different

Fig. 2

 2 Fig. 2 Temporal gain for a CW signal centered at 3.6 THz (solid lines) or a pulse with a chirp rate equals to 1.2 THz/ns centered at 3.6 THz (dashed lines) for (a) L=6 m or (b) L=13 m. The pump profiles are shown in Fig. 1.a.

Fig. 1 .

 1 Fig. 1. (a) Pump shapes: flat-top (blue), Gaussian (green) or asymmetric (red) profiles. (b-g) Temporal distributions of the spectral gain (dB) for the three pump profiles. The green and white dashed lines correspond to the chirped signal or a CW to calculate the temporal gain (see Fig. 2).

Fig. 4 (

 4 Fig. 4 (a) Injected signal (grey line), fluorescence spectrum (black line), amplified signal (blue line). (b)Spectral gain of the signal.

Fig. 3 .

 3 Fig.3. Experimental set-up. AOM acoustooptic modulator, EOM electrooptic modulator, AWG arbitrary wave generator, YDFA Ytterbium doped fiber amplifier, WDM wavelength division multiplexer

Fig. 5 (

 5 Fig. 5 (a) Amplified signal when a hole is applied in the pump temporal profile (inset). Delay=+200 ps (green), 0 ps (red), -200 ps (blue) (b) Spectral dip in the signal as a function of the time.

Fig. 6 (

 6 Fig. 6 (a) Flat top (black line), structured (blue line) and Gaussian profile (red line) of the pump (b-c) Corresponding spectral gain.

  3 W -1 km -1 . The simulations is

Franche-Comté, 25030 Besançon, France (Damien.bigourd@femto-st.fr). P.

Morin

, J. Dubertrand and E. Hugonnot are with the Commissariat à l'Energie Atomique et aux Energies Alternatives, Centre d'´Etudes Scientifiques et Techniques d'Aquitaine, 33116 Le Barp Cedex, France (emmanuel.hugonnot@cea.fr) Y. Quiquempois, G. Bouwmans are with CNRS-Université Lille 1, PhLAM/IRCICA USR 3380/UMR 8523, F-59655 Villeneuve d'Ascq Cedex, France (Geraud.Bouwmans@univ-lille1.fr ) Parametric gain shaping from a structured pump pulse D. Bigourd, P. Morin, J. Dubertrand, C. Fourcade-Dutin, H. Maillotte, Y. Quiquempois, G.


A Rjeb

J-P Gayon

S Norre

Fleet-sizing of robots in a logistics warehouse Transport operation between reception area and storage area

Keywords: Fleet-sizing, automated warehouse, internal logistics, transport operations, robots

We are interested in the problem of determining the size of a fleet of robots able to perform transport operations in a logistics warehouse. The operations are divided into several phases: loading, loaded travel, unloading and empty travel. We first determine the minimum number of robots needed to transport a set of homogeneous loads from point A to point B over a determined time interval. For the same number of robots, many assignments exist. We determine the assignment that minimizes the makespan, which is also the one which balances the loads between robots. If we consider an infinity of pickup stations, we obtain a simple analytical expression to find the optimal number of robots. With a single pickup station, we present a simple algorithm that computes the optimal fleet size. Finally, we show that the extension with heterogeneous loads can be formulated as a bin packing problem.

INTRODUCTION

The concept of multi-robot system (MRS) appeared in 1989 in order to replace a large robot by a set of small robots [START_REF] Mouad | Architecture de COntrôle/COmmande dédiée aux systèmes Distribués Autonomes (ACO 2 DA): application à une plate-forme multi-véhicules[END_REF]. For multi-robot systems in the context of logistics warehouses, the literature has since been enriched by focusing in particular on the following fields of study: guide-path design, fleet-sizing, vehicle scheduling, positioning, battery management, deadlock management.

One of the major problems of logistics warehouses is the problem of fleet-sizing which consists in determining the optimal number of vehicles able to perform all requested tasks in a given time interval with a minimum total cost [START_REF] Sinriech | An economic model for determining agv fleet size[END_REF]. Indeed, the number of vehicles greatly influences multi-robot systems performance and determining the right type and number of vehicles is essential.

The fleet-sizing problem has two aspects: an operational aspect and an economic aspect [START_REF] Sinriech | An economic model for determining agv fleet size[END_REF]. For the operational aspect, [START_REF] Egbelu | The use of non-simulation approaches in estimating vehicle requirements in an automated guided vehicle based transport system[END_REF] presents six strategic factors that influence strongly the estimation of the number of vehicles: (1) The system layout, (2) the location of load transfer points, (3) trip exchanges between workcenters per unit time, (4) vehicle-dispatching strategy, (5) system reliability, (6) travel speed.

For the economic aspect, the total costs of robots can be divided as follows [START_REF] Sinriech | An economic model for determining agv fleet size[END_REF]:

Thesis co-financed by the European Union within the framework of the European Regional Development Fund ERDF

• Operating costs including everything that is related to the operations such as maintenance costs, energy costs, etc. • Design costs including the software costs (scheduling, routing, dispatching, etc) and the hardware costs (batteries / chargers, communication links, etc). • Guide-path design cost related to the the guide-path configuration.

The literature considers the fleet-sizing problem in different ways (see e.g. [START_REF] Choobineh | Fleet sizing of automated guided vehicles: a linear programming approach based on closed queuing networks[END_REF] and [START_REF] Ganesharajah | Design and operational issues in agv-served manufacturing systems[END_REF] for reviews).

According to [START_REF] Choobineh | Fleet sizing of automated guided vehicles: a linear programming approach based on closed queuing networks[END_REF], the sizing models presented in the literature can be classified into four categories of approaches: (1) discrete and continuous event simulation models, (2) calculus approaches including a set of simple computations, (3) deterministic operations research models including linear and integer programming, (4) stochastic models including queuing theory.

Since the first research works dealing with this fleet-sizing problem, the analytical models used to determine the optimal size of the fleet are not entirely reliable in terms of optimization. Indeed, analytical models today are only used to find an initial number of robots which is not necessarily optimal before continuing with the sizing procedure.

In order to refine the analytical results, most of the work is based on the number of robots found analytically to launch an optimization study through algorithmic or simulation models.

Among the most used models to start the sizing procedure, we can cite [START_REF] Egbelu | The use of non-simulation approaches in estimating vehicle requirements in an automated guided vehicle based transport system[END_REF] model. This model calculates the number of robots in the fleet according to the following expression: m * = total operating time time interval Although this model takes into account several constraints (empty circulation of vehicles, number of pickup / delivery stations per zone ...), it does not take the integer values into consideration which strongly impacts the optimality of the solution.

Recently, [START_REF] Vivaldini | Integrated tasks assignment and routing for the estimation of the optimal number of agvs[END_REF] improved this model by taking into account integer values. On the other hand, this model remains limited from an optimization point of view since it does not take into account the waiting times caused by several problems (unavailability of the load to be transported, loading station occupied, etc). Indeed, [START_REF] Vivaldini | Integrated tasks assignment and routing for the estimation of the optimal number of agvs[END_REF] consider the number found by this model as a starting point for a second optimization step.

The model of [START_REF] Choobineh | Fleet sizing of automated guided vehicles: a linear programming approach based on closed queuing networks[END_REF] has been used by several researchers to find an initial value for the size of the fleet. This model was used recently by [START_REF] Chawla | Automatic guided vehicles fleet size optimization for flexible manufacturing system by grey wolf optimization algorithm[END_REF] and [START_REF] Chawla | Material handling robots fleet size optimization by a heuristic[END_REF] in order to show that the values found analytically can be optimized through optimization tools like (gray wolf optimization algorithm, particle swarm algorithm).

TRANSPORT OF HOMOGENEOUS LOADS

In this section, we present a transport model for loads with the same type.

Notations and assumptions

Assumptions

-The robot parking is located at point A -Transport is unitary: each robot carries at most a single load at a given time -The loads to be transported are always available -Deadlock during travel and idleness of robots are not taken into account -The battery charging problem is not taken into account

Objective function

The purpose of this study is to determine the minimum number of robots m * allowing to transport the n loads in the time interval [0, T ].

Many load assignments can exist for the same number of robots. We will denote by C j the time at which robot R j has terminated its last transport cycle. We will also denote by C max = max j C j the makespan which represents the time at which all loads have been transported and all robots are back in A. We will focus in particular on the assignment which minimizes the makespan but also on the one which distributes the loads in the most balanced way between robots.

Cycle time of a robot We will denote by p the cycle time of a robot for one load (loading, loaded travel, unloading and empty travel):

p = d AB • 1 v l + 1 v e + t l + t u

Many pickup / delivery stations

We assume in this part that there are at least as many pickup / delivery stations as there are robots at the optimum (n A ≥ m * and n B ≥ m * ). Under this assumption, there is no waiting time for loading or unloading.

We can notice that there is a solution to the problem if and only if the cycle time p is less than or equal to the time interval [0, T ]. We will assume in the rest of this paper that p ≤ T in order to place ourselves in the case where a solution exists.

Optimal number of robots

The maximum number of loads n max transported by a robot on the time interval [0, T ] must satisfy the inequality (n max p ≤ T ). The number n max must be integer, which implies that

n max = T p .
where x is the floor function of x. We will also denote by x the ceiling function of x.

The maximum number of loads transported by m robots, denoted by N max (m), is then

N max (m) = m T p .
In order to be able to transport the set of n loads, it is necessary that

N max (m) ≥ n, so m ≥ n T p
.

(1)

The optimal number of robots is the smallest value of m which satisfies (1), so we have

m * =     n T p     . ( 2
)
Robots load balancing There can be several ways of assigning loads between the m * robots. The most balanced solution is the one that assigns the load with the objective of balancing the use of robots as much as possible. More precisely, it is the solution which minimizes the maximum difference between the different end dates max i,j |C i -C j |. In the most balanced solution, the end date of the latest robot finishing its cycle is C max where

C max = n m * p
We notice that the most balanced solution is also the one which minimizes the C max .

Let us now summarize our main results in the following theorem.

Theorem 1. When considering an infinity of pickup stations, the optimal number of robots is

m * =     n T p     .
In the most balanced solution, noting s = n m * m * -n, we have

• s robots carrying n m * loads • (m * -s) robots carrying n m * loads. Finally, the most balanced solution is also the one that minimizes the makespan where

C max = n m * p.
Figs. 2 and Fig. 3 present two solutions for a transport model with an infinity of pickup stations and an optimal number of robots m * = 4. In the most balanced optimal solution (Fig. 2), 1 robot carries 4 loads and 3 robots carry 3 loads. In the unbalanced solution (Fig. 3), 3 robots carry 4 loads and 1 robot carries 1 load.

We can notice that the two solutions have the same makespan C max = 28. On the other hand, the most balanced solution minimizes the differences between all the C j which guarantees the balanced use of robots for this model. Continuous approximation If we allow a non-integer number of robots and loads, the optimal number of robots can be simplified to m * = np/T and all the robots carry T /p loads with a makespan C max = T . These results are consistent with [START_REF] Egbelu | The use of non-simulation approaches in estimating vehicle requirements in an automated guided vehicle based transport system[END_REF]. This approximation is good if p is much shorter than T and the optimal number of robots is large, i.e. if np/T is large.

Single pickup station

We assume in this section that the system contains a single pickup station at point A and at least one delivery station at point B, thus, n A = 1 and n B ≥ 1.

With a single pickup station, only one robot can load at a given time. We will denote by R j the robot starting in j-th position with (1 ≤ j ≤ m * ). During the first load, R 1 does not wait, R 2 waits t l and more generally R j waits (j -1)t l . For the following loadings, we must distinguish three cases:

• Case 1 : Waiting for first loading only (p ≥ m • t l )

In this situation, when R 1 returns for the first time to A, robots R 2 , • • • , R m have completed their first loading. In this case, R 1 never waits to load and the others wait only for their first loading (see example in Fig. 4). In this situation, when R 1 returns for the first time to A, R m is being loaded. Note that the robots wait for (m • t l -p) from their second loading. We observe for example in Fig. 5 that R 1 is forced to wait, on its return to A, for the end of loading of R 5 . Here, when R 1 returns for the first time in A, R m has not yet started its first loading. Hence, there is at least one useless robot. We will therefore not consider this case subsequently.

Optimal number of robots

We denote by n max j the maximum number of loads transported by R j . In order to find n max j , we will distinguish the two first cases identified previously.

-Case 1 : The number of robots m is integer, so, the condition m • t l ≤ p is equivalent to m ≤ p t l In addition, the time required for robot R j to transport n j loads must be less than T , so we can write this condition as p • n j + (j -1) • t l ≤ T .

As n max j is an integer, we have

n max j = T -t l • (j -1) p -Case 2 :
As m is an integer, the case p t l < m < p t l + 1 is only possible if p t l is not an integer. It follows that the condition is equivalent to m = m 0 with m 0 = p t l . The time required by R j to transport n j loads must be less than T , i.e.

n j • p + (j -1) • t l + (n j -1)(m 0 • t l -p) ≤ T.

It follows that

n max j = T -t l (j -1) -p + m 0 • t l m 0 • t l
Then, we get the expression of N max (m) = m j=1 n max j , the maximum number of loads transported by m robots :

N max (m) =                m j=1 T -t l (j -1) p if m < m 0 m0 j=1 T -t l (j -1) -p + m 0 t l m 0 • t l if m ≥ m 0
We conclude that the optimal number of robots is :

m * = min [m ∈ {1, .., m 0 } : N max (m) ≥ n]
Robots load balancing In this section, we use the found minimum number of robots m * and the balanced distribution of tasks presented previously (Theorem 1) which always remains valid for this case.

Let n j be the number of loads transported by robot R j during the interval [0, T ]. We have

n j =        n m * if j ∈ {1, • • • , m * -s} n m * if j ∈ {m * -s + 1, • • • , m * } et s = n m * m * -n
We can express the end date of the last transport cycle C j for any robot R j :

-Case 1 :

The operating time of R j includes the transport time of n j loads as well as the waiting time for the first loading (see Fig. ( 4)).

It follows that

C j = n j • p + (j -1) • t l -Case 2 :
The operating time of R j contains the transport time of n j loads, the waiting time for the first loading (except the first robot R 1 ) which is proportional to the loading time t l and the waiting time for any other loading which is proportional to the quantity of loads (m * • t l -p) (see Fig. 5). Hence, we have

C j = n j • p + (j -1) • t l + (n j -1)(m * • t l -p)
In the end, we can summarize our results in the following theorem.

Theorem 2. When considering a single pickup station, the optimal number of robots is

m * = min [m ∈ {1, • • • , m 0 } : N max (m) ≥ n] with m 0 = p t l and N max (m) =                m j=1 T -t l (j -1) p if m < m 0 m0 j=1 T -t l (j -1) -p + m 0 t l m 0 t l if m ≥ m 0
The most balanced solution has a makespan

C max = max 1≤j≤m * (C j )
where

C j = n j p + (j -1) t l if m * < m 0 n j p + (j -1) t l + (n j -1)(m * t l -p) if m * ≥ m 0
When applying Theorem 2 to the case with zero loading time (t l = 0), we re-obtain Theorem 1 as the condition (p ≥ m • t l ) is always true.

TRANSPORT OF HETEROGENEOUS LOADS

In this section, we extend our results to heterogeneous loads in the case of many stations. Fig. 6 presents an example with two types of loads. We assume that the weight of the load impacts the speed of a loaded robot. Without this assumption the problem remains the same as the one of section 2.2 since we use the same type of robots.

We will use the following notations:

-n : total number of loads (denoted 1, • • • , n) -v li : travel speed of a robot carrying load i

The cycle time of the load i is

p i = d AB • 1 v li + 1 v e + t l + t u
We assume that p i ≤ T for all i, otherwise there is no solution to the problem.

This problem can be formulated as a bin packing problem. We can consider that we have n objects (loads). Each object i has a size p i . The problem is to put all the objects in bins (robots) of size T (time limit). The smallest number of bins is the optimal number of robots m * .

This problem can be formulated as an integer linear programming. Denote by x ij the binary variable equal to 1 if load i is transported by robot R j , 0 otherwise. Denote also by y j the binary variable equal to 1 if bin j is used, 0 otherwise.

As there are at most n needed robots, we therefore seek to minimize the number of bins used min n j=1 y j under the following constraints:

n i=1 p i • x ij ≤ T • y j j = 1, • • • , n n j=1 x ij = 1 i = 1, • • • , n y j ∈ {0, 1} j = 1, • • • , n x ij ∈ {0, 1} i, j = 1, • • • , n
There are different methods to solve this problem (see for example [START_REF] Toth | Knapsack problems: Algorithms and computer implementations[END_REF]).

Continuous approximation If we allow a non-integer number of robots and loads, we can simplify the solution of our problem to m * = n i=1 p i /T

CONCLUSION

This paper considers a fleet-sizing problem of robots in a logistics warehouse. An analytical approach allowed us to deal with this problem for a simple model of transporting homogeneous loads between two storage areas in the warehouse, taking into account the number of pickup and delivery stations in the system. For the case with many pickup stations (n A ≥ m * and n B ≥ n A ), the size of the fleet is given by a simple analytical expression (see Theorem 1). For the case with a single pickup station (n A = 1 and n B ≥ n A ), the problem is solved using a simple algorithm (see Theorem 2). We also show how to formulate the problem with heterogeneous loads and many stations by a bin packing problem and an integer linear programming.

It would be interesting to extend our model to the case with heterogeneous robots in capacity or speed.

Fig. 1

 1 Fig. 1 illustrates an homogeneous fleet of robots that transport homogeneous loads from zone A (reception) to zone B (storage). A robot moves loaded from A to B before crossing the same empty distance from B to A.

Fig. 1 .

 1 Fig. 1. Transport model between A and B

  If n m * is an integer, then it is possible to distribute the loads in a perfectly balanced way between the robots by assigning n m * loads to each one. If n m * is not an integer, the robots will transport, in the most balanced solution, either n m * loads, or n m * = n m * + 1 loads. Let s be the number of robots carrying n m * loads. Thus we have (m * -s) other robots carrying n m * = n m * + 1 loads. The total number of loads to be transported being n, we have n = s n m * + (m * -s) n m * + 1 . which implies that s = n m * m * -n. We notice that if n m * is an integer, this last formula remains correct by giving s = 0.

Fig. 2 .

 2 Fig. 2. Balanced solution (n = 13, t l = 1, p = 7, T = 30, m * = 4)

Fig. 4 .

 4 Fig. 4. Waiting for first loading only (p ≥ m • t l ) (n = 9, t l = 1, p = 7, T = 21 , m * = 4)• Case 2 : Waiting for every loading (m-1)t l < p < mt l

Fig. 5 .

 5 Fig. 5. Waiting for every loading (m -1) • t l < p < m • t l (n = 7, t l = 2, p = 9, T = 21, m * = 5) • Case 3 : Unnecessary robots p ≤ (m -1) • t l

Fig. 6 .

 6 Fig. 6. A transport model with two types of loads