Sébastien Grot
email: sgrot@keopsys.com

Lew Goldberg

Pascal Besnard

Yves Jaouën

High performance pulsed coherent fiber laser emitting >1.7kW at 1060nm in ns regime

Keywords: fiber laser, optical pulses, double clad fibers, non linear effects, Yb 3+

Fiber lasers and amplifiers are used in a variety of applications either for scientific (spectroscopy, medicine …) or industrial applications (free space communications, laser marking and drilling …). The combination of doped double clad fibers (DCF) and high power multimode semiconductors laser diodes technologies allows to achieve very high output power in very compact, robust and maintenance free systems. Yb 3+ doped DCF are well suited for 1µm wavelength amplification. In pulsed regime, achievable peak power can be strongly limited by nonlinear effects such as Kerr effect, Stimulated Raman Scattering (SRS) or Stimulated Brillouin Scattering (SBS). Consequently, the optimisation of optical amplifier architecture is required. In this paper, we demonstrate performances obtained for the generation of 2ns optical pulses up to >1.7kW peak power in a Master Oscillator Power Fiber Amplifier (MOPFA) configuration. The laser seed signal at 1060nm is emitted out of a single longitudinal mode source with spectral linewidth <0.2nm. The pulse repetition rate can be changed between 3 and 30MHz. The high power stage, based on a 2-stages architecture, allows to deliver >10W average output power with a good beam quality (M 2 <1.2). No significant limitation due to nonlinear effects of the type of the Kerr effect or SRS appears by means of the optimisation of the final stage's fiber parameters. Results, such as a concentration of more than 80% of the total output power in a 1nm window around the central wavelength and above all an excellent conservation of the spectral properties of the seed source are demonstrated for a peak power of >1.7kW. These high performances are obtained in a fully-integrated device.

1-INTRODUCTION

The combination of the DCF technology with high power semiconductor laser diodes allows to develop high power optical fiber sources in the 1-100W range. Yb 3+ doped DCF are very efficient to develop lasers and amplifiers in the 1µm wavelength range [START_REF] Pask | Ytterbiumdoped silica fiber lasers : versatile sources for 1-1.2 µm region[END_REF][START_REF] Bordais | Double-clad 10W Yb 3+ -doped fiber master oscillator power fiber amplifier for He 3+ optical pumping[END_REF]. Intrinsic properties of DCF authorize a good coupling efficiency of the light emitted out of multimode high power laser diodes. Pulsed regime allows generation of peak power >1kW. One of the most important limitation concerns nonlinear effects [START_REF] Jaouën | High power cladding-pumped Er 3+ /Yb 3+ fiber amplifiers: technologies, performances and impact of nonlinear effects[END_REF]. To achieve high energy extraction and increase the nonlinear effects threshold the use of large mode area (LMA) fibers are required. In this paper, we present results obtained while amplifying a very narrow laser signal up to >1.7kW in 2ns pulses. The most difficult here was to achieve such high peak power while maintaining the spectral source properties (linewidth <0.2nm) and spatial properties (M 2 <1.2).

2-POWER FIBER AMPLIFIER OPTIMISATION

The architecture of the device is presented on fig. 1. To achieve output power higher than few watts average, a 2-stages amplifier is required. The preamplifier stage is pumped by a broad stripe laser diode at 975nm using the V-groove side pumping technique [START_REF] Ripin | High efficiency side-coupling of light into optical fibres using imbedded V-grooves[END_REF]. The output power of the preamplifier stage (~200mW), corresponding to ~33W peak power, is high enough to efficiently saturate the power amplifier stage. To increase the nonlinear effects threshold, we have retained a slightly multimode Yb 3+ doped fiber for the booster stage. This DCF has a core diameter of 13µm and a numerical aperture of 0.13. The length of DCF used is 5m and corresponds to a pump power absorption >10dB (absorption ~2.5dB/m at 920nm). A fiber pigtail of quasi single mode fiber at the signal wavelength (short piece of SMF28, mode field diameter <9µm) is used to preserve the beam quality parameter (M 2 <1.2). The source generates 2ns pulses at a repetition rate adjustable between 3 and 30MHz. The limitation due to nonlinear effects can appear for peak power higher than few watts/m at 1µm (for fibers of ~20-40µm 2 effective mode area). For optical pulses >10 ns, the Stimulated Brillouin Scattering (SBS) effect is the most important non linear effect [START_REF] Teodoro | Diffraction limited, 300-kW peak-power pulses from a coiled multimode fiber amplifier[END_REF][START_REF] Grot | High power Yb 3+ -doped double clad fibre amplification at 1127nm: 977 against 920nm pumping comparison[END_REF]. In our case, the dominant non linear effects are Kerr effect and SRS respectively. Even if there is a penalty in terms of non linear effects, we had to consider the use of commercially available 920nm high power fiber pigtailed laser diodes. The penalty in terms of non linear effects threshold is balanced by the possibility of using uncooled laser diodes with no wavelength stabilization. The booster stage is pumped by a combination of high power laser diodes using a fiber bundle for a total pump power budget of 13W. A contra-pumping scheme allows increasing non linear effects threshold [START_REF] Kulcsar | Multiple-Stokes stimulated Brillouin scattering generation in pulsed high power double-cladding Er 3+ /Yb 3+ co-doped fiber amplifier[END_REF]. However, we have retained a co-pumping scheme to avoid the penalty of the insertion loss of the fiber bundle at the output of the booster amplifier. This configuration allows to reduce pump power budget.

3-SPECTRAL BROADENING AND OPTIMISATION OF THE SOURCE PARAMETERS

In pulsed regime, for peak powers higher than a few watts/m, an optical spectrum broadening can occur due to Kerr effect. This nonlinear effect results in a spectral broadening of the seed source in dependence with output peak power. The seed source is an external cavity laser with a linewidth <0.2nm. The most important contribution to linewidth broadening is the SPM effect (Fig. 2). The spectrum broadening can be estimated using the expression [START_REF]Nonlinear fiber optics[END_REF]

eff L t peak P Yb SPM ∂ ∂ ≅ ∆ γ ω (1)
where γ Yb is the nonlinear Kerr parameter and L eff the effective length. The nonlinear Kerr parameter γ Yb ~10 W -1 .km -1 has been determined using by the measurement of the four-wave mixing (FWM) efficiency between two co-propagating signals with co-linear polarization. The FWM phase matching condition is realized for <1nm wavelength spacing. Assuming an uniform amplifier gain distribution along the Yb 3+ doped fiber length (L=5m), the effective length is L eff ~1.3m [START_REF] Grot | High power Yb 3+ -doped double clad fibre amplification at 1127nm: 977 against 920nm pumping comparison[END_REF]. A ~0.08 nm spectral broadening corresponds to an output peak power of 1.7 kW. In practise more than 80% of the overall power at 1.7kW is contained in a 1nm window around the central wavelength and even more than 90% in 5nm around it. In the meantime, the laser linewidth and spectral purity of the source are preserved (Fig. 3).

CONVERSION EFFICIENCY AND PEAK POWER

The booster amplifier stage architecture has been designed to maximize power conversion efficiency. We have measured better than 75% efficiency: 10W output power at 1060nm is obtained for a pump power of 13W (Fig. 4). Peak power is higher than 1.7kW in 2ns pulses at a repetition rate of 3MHz (Fig 5). Thanks to the geometry of the 2 nd clad of the DCF and its numerical aperture, we have achieved a quasi-conservation of the fundamental mode propagation throughout the booster. The signal to noise ratio is not limited by SRS and remains >28dB even when peak power is >1.7kW. It is only limited by the amplified spontaneous emission around 1075nm.

Peak power (kW) At such a high level of peak power, SRS nonlinear effect appears. In case of standard Yb 3+ -doped DCF (i.e. mode field diameter of ~6µm) are used, this phenomenon is visible at peak powers as low as few hundreds of watts [START_REF] Jaouën | High power cladding-pumped Er 3+ /Yb 3+ fiber amplifiers: technologies, performances and impact of nonlinear effects[END_REF][START_REF] Champert | Multiwatt average power continua generation in holey fibers pumped by kilowatt peak power seeded ytterbium fiber amplifier[END_REF]. The large mode area fiber used in this power amplifier allows to increase the SRS power threshold up to 1.5kW. Even at 1.7kW most of the total average power is contained in a 1nm window around the central wavelength. At 30MHz repetition rate, SRS disappears as peak power is <160W.

CONCLUSION

We have demonstrated the feasibility of a high power pulsed Yb 3+ doped double-clad fiber laser in a MOPFA configuration capable of delivering >10W average power in 2ns pulses. Booster stage fiber parameters optimization allows to achieve >1.7kW at a repetition rate of 3MHz. No pulse distortion in temporal domain is observed even at this level of peak power. The beam quality (M 2 <1.2) proves that the preservation of the propagation of a quasi-gaussian mode is achieved. The performance also results in the preservation of the spectral purity of the source even at high peak power. To our knowledge, these good specifications constitute the best performances with the constrains of the 920nm pumping wavelength and co-pumping scheme of the booster stage and are achieved in a fully integrated device.

Figure 1 :

 1 Figure 1: Architecture of the MOPFA fiber laser designed to achieve >10W average output power in 2ns pulses at 3-30MHz repetition rate

Figure 2 :

 2 Figure 2: Spectral broadening due to SPM

Figure 3 :

 3 Figure 3: Optical spectra of the external cavity laser diode before and after amplification

Figure 4 :

 4 Figure 4: Average output power and peak power over +/-0.5nm around central wavelength versus injected pump power

Figure 5 :

 5 Figure 5: Optical spectra obtained at 10W average output power at 3 and 30MHz repetition rate.


Nicolas Olivier

Miguel Angel Luengo-Oroz

Louise Duloquin

Emmanuel Faure
email: emmanuel.beaurepaire@polytechnique.edu

Thierry Savy

Israël Veilleux

Xavier Solinas

Delphine Débarre

Paul Bourgine

Andrés Santos

Cell lineage reconstruction of early zebrafish embryos using label-free nonlinear microscopy

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Cell Lineage Reconstruction of Early Zebrafish Embryos Using Label-Free Nonlinear Microscopy

Nicolas Olivier, 1 * Miguel A. Luengo-Oroz, 2 * Louise Duloquin, 3 * Emmanuel Faure, 4 Thierry Savy, 4 Israël Veilleux, 1 Xavier Solinas, 1 Delphine Débarre, 1 Paul Bourgine, 4,5 Andrés Santos, 2 Nadine Peyriéras, 3,6 † Emmanuel Beaurepaire 1 † Quantifying cell behaviors in animal early embryogenesis remains a challenging issue requiring in toto imaging and automated image analysis. We designed a framework for imaging and reconstructing unstained whole zebrafish embryos for their first 10 cell division cycles and report measurements along the cell lineage with micrometer spatial resolution and minute temporal accuracy. Point-scanning multiphoton excitation optimized to preferentially probe the innermost regions of the embryo provided intrinsic signals highlighting all mitotic spindles and cell boundaries. Automated image analysis revealed the phenomenology of cell proliferation. Blastomeres continuously drift out of synchrony. After the 32-cell stage, the cell cycle lengthens according to cell radial position, leading to apparent division waves. Progressive amplification of this process is the rule, contrasting with classical descriptions of abrupt changes in the system dynamics.

Although classical developmental biology is characterized by qualitative descriptions, recent work underlines the requirements for precise measurements to enable formal reconstruction integrating the genetic, molecular, and cellular levels of organization (1)(2)(3). The optimization of microscopy imaging techniques and improved data algorithmic processing are key issues in such reconstructions. Parallelized linear microscopy such as light-sheet fluorescence mi-croscopy provides fast imaging but suffers from loss of information with depth (4). Point-scanning two-photon microscopy provides deeper imaging (5) but exhibits slower frame rate, compromising automated individual cell tracking in whole organisms (6). Furthermore, the usual implementation of these two paradigms does not allow homogeneous illumination in spherical samples, leading to a difficult tradeoff between the detection of deep structures and illumination-induced perturbation in outer layers. Finally, relying on fluorescent staining of biological structures brings additional artifacts and limitations. Exploiting the intrinsic optical nonlinear properties of the sample is a valuable, although challenging, alternative. Second-harmonic generation (SHG) is obtained from dense noncentrosymmetric structures such as oriented microtubule assemblies (7-9), including mitotic spindles (8,[START_REF] Sun | [END_REF]. Third-harmonic generation (THG) is obtained from optical heterogeneities (11)-such as the interface between an aqueous medium and a lipidic, mineralized, or absorbing structure (12)-and allows morphological imaging of small organisms [START_REF] Sun | [END_REF]13).

Here, we show that combining SHG and THG imaging of unlabeled embryos with a scanning scheme matching embryo morphology provides adequate three-dimensional (3D) imaging over time for the automated reconstruction of cell behavior during zebrafish embryo cleavage stages (14). Ad hoc image analysis strategies for cell position, division, and shape identification were used to produce a complete and validated lineage tree for a cohort of six zebrafish embryos until the 1000-cell stage, annotated with minute-level division timing, micrometer-accuracy cell coordinates, and shape characteristics. These data provided a quantitative spatiotemporal description of the wavelike division cycles and allowed the construction of a prototypic digital blastula. The cycle duration of sister cells exhibited variability that did not correlate with cell volume, revealing unexpected cell division asynchrony and asymmetry from the first division cycles and leading to increasing cell heterogeneity by the time of midblastula transition (MBT) (14).

An appropriate image acquisition scheme was devised to provide high-resolution time-lapse imaging of intrinsic SHG and THG signals (Fig. 1 and supporting online material). Excitation in the 1.2-mm range reduced nonlinear endogenous absorption by the sample and allowed simultaneous two-photon-excited fluorescence (2PEF) imaging of red fluorescent proteins (Fig. 1, B andC) for control experiments. When imaging a spherical embryo, scattering and aberrations typically result in reduced signal at the center of each plane (Fig. 1E). We therefore scanned each plane of a half-sphere along a spiral trajectory with variable speed to spend more time imaging the innermost cells (Fig. 1, D to F, and fig. S1). This conformal strategy provided optimal acquisition time and minimal photoperturbation (fig. S2). SHG and THG signals were co-optimized by using rotating linear incident polarization. In addition, because THG contrast from a specific structure depends on its size relative to the focal volume (15), moderate focusing (3.5-mm Z-resolution) was used to highlight cell interface compared with smaller subcellular structures (Fig. 2D).

Combining the conformal scanning scheme described above, sensitive detection, infrared excitation wavelength, and appropriate focusing and polarization conditions allowed homogenous detection of mitotic spindles and cell and tissue phenotypic features in the whole unlabeled zebrafish embryo during cleavage stages. The blastoderm was contained in a half sphere of 440-mmradius imaged with a temporal resolution of 80 s and a volumetric pixel size of 2 by 2 by 4 mm, suitable for further automated reconstruction of the cell lineage tree.

The intrinsic THG signal revealed a number of structures and dynamic processes (Fig. 2, A to J, and movies S1 to S6) and highlighted cell contours even better than membrane staining by farnesylated mCherry, indicating the presence of a sizeable intercellular space as corroborated by numerical simulations (Fig. 2, K to M). Highernumerical aperture THG imaging produced signal discontinuities or double interfaces at cell boundaries (Fig. 2D), indicating local distance variation between cells. This interpretation is reinforced by the absence of THG signal from the lateral junction of the enveloping layer (EVL) (fig. S3). THG also allowed visualizing trafficking phenomena, including cell membrane completion at the 32-cell stage (movie S2) and nuclear envelope dynamics during mitosis (movie S3). In addition, when embryos were overilluminated, photo damage correlated with enhanced intracellular THG signal, serving as a control of embryo physiological condition (fig. S2). THG also revealed yolk lipid platelets interface and allowed visualizing the phenomenology of cytoplasmic streams during the first cell cycles (movie S1). Combining THG/SHG/2PEF in ßactin:H2B/mcherry transgenic embryos ( 16) allowed correlating SHG and THG signals with mitosis phases (Fig. 2, G to I, and movie S4). As shown in Fig. 2O, SHG intensity exhibited a nearly Gaussian temporal profile peaking at metaphase, providing measurement for division timing and cell cycle duration with subminute accuracy.

We designed an original image processing methodology for automatically extracting cell position, membrane geometry, and cell lineage, taking advantage of mitosis metasynchrony and limited cell displacements during the first 10 division cycles (fig. S4 and movies S5 and S6). daughter cells after mitosis was then used to link consecutive sequences. Finally, accurate cell cycle timing was obtained by fitting time-dependent SHG intensity to a Gaussian function (Fig. 2O).

The lineage tree validation is a critical step as tracking errors propagate along the tree. Validation was performed using an interactive visualization interface that allowed comparing raw and reconstructed data at each cell cycle before processing the next one (movie S8). This strategy provided error-free lineage trees with minimal human intervention (fig. S5). Cell shapes were then extracted from THG images by using cells' spatial coordinates as seeds to perform a region growing-based contour detection. A rough shape approximation was first obtained by building a Voronoi diagram, and membrane detection was refined using a viscous watershed algorithm (19) (movie S9 and Fig. 3, B to E). Data sets from six different embryos were processed to obtain the corresponding digital blastulas (Fig. 3, fig. S5, and movies S10 to S13). The phenomenological reconstruction allowed a systematic analysis of spatiotemporal correlations between cell position, cell cycle duration, mitosis duration, cell volume, and cell division orientation. Cells during cycles 1 to 4 were constrained by the absence of basal plasma membrane contacting the yolk, and divisions displayed the known stereotyped orientation with little fluctuation from one embryo to the other (Fig. 3M). Division orientation at cycle 5 started exploring one more spatial direction at 90°relative to the previous division, thus allowing the formation of a partially double-layered blastoderm (Fig. 3, L andN). This quantitative observation contrasts with the classical description presenting the 32-cell stage as a single cell layer (14). From divisions 5 to 8, larger fluctuations occurred in division orientation, resulting in progressively more random orientations (Fig. 3N and fig. S6). By averaging and registering six different embryos, we constructed a prototype of the lineage spatiotemporal deployment (fig. S6 and movie S10). From cycle 3 to cycle 9, cell displacements were limited to cytokinesis, and displacement speed from one division to the other remained constant (fig. S7), suggesting that biochemical dynamics (20,21) do not impact cell displacements until MBT, quantitatively supporting previous observations (18). Consistently, the clones corresponding to the first eight blastomeres observed at the 512-cell stage remained clustered (Fig. 3, F to J, and movie S10).

In the temporal domain, the main dynamical change was the evolution of cell cycle duration as previously observed (22). The cell cycle shortened until division 5 (18.5 T 1.5 min to 15.5 T 1min), followed by lengthening until division 10 (up to 18.5 T 1.5 min) (Fig. 4A). The minimum value correlated with plasma membrane completion at the 32-cell stage, suggesting a specific regime for the division of the large early blastomeres (23,24). Continuous cell cycle lengthening thereafter is hypothesized to result from the degradation of maternal products and S-phase lengthening (21). Furthermore, as early as the 4-cell stage, cell divisions did not occur synchronously, and the delay measured between the first and last division at each cell cycle increased continuously through cleavages (Fig. 4, B andC).

In addition, after the 32-cell stage, cell cycle lengthening increased with the physical distance from the blastoderm surface (Fig. 4, D andE, and figs. S8 to S11), so that in the particular case of sister cells, the deeper tended to divide later than her sister located closer to the surface (fig. S13). We hypothesize that this feature correlates with deposition of maternal components by cytoplasmic streams (movie S1) that might have produced a concentration gradient. In any case, fluctuation in sister cells' division timing and steady lengthening in cell cycle duration as a function of cell position relative to the blastoderm surface are sufficient to explain the progressive appearance of a global wave-like pattern of radial division (Fig. 4, D andE). Calculating the position of this pseudo-wave origin relative to the animal pole indicated an early symmetry breaking variably amplified in some cases (fig. S14). A discussion of the cell division patterns during cleavage stages must also take into account the formation of the yolk syncitial layer (YSL) (25). We observed that YSL nuclei typically individualized at the 512-cell stage and underwent three division cycles with an average timing of 20 to 30 min and displaying a peripheral pattern, thus differing from the radial pattern observed in the blastoderm (movies S14 to S16 and figs. S15 and 16). YSL formation is known to display phenotypic variability (26), but extreme cases of artifactual YSL formation were observed in embryos mechanically constrained by the mounting medium (movie S16 and fig. S17). Altogether, our observations suggest that there is no switch from a radial toward a peripheral pseudo-wave pattern in the blastoderm, in contrast to a recent proposal by Keller et al. (4).

Cell volume measurements indicated that up to cycle 10, the nucleo-cytoplasmic ratio (N/C) was unlikely to trigger cell cycle lengthening as we observed cell volume variability, including random fluctuation between sister cells, with no correlation to cycle duration (figs. S11 and S12). According to other studies, N/C becomes an important factor at later developmental stages (18). Indeed, a pause in cell cycle, regulated by N/C, is described as one of the most prominent features of the MBT (27). This means that by MBT, the N/C variability generated earlier could act as a local fluctuation factor, breaking the global pseudo-wave pattern. Cell cycle lengthening, both globally and as a function of cell position, is the main parameter that should be taken into account for further modeling zebrafish embryo cleavage morphogenesis. Continuity and amplification of the process are the rule, rather than the classically described abrupt changes from synchrony to metasynchrony and asynchrony (14).

Harmonic time-lapse microscopy provides in toto imaging, enabling automated and validated reconstruction of the zebrafish embryo lineage tree during the cleavage period. Because the reconstruction scheme is based on cell cycle metasynchrony and limited cell displacements, we provide a quantitative assessment that cell intrinsic motility does not arise before cell cycle 10, described as the MBT stage (18). The framework presented here automates the phenomenological reconstruction of animal early embryogenesis without the use of fluorescent staining. Conformal acquisition with excitation power and scanning speed matching specimen morphology makes point-scanning multiphoton microscopy an excellent strategy for deep embryo imaging. We anticipate that combining harmonic and fluorescent signals in conformal scanning time-lapse microscopy will provide optimal data sets to achieve the automated tracking of cell trajectories and cell divisions throughout vertebrate embryogenesis. Single-cell analysis along the lineage tree is necessary to provide relevant spatial and temporal correlations underlying emerging patterns. Most of developmental biology still relies on visual inspection and manual work to achieve limited cell tracking, and recent studies have provided only a global estimation of cell movements. In this context, we expect our framework to establish standards in terms of measurements precision and accuracy of automated algorithmic segmentation and tracking procedures. The quantitative phenomenological reconstruction of the zebrafish embryo cleavage period provided here should serve as a reference for further analysis of the interplay between nano-(molecular and genetic) and macro-(biomechanics) level dynamics.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Conformal scanning harmonic microscopy scheme. (A) Schematized zebrafish embryo imaged from the animal pole. (B) Sample mounting for upright microscopy; excitation (Exc.) and detection geometry (2PEF epidetected, THG-SHG transmitted). (C) Energy diagrams and wavelengths involved. (D) Lateral position (left) and depth (right) o ft h ef o c a lp o i n ti n s i d eas p h e r i c a ls a m p l ew h e n using raster scanning (top) versus conformal scanning (bottom). (E) THG imaging with raster (top) and spiral scanning (bottom). Acquisition time is 5 s. Scale bar, 100 mm. (F) XZ two-photon excitation using raster (top) and conformal (bottom) scan patterns with preferential energy delivery to the deepest regions in the conformal scheme.A C B

JFig. 3 .

 3 Fig. 3. Automated reconstruction of the cell lineage tree for a cohort of zebrafish early blastulas. (A) SHG spot detection (red) overlaid to XY projections of SHG images; A1 to A5 encompass cell cycle 5 (with cell division 1at45min);seemovieS8.(B and C) THG image segmentation at the 8-cell stage (B) and 128-cell stage (C); see movie S9. (D and E)Reconstructionof4cell stage (D) and 256-cell stage (E) embryos. In white, raw THG signal; orange, detected SHG spots; colors, segmented cell contours. (F to J)Digi talembr yo from 1-cell stage to 512-cell stage. See movie S10. (K)Flatrepresentationof the lineage tree with eight-cell stage clones; same colors as in (G) to (J). (L) Spatial deployment of the cell lineage up to cycle 5. (Top) Animal pole view. (Bottom) Lateral views. See movie S11 and fig. S6. (M) Angle between successive sister cells dipoles, in degrees. 55°(black line) corresponds to random orientation. (N) Angle between sister cells dipoles and the normal to the embryo surface. Cells divided tangentially to the embryo surface until cycle 5 (16 to 32 cells). At cycle 5, the four central blastomeres [shown in black in (L) and (N)] divided orthogonally to the surface.

Fig. 4 .

 4 Fig. 4. Data statistical homogeneity, local features, and global patterns. (A) Mean cell cycle duration (metaphase-to-metaphase) from cycles 2 to 10 for six embryos (black line, average). (B) Mean cell cycle duration (green) compared with cycle duration at the beginning (blue) and at the end (red) of the division cycle, averaged for six embryos. (C) Temporal histogram of mitoses in six embryos. (D) Spatiotemporal representation of the mitotic pseudo-wave at cycle 8; time ranges from 0 min (blue, cycle onset) to 8 min (red, cycle end); white arrow points to the pseudo-wave origin. (E) Scatter plot of mitosis time (SHG peak at metaphase) as a function of distance to the pseudo-wave origin. (F) XY projection of the position of the pseudo-wave origin for cycles 5 to 10 with respect to the animal pole (red cross) and other mitoses at cycle 10 (black crosses). Data for five other embryos are in the supporting online material.


Karol Mikula

Róbert Špir

Michal Smíšek

Emmanuel Faure

Nadine Peyriéras

Nonlinear PDE based numerical methods for cell tracking in zebrafish embryogenesis

The paper presents numerical algorithms leading to an automated cell tracking and reconstruction of the cell lineage tree during the first hours of animal embryogenesis. We present results obtained for large-scale 3D+time two-photon laser scanning microscopy images of early stages of zebrafish (Danio rerio) embryo development. Our approach consists of three basic steps -the image filtering, the cell centers detection and the cell trajectories extraction yielding the lineage tree reconstruction. In all three steps we use nonlinear partial differential equations. For the filtering the geodesic mean curvature flow in level set formulation is used, for the cell center detection the motion of level sets by a constant speed regularized by mean curvature flow is used and the solution of the eikonal equation is essential for the cell trajectories extraction. The core of our new tracking method is an original approach to cell trajectories extraction based on finding a continuous centered paths inside the spatio-temporal tree structures representing cell movement and divisions. Such paths are found by using a suitably designed distance function from cell centers detected in all time steps of the 3D+time image sequence and by a backtracking in the steepest descent direction of a potential field based on this distance function. We also present efficient and naturally parallelizable discretizations of the aforementioned nonlinear PDEs and discuss properties and results of our new tracking method on artificial and real 4D data.

Introduction

The development of the modern microscopy techniques allows the in vivo imaging of organisms at cell level at very early stages of development without corrupting the cell integrity and normal development of the embryo. The multi-photon laser scanning microscopy is able to deliver 3D+time images of long periods of the embryonic development with a relatively short time step. With such data, one can perform various analyses of the embryogenesis and the development of organisms. Due to a similarity with human in many aspects and due to transparency for the microscope, the zebrafish (Danio rerio) embryogenesis is studied extensively and the results are used both in basic and applied biology and medicine research, e.g. in the anticancer drug design.

Although a large amount of work has been already done, e.g. by a combination of various image processing techniques and manual inspection, for the zebrafish developmental stages up to about one thousand cells [START_REF] Mikut | Automated processing of zebrafish imaging data: as u r v e[END_REF], there is a great challenge to study embryogenesis at very complex stages of development with thousands of cells present. In [START_REF] Krivá | 3D early embryogenesis image filtering by nonlinear partial differential equations[END_REF][START_REF] Zanella | Cells segmentation from 3-D confocal images of early zebrafish embryogenesis[END_REF][START_REF] Mikula | Segmentation of 3D cell membrane images by PDE methods and its applications[END_REF][START_REF] Bourgine | 4D embryogenesis image analysis using PDE methods of image processing[END_REF] new efficient and robust 3D filtering and segmentation algorithms together with a workflow for performing various analyses of such complex developmental stages were suggested and studied. Among other goals, the image processing algorithms aim to deliver the so-called cell lineage tree, the spatio-temporal branching process giving topological description of movement and division of the cells. This problem is related to tracking of cells during the embryogenesis. Having the tree, one can track the cells by following the links in the tree, either forward or backward. On the other hand, having the good cell tracking algorithm and thus individual cell trajectories, one can reconstruct the tree.

The tracking of cells and construction of the cell lineage tree for such complex stages of embryogenesis is a very difficult problem unsolved satisfactorily yet. A seminal work towards building the cell lineage tree for the complex stages of zebrafish embryo development based on stochastic simulated annealing minimization of a heuristic energy functional has been done in [START_REF] Melani | Algoritmos de procesamiento de imagenes para la reconstruccion del desarrollo embrionario del pez cebra[END_REF]. After the construction of the tree, the individual cells or cell populations are tracked. In this paper we suggest new PDE based cell tracking method for such complex stages of zebrafish embryogenesis. In the present work we continue the effort started in [START_REF] Bellaïche | New robust algorithm for tracking cells in videos of drosophila morphogenesis based on finding an ideal path in segmented spatio-temporal cellular structures[END_REF][START_REF] Mikula | 4D numerical schemes for cell image segmentation and tracking[END_REF] devoted to 2D+time and 3D+time cells tracking. In contrast to [START_REF] Melani | Algoritmos de procesamiento de imagenes para la reconstruccion del desarrollo embrionario del pez cebra[END_REF], in our approach we first extract all possible cell trajectories inside the 3D+time data and then the cell lineage tree is reconstructed by finding merging trajectories going backward in time.

Our overall approach to cell tracking consists of three basic steps -the image filtering, the cell centers detection and the cell trajectories extraction yielding the cell lineage tree reconstruction. In all three steps we use suitably designed nonlinear partial differential equations -for filtering the geodesic mean curvature flow in level set formulation is used, for cell center detection the motion of level sets by a constant speed regularized by mean curvature flow is used and solution of the eikonal equation is essential for the cell trajectories extraction. We also present efficient and naturally parallelizable discretizations of the aforementioned nonlinear PDEs and discuss properties and results of our new tracking method.

The first step of our technique is the image filtering. This is an important step because noise is always present in microscopy images and its level increases with the speed of image acquisition. We perform the filtering of 3D images of the 3D+time image sequence by the numerical solution of the so-called geodesic mean curvature flow (GMCF) model [START_REF] Caselles | Geodesic active contours[END_REF][START_REF] Kichenassamy | Conformal curvature flows: from phase transitions to active vision[END_REF][START_REF] Chen | Image denoising and segmentation via nonlinear diffusion[END_REF] which was chosen from several available methods by careful testing [START_REF] Krivá | 3D early embryogenesis image filtering by nonlinear partial differential equations[END_REF].

In the second step of our approach we use the level set center detection (LSCD) method [START_REF] Frolkovič | A counting number of cells and cell segmentation using advection-diffusion equations[END_REF] to extract an approximate position of the cell nuclei centers (which we also call the cell identifiers) in every 3D image of the 3D+time image sequence. LSCD method is represented by an evolutionary process based on the motion of level sets of 3D image intensity in normal direction by a constant speed regularized by the mean curvature flow. Following the local maxima of the solution during the evolution one can obtain a reasonable approximation of cell nuclei center positions and its number in every 3D image.

The cell identifiers represent the basic input for the third step of the overall tracking procedure, the cell trajectories extraction and the cell lineage tree reconstruction. The cell trajectories are extracted by finding centered paths inside a specifically constructed 4D spatio-temporal tree structure. Although we call it tree structure, it represents a subset of nonzero measure of the Euclidean space R 4 and is given by a 4D segmentation step. Thus, in the first step of the trajectories extraction we construct a 4D segmentation yielding the desired 4D spatio-temporal tree structure. Then a computation of constrained distance functions inside this 4D segmentation is performed by solving a spatially 4D eikonal equation. By a proper combination of computed distance functions we build a potential field which is backtracked in the steepest descent direction in order to get the cell trajectories. Finally we construct the cell lineage tree by detecting trajectories which merge together when going backward in time indicating mitosis and thus a branching node of the cell lineage tree.

The paper is organized as follows. In the next subsection we present the zebrafish embryo development data acquired by the two-photon laser scanning microscope. Then we present all three steps of our approach, the image filtering, the cell center detection and the cell trajectories extraction. Together with the PDE based models, we outline the numerical approaches used in all three steps of our tracking method. Finally we discuss the numerical experiment devoted to tracking cells in a representative 3D video of the early zebrafish embryogenesis.

3D+time = 4D image data of zebrafish embryogenesis

Data for processing comes from two-photon laser scanning microscopy and they are of general interest, so we present some illustrative examples how it looks like. The data we deal with represents first hours of zebrafish embryo development, approximately from the 4th until 14th-20th hour. Usually, the data comes in two channels (two colors), in one channel there is an acquisition of the cell nuclei and in the second channel the acquisition of the cell membranes. The embryo labeling is obtained through RNA injection performed at the one-cell stage to obtain expression of fluorescence proteins staining nuclei and membranes, respectively. The embryo imaging starts about 4 hours post-fertilization and take next 10-15 hours. The 3D images of both channels are obtained by moving the focal plane from the top more deeply inside the embryo and their quality depend on the speed of scanning in one plane. Many various quality datasets are available thanks to EC Embryomics and BioEmergences projects (http :/ /bioemergences .iscpif .fr). The data quality is related to a level of noise depending mainly on the size of the time step between acquisition of subsequent 3D images. The acquisition step ranges from 50 seconds to 5 minutes. A longer time step produces better 3D image quality, one has a good visual impression of what is happening during embryogenesis and such data is well suited for segmentation purposes e.g. for obtaining a shape of cells and other their characteristics during the embryogenesis [START_REF] Zanella | Cells segmentation from 3-D confocal images of early zebrafish embryogenesis[END_REF][START_REF] Mikula | Segmentation of 3D cell membrane images by PDE methods and its applications[END_REF]. On the other hand, such data is not suitable for tracking since the cells move too far between single 3D images and consequently mother-daughters cell correspondences can be lost.

In Fig. 1 we see the cell membranes visualization in the beginning and at the end of the imaging. The usage of membrane images for the cell tracking would be very natural, cf. [START_REF] Bellaïche | New robust algorithm for tracking cells in videos of drosophila morphogenesis based on finding an ideal path in segmented spatio-temporal cellular structures[END_REF], but the quality of currently produced 3D membrane images is not sufficient for our purposes. In Fig. 2 we plot the cell nuclei visualization, again at the beginning and at the end of the imaging. Although nuclei images contain certain level of noise they are suitable for application of our tracking technique. In both figures we can observe how complex the development of organism in few first hours of embryogenesis is. It grows from an unorganized set of cells to various morphogenetic structures corresponding to future organs of the zebrafish adult.

For the testing of the method and for the extraction of cell trajectories we use the representative dataset with acquisition step dθ = 67 seconds, N θ number of time steps and with dimension of every 3D image 512 × 512 × 104 voxels. The real voxel size is dx 1 = dx 2 = dx 3 = 1.37 micrometer in every spatial dimension and the whole 3D+time sequence has 760 time steps. In this paper we present its processing up to time step N θ = 479 in which the biologists selected manually cells forming presumptive organs during the zebrafish brain early embryogenesis. Using developed approach we can track those cell populations, see Section 5, and thus follow their dynamics and clonal history. Since the size of the single 3D image is about 27 MB, for 479 time steps we deal with about 13 GB of raw data which requires the usage of high-performance parallel computing facilities especially when designing spatially 4D numerical algorithms.

From a mathematical point of view, the 3D+time image sequence is understood as a function u(x 1 , x 2 , x 3 , θ), u : Λ → [0, 1], where Λ is a bounded spatio-temporal (rectangular) subdomain of R 4 . When we consider processing of 3D images, in Sections 2 and 3, we consider u at a corresponding time step θ , but for simplicity we denote it as u(x) = u(x 1 , x 2 , x 3 ) omitting the θ variable. Similarly, the subset of Λ for a given θ in 3D image processing will be denoted by Ω. Since our image processing algorithms are based on evolutionary PDEs we add the artificial time (filtering "time", center detection "time", "time" relaxation for the eikonal equation) variable t. This variable should not be confused with the real image sequence time θ .

Nonlinear diffusion filtering

The image filtering is the first step of our overall tracking procedure. The noise is intrinsically linked to the microscopy technique and its level increases with decreasing the time step dθ of the scanning. The careful study of various filtering methods was performed in [START_REF] Krivá | 3D early embryogenesis image filtering by nonlinear partial differential equations[END_REF] and the geodesic mean curvature flow (GMCF) was chosen as an optimal one with respect to smoothing of spurious noisy structures and sharpening the edges of 3D+time cell nuclei images during zebrafish embryogenesis. The GMCF filtering model is based on a discretization of the following nonlinear diffusion equation

u t -|∇u|∇. g |∇G σ * u| ∇u |∇u| = 0, (1)
where u(t, x), t > 0, represents the filtered image intensity function. We start from the initial condition u(0, x) = u 0 (x),

where u 0 is an original 3D image and we consider few discrete time steps of the discretized model. We consider zero

Neumann boundary conditions on the boundary ∂Ω of the 3D image domain Ω. In this model, the mean curvature motion of the level sets of function u is determined by the edge indicator function

g(s) = 1 1 + Ks 2 , K ≥ 0 ( 2
)
that is applied to the image gradient presmoothed by the Gaussian kernel G σ with a small variance σ . The essential property of this function is that its negative gradient points towards the edges in the image and overall nonlinear diffusion process given by (1) causes accumulation of the level sets of u along the boundaries of objects in the image and therefore the filtering is also edge preserving. The filtered image u F is obtained as the solution of (1) at a time t = T F . The optimal choice of the model and discretization parameters was studied in [START_REF] Krivá | 3D early embryogenesis image filtering by nonlinear partial differential equations[END_REF] on the basis of the mean Hausdorff distance of the level sets of the filtered image and the manually segmented nucleus shapes. In Fig. 3 we present the filtering result by using optimal parameters of the GMCF model.

Discretization of the nonlinear diffusion filtering equation

For discretization of (1) in time we use the semi-implicit approach that guarantees unconditional stability of the numerical scheme. Let us suppose that we solve the PDE (1) in the interval [0, T F ] and let N F be a number of uniform time steps with size τ F = T F /N F . The semi-implicit scheme for solving Eq. ( 1) is given by

u n -u n-1 τ F -∇u n-1 ∇. g ∇G σ * u n-1 ∇u n |∇u n-1 | = 0, (3)
where u n represents the solution in the nth filtering step, n = 1 ...N F .

For spatial discretization of (3) we use the finite volume method. First, we identify the finite volumes of the mesh T h with the voxels of the 3D image. We denote each finite volume by V ijk , i = 1 ...N 1 , j = 1 ...N 2 , k = 1 ...N 3 . Such a finite volume grid is rectangular and let h 1 = dx 1 , h 2 = dx 2 , h 3 = dx 3 represent the size of the voxels in x 1 , x 2 , x 3 direction, respectively. Let m(V ijk ) denote the volume of V ijk and c ijk its barycenter. By u n ijk we denote the approximate value of u n in c ijk .

For all volumes V ijk , we define two index sets. First, let N ijk denote the set of all (p, q, r) such that p, q, r ∈{-1, 0, 1}, |p| +|q| +|r| = 1. Then, let P ijk represent the set of all (p, q, r) satisfying p, q, r ∈{-1, 0, 1}, |p| +|q| +|r| = 2. Let us first consider (p, q, r) ∈ N ijk . The line connecting the center of V ijk and the center of its neighbor V i+p, j+q,k+r is denoted by 3) by |∇u n-1 | and integrate it over every finite volume V ijk . After application of the Green theorem on the right hand side we get

V ijk 1 |∇u n-1 | u n -u n-1 τ F dx = N ijk e pqr ijk g ∇G σ * u n-1 ∇u n |∇u n-1 | ν pqr ijk dγ . ( 4
)
As one can see from Eq. ( 4) in order to proceed further we have to approximate |∇u n-1 | on voxel faces e pqr ijk and on voxel

V ijk itself. Similarly, the term g(|∇u n-1 σ |), u n-1 σ = G σ * u n-1
must be approximated on the voxel faces. Moreover, the norm of the gradient |∇u n-1 | appears in denominator. Due to Evans and Spruck regularization approach [START_REF] Evans | Motion of level sets by mean curvature[END_REF], such term is substituted by the regularized one ε 2 +|∇u n-1 | 2 , where ε is the regularization parameter, usually ε ≪ 1. In order to approximate (regularized) values of the gradient norm we define, cf. [START_REF] Mikula | Finite volume schemes for the generalized subjective surface equation in image segmentation[END_REF], values of u n-1 in the midpoints y pqr ijk of the voxel edges which are approximated for any (p, q, r) ∈ P ijk by

u pq0 ijk = 1 4 u n-1 ijk + u n-1 i+p, j,k + u n-1 i, j+q,k + u n-1 i+p, j+q,k , u p0r ijk = 1 4 u n-1 ijk + u n-1 i+p, j,k + u n-1 i, j,k+r + u n-1 i+p, j,k+r , u 0qr ijk = 1 4 u n-1 ijk + u n-1 i, j+q,k + u n-1 i, j,k+r + u n-1 i, j+q,k+r .
Let us denote by ∇ pqr u n-1 ijk the approximation of the gradient in the barycenter x pqr ijk of the face e pqr ijk , (p, q, r) ∈ N ijk , of the voxel V ijk . Using the above definitions, we can define gradient approximations on faces by

∇ p00 u n-1 ijk = p u n-1 i+p, j,k -u n-1 ijk /h 1 , u p10 ijk -u p,-1,0 ijk /h 2 , u p01 ijk -u p,0,-1 ijk /h 3 , ∇ 0q0 u n-1 ijk = u 1q0 ijk -u -1,q,0 ijk /h 1 , q u n-1 i, j+q,k -u n-1 ijk /h 2 , u 0q1 ijk -u 0,q,-1 ijk /h 3 , ∇ 00r u n-1 ijk = u 10r ijk -u -1,0,r ijk /h 1 , u 01r ijk -u 0,-1,r ijk /h 2 , r u n-1 i, j,k+r -u n-1 ijk /h 3
and then the terms (including regularization) which we need in the scheme will be denoted by

Q pqr;n-1 ε,ijk = ε 2 + ∇ pqr u n-1 ijk 2 , Q n-1 ε,ijk = ε 2 + 1 6 N ijk ∇ pqr u n-1 ijk 2 , g pqr;n-1 ijk = g ∇ pqr u n-1 σ ;ijk , (5)
where u n-1 σ ;ijk denotes the result of convolution of u n-1 with Gaussian kernel in voxel V ijk . If we approximate the left hand side of (4) by and using all further definitions in (5) we get discretization of (1) in the following form

V ijk 1 |∇u n-1 | u n -u n-1 τ F dx ≈ m(V ijk ) 1 Q n-1 ε,ijk u n ijk -u n-1 ijk τ F •
m(V ijk ) u n ijk -u n-1 ijk τ F = Q n-1 ε,ijk N ijk m(e pqr ijk ) m(b pqr ijk ) g pqr;n-1 ijk Q pqr;n-1 ε,ijk u n i+p, j+q,k+r -u n ijk . (6)
The discrete equation ( 6), together with boundary conditions (the zero Neumann boundary conditions are treated by the reflection of inner values to additional finite volumes along the boundary) represent the linear system of equations with unknowns

u n ijk , i = 1 ...N 1 , j = 1 ...N 2 , k = 1 ...N 3 .
The system is solved in each step n by the parallel Red-Black SOR iterative method implemented in the MPI framework. Since we start the iterative process using the result from the previous time step n -1, we have a good initial approximation and parallel SOR method is fast approaching the linear system solution.

Cell center detection

The second step of our approach is the detection of cell nuclei centers, we call them also cell identifiers. The cell center detection method is based on the fact that objects visible in the image can be seen as humps of relatively higher image intensity, bright nuclei on dark background is a typical example, see figures in Section 1.1. Any such hump is represented by certain image intensity level sets. The diameter of these level sets allows us to distinguish between significant objects, e.g. cell nuclei, and spurious inner structures which still remain after GMCF filtering. For cell nuclei, the diameter d is relatively large, 0 ≪ c 1 ≤ d ≤ c 2 , while the diameter of the spurious inner structures is much smaller, 0 < d ≪ c 1 . If the level sets are moving (advected) at a constant speed in the direction of the inner normal, the encompassed volume is decreasing and finally the hump disappears. It is also well-known that if the evolution of the level set depends on its mean curvature (such evolution represents an intrinsic diffusion of the isosurfaces) then the speed of shrinking tends to infinity as the diameter of the level set tends to zero. So our model is based on the fact that the level sets with small diameter corresponding to spurious structures disappear quickly (by the above mentioned nonlinear advection-diffusion mechanism) while level sets representing cell nuclei are observable in a much longer time scale, in Fig. 4 we present example of such evolution computed by our model. Since the motion of every level set is given by the normal velocity V = δ + µk where δ and µ are constants (model parameters) and k is the mean curvature, we formulate our level set center detection (LSCD) method in the form of following nonlinear advection-diffusion equation [START_REF] Frolkovič | A counting number of cells and cell segmentation using advection-diffusion equations[END_REF] u t + δ ∇u |∇u|

.∇u -µ|∇u|∇.

∇u |∇u| = 0 (7)
which is applied to the initial condition given by u F , the result of GMCF filtering. Again, we consider the zero Neumann boundary condition and the equation is solved in time interval [0, T C ]. Due to the shrinking and smoothing of all (real and spurious) structures in the evolutionary process represented by [START_REF] Corsaro | Semi-implicit co-volume method in 3D image segmentation[END_REF], we observe a decrease of the number of local maxima of the solution u as time proceeds. This decrease is fast in the beginning and much slower later and we stop this process when the speed of decrease is below a certain threshold. The cell centers are detected as points of local maxima of the function u at the stopping time T C . This is done for every 3D image of the sequence and usually 10 discrete time steps of LSCD model are sufficient to get high quality estimate of the cell identifiers. Moreover, if we still need to remove some false positives (more identifiers inside one nucleus) we can use the nucleus segmentation approach by generalized subjective surface (GSUBSURF) method [START_REF] Sarti | Subjective surfaces: a method for completing missing boundaries[END_REF][START_REF] Corsaro | Semi-implicit co-volume method in 3D image segmentation[END_REF][START_REF] Mikula | 3D embryogenesis image segmentation by the generalized subjective surface method using the finite volume technique[END_REF][START_REF] Zanella | Cells segmentation from 3-D confocal images of early zebrafish embryogenesis[END_REF] as described in [START_REF] Bourgine | 4D embryogenesis image analysis using PDE methods of image processing[END_REF], Section 6.3.1. The set of all cell identifiers detected in all time steps of the image sequence is denoted by s l m , m = 1, ..., n l C , l = 1, ..., N θ . In [START_REF] Bourgine | 4D embryogenesis image analysis using PDE methods of image processing[END_REF][START_REF] Zanella | Cells segmentation from 3-D confocal images of early zebrafish embryogenesis[END_REF][START_REF] Mikula | Segmentation of 3D cell membrane images by PDE methods and its applications[END_REF][START_REF] Campana | Methods toward in vivo measurement of zebrafish epithelial and deep cell proliferation[END_REF][START_REF] Melani | Algoritmos de procesamiento de imagenes para la reconstruccion del desarrollo embrionario del pez cebra[END_REF] the centers detected by LSCD were used for various purposes including cell shape, epithelium and whole embryo segmentation and tracking.

Discretization of the nonlinear level set center detection equation

The time discretization of LSCD equation ( 7) is explicit in advective and semi-implicit in diffusion parts and it uses the finite volume method together with upwind principle for space discretization [START_REF] Frolkovič | Flux-based level set method: a finite volume method for evolving interfaces[END_REF][START_REF] Frolkovič | A counting number of cells and cell segmentation using advection-diffusion equations[END_REF][START_REF] Mikula | Finite volume schemes for the generalized subjective surface equation in image segmentation[END_REF]. We chose a number of center detection steps N C and define a uniform time step size τ C = T C /N C . Let us define integrated advective fluxes for (p, q, r) ∈

N ijk by v pqr ijk = m e pqr ijk δ u n-1 i+p, j+q,k+r -u n-1 ijk Q pqr;n-1 ε,ijk m(b pqr ijk ) (8)
and identify inflow boundaries by the index set N in ijk := {(p, q, r) ∈ N ijk , v pqr ijk ≤ 0}. Using these definitions and definitions (5) we can write a discretization of Eq. ( 7) in the following form m(V ijk )

u n ijk -u n-1 ijk τ C = N in ijk v pqr ijk u n-1 i+p, j+q,k+r -u n-1 ijk + Q n-1 ε,ijk N ijk m e pqr ijk u n i+p, j+q,k+r -u n ijk Q pqr;n-1 ε,ijk m(b pqr ijk ) . ( 9
)
The method is parallelized in MPI framework and the Red-Black SOR parallel solver is used for solving the linear system [START_REF] Frolkovič | Flux-based level set method: a finite volume method for evolving interfaces[END_REF] in every time step n.

Cell trajectories extraction and lineage tree construction

The main new idea of our cell tracking method is presented in this section. The tracking is based on the cell trajectories extraction which is done by finding centered paths inside a specifically constructed 4D spatio-temporal tree structure. The 4D spatio-temporal tree structure is a subset of all doxels (4 dimensional pixels) of the 3D image sequence which represents the 3D cell shapes moving (and splitting) in time. The method is composed by five basic steps:

• construction of a 4D segmentation yielding the 4D spatio-temporal tree structure, • computation of constrained 4D distance functions inside this 4D segmentation and designing of a proper combination of them in order to build a potential field for tracking,

• extraction of cell trajectories using the steepest descent direction of the potential field,

• centering the extracted trajectories inside the 4D spatio-temporal trees in order to eliminate duplicates,

• construction of the cell lineage tree by detecting trajectories which merge together when going backward in time indicating mitosis and thus a branching node of the lineage tree.

4D segmentation

We call 4D segmentation a spatio-temporal structure which approximates the space-time movement of nuclei in 3D+time image sequence. As one can see from Fig. 5 the shape of cell nuclei in 3D images is reasonably approximated by spheres or ellipsoids. Thus, in order to construct the 4D segmentation we use cell identifiers detected in all time steps,

., n l

C , l = 1, ..., N θ , and create 4D ellipsoids around all these points. Currently we are using the same halfaxes in space equal to 2dx, where dx = dy = dz is voxel size (same in all three spatial dimensions) and with halfaxis equal to dθ in temporal dimension. The nonzero temporal halfaxis is important due to the time overlap which we create and thus we achieve better connectivity of 4D spatio-temporal tree structures, Fig. 6 shows our 4D segmentation projected to one time step. Thanks to the time overlap we also make connected such branches of the 4D spatio-temporal tree where a cell center is not detected in one frame but it is detected in two neighboring frames and thus we correct false negative errors of the cell center detection algorithm.

We also have to note that our current approach is not optimal in the sense of getting real 4D image segmentation of cell nuclei movement. Such approach using 4D generalized subjective surface (GSUBSURF) method was originally developed in [START_REF] Mikula | 4D numerical schemes for cell image segmentation and tracking[END_REF] but used only locally in time (for short time sequence) and space (for a detail of the full 3D volume) and utilization of such type of method for long 3D+time image sequences remains a great challenge from modeling, discretization and (parallel) computation point of view and it will be objective of our further study. On the other hand, the simplified approach suggested and used in this paper also gives, thanks to the time overlap, the reliable 4D segmentation which can be used in further steps of cell trajectories extraction.

Constrained distance functions

Our 4D segmentation containing 4D spatio-temporal tree structures can be represented by a 4D piecewise constant function, with some BIG value outside of the segmentation and with zero value inside it. Using this information, we compute two types of distance functions inside the 4D spatio-temporal tree structures. We call them constrained because all the calculations are constrained by the boundaries of the 4D segmentation. Due to that fact, the computed distances between 4D points of the 3D+time image sequence are not a standard Euclidean distances in R 4 but they represent a minimal Euclidean paths between the points inside the 4D segmentation.

The first type of distance functions will be denoted by D l (x 1 , x 2 , x 3 , θ), l ≥ 1. The function D l (x 1 , x 2 , x 3 , θ) represents the distance of any inner point (x 1 , x 2 , x 3 , θ) of the 4D segmentation to the set of cell identifiers {s l m , m = 1, ..., n l C } detected in the lth time step. Moreover we impose the constraints that D l (x 1 , x 2 , x 3 , θ) = BIG if θ<l, which means we look only forward in time when computing D l , and D l (x 1 , x 2 , x 3 , θ) = BIG for all (x 1 , x 2 , x 3 , θ) outside of the 4D segmentation. It is worth to note that we compute D l for all time steps l = 1, ..., N θ of the image sequence so it represents the most computationally demanding part of our approach. The second type of distance function, D B (x 1 , x 2 , x 3 , θ), is computed just once and it represents the distance of any inner point of the 4D segmentation to the boundary of the 4D segmentation.

In this case, D B (x 1 , x 2 , x 3 , θ) = 0f o r all (x 1 , x 2 , x 3 , θ) outside of the 4D segmentation. Both types of distance functions are computed numerically on the doxel structure of the 3D+time image sequence using the 4D Rouy-Tourin scheme. The scheme is presented in the next section.

Together with computing constrained distance functions D l (x 1 , x 2 , x 3 , θ), l ≥ 1, we subsequently build their combination denoted by D(x 1 , x 2 , x 3 , θ). The combined distance function D(x 1 , x 2 , x 3 , θ) is built in such a way that it gives the maximal shortest distance of a point (x 1 , x 2 , x 3 , θ) of the 4D segmentation to a cell identifier s l m among all m and all l ≤ θ . Simply saying, it represents the distance to the most far cell identifier to which the point (x 1 , x 2 , x 3 , θ) is continuously connected inside the 4D segmentation. Such construction of the combined distance function D will help us to get all, also partial, cell trajectories contained in our 4D spatio-temporal tree structures. Its utilization represents a main important novelty with respect to [START_REF] Bellaïche | New robust algorithm for tracking cells in videos of drosophila morphogenesis based on finding an ideal path in segmented spatio-temporal cellular structures[END_REF][START_REF] Mikula | 4D numerical schemes for cell image segmentation and tracking[END_REF] where just D = D 1 was used so only cell trajectories starting in the first frame could be detected. In order to construct the combined distance function D, first we define a set of values

D(x 1 , x 2 , x 3 ,θ)= D l (x 1 , x 2 , x 3 ,θ); l ∈{1, 2,...}, D l (x 1 , x 2 , x 3 ,θ)<BIG
and then we set

D(x 1 , x 2 , x 3 ,θ)= max D(x 1 , x 2 , x 3 ,θ), D(x 1 , x 2 , x 3 ,θ) = ∅ BIG otherwise.
In Fig. 7 we show how the combined distance function contains information about all simply connected components of the 4D segmentation. Inside the valleys encompassed by the BIG values, the value of D is growing from zero, in cell identifier where any simply connected component begins, up to a maximal value, where the simply connected component ends. Now, looking at Fig. 7 and also Fig. 8 left, we can think about D as a potential field and traverse it in the steepest descent direction from the highest value in every simply connected component to the zero value. The paths obtained in such way would represent good approximation of the space-time cell trajectories. And, if the 4D segmentation would contain only perfectly separated 4D spatio-temporal tree structures, we would obtain correctly all (also partial) cell trajectories which can be extracted from the data. Unfortunately, in the real 4D data it is not always the case and we must deal with imperfections given mainly by a cells overlapping. In order to overcome this difficulty we have to keep the extracted paths at a certain distance from the spatio-temporal cell boundaries or, in other words, they should be more centered inside the 4D segmentation. This can be achieved by using the constrained distance function D B (x 1 , x 2 , x 3 , θ) [15,1,17] values of which grow from boundaries to the center of the 4D spatio-temporal trees, see Fig. 8 right. Based on these facts we build the new potential field [START_REF] Frolkovič | A counting number of cells and cell segmentation using advection-diffusion equations[END_REF] which will be used in our algorithm for the extraction of cell trajectories.

V (x 1 , x 2 , x 3 ,θ)= D(x 1 , x 2 , x 3 ,θ)-D B (x 1 , x 2 , x 3 ,θ)

Numerical method for computing the constrained distance function

For the calculation of the constrained distance functions we solve numerically the so-called time relaxed eikonal equation d t +|∇d|=1 [START_REF] Kichenassamy | Conformal curvature flows: from phase transitions to active vision[END_REF] for the unknown function d(x 1 , x 2 , x 3 , θ, t). In contrast to Sections 2 and 3 we solve here spatially 4D problem, so ∇d is a 4D gradient of the function d, i.e. the vector of partial derivatives with respect to x 1 , x 2 , x 3 and θ variables. For the discretization of Eq. ( 11) we use spatially the 4D Rouy-Tourin scheme [START_REF] Rouy | Viscosity solutions approach to shape-from-shading[END_REF]. The constrained distance function, which we look for, is obtained as an equilibrium of the numerical solution, i.e. it satisfies numerically the classical eikonal equation |∇d| = 1. When solving Eq. ( 11), the values of the solution at points outside the 4D segmentation are fixed to a BIG value which is arbitrary but bigger than any distance which can be obtained for points inside the 4D segmentation. For Eq. ( 11), we have to prescribe zero Dirichlet condition on the set from which we compute the distances. In case of computing the constrained distance function D l we prescribe the zero values to the set of cell identifiers {s l m , m = 1, ..., n l C } detected in the lth time step. In case of computing the constrained distance function D B we prescribe the zero values for all boundary points of the 4D segmentation.

In analogy with Sections 2 and 3, where we solve spatially 3D problems, we identify here the 4D doxels with finite volumes V ijkl having four indices. Without losing generality, we rescale the time step dθ to be equal to dx 1 = dx 2 = dx 3 and denote their common value by h D (standardly we set h D = 1). Let d n ijkl denote the approximate value of solution d in barycenter of V ijkl in a discrete step t n = nτ D where τ D is the length of step discretizing t variable. Then, for every V ijkl we define the index set N ijkl of all (p, q, r, s) such that p, q, r, s ∈{-1, 0, 1}, |p| +|q| +|r| +|s| = 1. In order to build the scheme, for any (p, q, r, s) ∈ N ijkl , we define (see formula [START_REF] Frolkovič | A counting number of cells and cell segmentation using advection-diffusion equations[END_REF] in [START_REF] Rouy | Viscosity solutions approach to shape-from-shading[END_REF])

D pqrs ijkl = min d n-1 i+p, j+q,k+r,l+s -d n-1 ijkl , 0 2 (12)
and then also

M 1000 ijkl = max D -1,0,0,0 ijkl , D 1,0,0,0 ijkl , M 0100 ijkl = max D 0,-1,0,0 ijkl , D 0,1,0,0 ijk M 0010 ijkl = max D 0,0,-1,0 ijkl , D 0,0,1,0 ijkl , M 0001 ijkl = max D 0,0,0,-1 ijkl , D 0,0, 0,1 ijkl . (13)
Using these notations, the 4D Rouy-Tourin scheme for solving Eq. ( 11) has the following form

d n ijkl = d n-1 ijkl + τ D - τ D h D M 1000 ijkl + M 0100 ijkl + M 0010 ijkl + M 0001 ijkl . (14)
Due to stability reasons, the coupling τ D = h D /2i s used. We also speed-up the basic algorithm ( 14) by subsequently omitting calculations in points which are already in equilibrium using the fixing strategy proposed in [START_REF] Bourgine | Extraction of the intercellular skeleton from 2D microscope images of early embryogenesis[END_REF], see also [START_REF] Mikula | Segmentation of 3D cell membrane images by PDE methods and its applications[END_REF]. The method [START_REF] Mikula | Finite volume schemes for the generalized subjective surface equation in image segmentation[END_REF] with the fixing was parallelized using OpenMP programming interface for shared memory parallel servers. Such parallel implementation can be used for processing the long-time image sequences, e.g. the one with 479 3D volumes discussed in Section 5. In that case we needed almost 200 GB of shared memory and we obtained the speed-up 23.5 for 32 threads. Such speed-up was obtained thanks to the utilization of the NUMA library functions which allow the optimal NUMA nodes memory usage on the nowadays standard servers with NUMA (non-uniform memory access) architecture.

Extraction of the cell trajectories

The cell trajectory is represented by a series of points in space-time (discrete spatio-temporal curve) for which we prescribe the condition that there exists exactly one point P l of the trajectory in time step

l = N b , ..., N e , 1 ≤ N b < N e ≤ N θ .
The extraction of cell trajectories is realized in two steps:

• first, we use backtracking in time by the steepest descent direction of the potential V built in (10) starting from all cell identifiers s l m , m = 1, ..., n l C detected in all time steps l = 2, ..., N θ , • then, we center all the extracted trajectories inside the 4D spatio-temporal trees, by using the constrained distance function D B , in order to eliminate duplicates.

The first step is realized as follows: Let s e m be one of the cell centers detected in the eth time step. We search recursively in the nearest vicinity of s e m , but only in the current time step e and the previous time step e -1, for a doxel with the minimal value of the potential V strictly less than its value at s e m . The cell trajectory point P e for the time step e is defined as the (last in search) doxel from which we move to a point in the previous time step e -1. Then from this point we continue the same search, etc. We end the process when we cannot move from the time step b to the previous time step b -1b y decreasing the value of potential V . Then the last point of the search in time step b, which must be some detected center, becomes the first point of the cell trajectory starting in the time step b and ending in time step e (when going forward in time). As an output of this first step, we get as many trajectories as is the number of cell identifiers in all time steps except the first one, which means that we have

N θ l=2 n l C trajectories.
The first step allows us to extract all cell trajectories starting and ending in any time step l between l = 1 and l = N θ .

On the other hand, there exist many duplicates, i.e. trajectories which represent the same cell space-time movement but for shorter time. To illustrate the above fact, let us consider a long cell trajectory going from the first to the last 3D volume of the image sequence inside a branch of the 4D spatio-temporal tree structure. Since we start the descent from centers detected at every l = N θ , ..., 2, and they all lay inside the branch corresponding to that cell, we obtain N θ -1( s h o r t e r ) trajectories laying inside the same branch of the 4D spatio-temporal tree. These trajectories can slightly differ because the steepest descent search does not give necessarily the same trajectory points when starting in different time steps. However, we know that the distance function D B is growing from the boundary to the centerline of the branch. If we restrict D B to a particular 3D volume, it grows from boundary of the branch to exactly one maximum located in the detected cell center. This holds generically due to the construction of our 4D segmentation by interconnecting the spheres around the detected cell centers. Thus, if we move a point of trajectory which is not centered using the steepest ascend direction of function D B iside the 3D volume, we end up in the cell center. After this procedure applied to all trajectories, the points of any shorter trajectory become a subset of points of the longer trajectory and we can remove the shorter one just by comparing the trajectories points. At the end, we obtain a set of unique trajectories in the sense that a mother cell life is represented by as many trajectories as is the number of her descendant cells. By experiments presented in Figs. 9-11 we illustrate the behavior of our cell trajectories extraction method and test its implementation.

First, we created a 4D dataset (3D volumes with 100 × 100 × 100 voxels in 40 time steps) containing few cells that are moving in space and time and performing several cell divisions. The dataset contains also partial cell movements starting and ending inside the 4D data. In fact, here we created directly the 4D segmentation, described in Section 4.1, which gives us the 4D spatio-temporal tree structures corresponding to cells movements and divisions. The 3D projection of that 4D segmentation is given in Fig. 9. Then we performed computations of constrained distance functions D l , l = 1, ..., 40 together with building of the combined distance function D and computation of the constrained distance function D B . After, we built the potential V and performed both steps (the steepest descent of V and centering by D B ) of cell trajectories extraction. The resulting extracted trajectories are presented in Fig. 10. By different colors we plot 3D projection of 4D points (represented by small 4D spheres) of different trajectories. As we can see, all trajectories, including partial and dividing, are correctly extracted. Since the trajectories are overlapping we can see only the color of the last drawn one, e.g. the red trajectory is the same as the blue one before the first division.

In the second testing example, we created two touching branches of a 4D segmentation imitating situation that the cell nuclei are not perfectly separated due to acquisition errors. And, we tested our method in an extremal case when the branches are fully overlapped in space in two subsequent time steps (we used the same cell identifiers for those two time steps when building the 4D segmentation). In Fig. 11 we present the result obtained after the first step of our method, i.e. after the backtracking in time by the steepest descent direction of V . Since we perform only the steepest descent and not the centering in this step, we can obtain two separated paths. When we center them in the second step they will remain different (although with two equal points) and thus they will represent two correctly extracted cell trajectories.

Cell lineage tree reconstruction

Let us consider a mother cell which is going to divide at time step l, i.e. at time step l + 1i t has two descendants and in later times maybe more due to further divisions. Without losing generality, let the number of descendant of this cell in the whole image sequence be N d = 2 m . As noticed already in the previous section, up to time l the life of the cell is represented by exactly N d trajectories which are, however, until time l composed by the same spatio-temporal points. From the time l + 1, the half of trajectories differs from the second half, but every half is again composed by the same points until the next division. The representation of cell life by the multiple trajectories which are partially the same does not cause any problem in visualization and/or reconstruction of single cell or cell population movements and divisions. However, the reduction of the equal multiple parts of trajectories is necessary for the reconstruction of the cell lineage tree and is explained below. The cell lineage tree is stored in a structure with a specific format. For each node of the cell lineage tree we store the spatio-temporal coordinates of the corresponding point, its unique global ID and the global ID of its mother in the extracted trajectory. In case of the first point of the trajectory, the mother ID is set to zero. To fill in the lineage tree structure, we subsequently take all points of the extracted trajectories, starting with the longest trajectories, and add subsequently the node representation of those points to the structure. For every trajectory we start by the first point and check whether the node corresponding to this point already exists in the structure. If it exists, it means that we have already added a trajectory which has some part equal to the current one. And also that there exists a later time after which the trajectories differ. We skip all equal points which are already represented as nodes in the lineage tree structure. Only the first different point of the current trajectory is added to the structure, together with the mother ID of the last equal point from the previously added trajectory. Using such approach we obtain the whole cell lineage tree where each node exists only once and the nodes are logically linked together in the same manner as it is in the real cell mother-daughter relation. For illustration, in Fig. 12 we present a detail of the cell lineage tree reconstructed for a data set discussed in the next section. The visualization of the cell lineage tree is done by MovIT, the software built by Thiery Savy, École Polytechnique, Paris.

Numerical experiment on real zebrafish embryogenesis data

We performed the cell tracking by the approach designed in Sections 2-4 on the dataset described in Section 1.1 representing the zebrafish brain early embryogenesis. All 3D images were filtered by 10 steps of GMCF model and the cell nuclei identifiers were detected by 10 steps of LSCD algorithm. From several millions of cell identifiers we built the 4D segmentation and then the cell trajectories were extracted by the approach developed in Section 4. We present here Figs. [START_REF] Melani | Algoritmos de procesamiento de imagenes para la reconstruccion del desarrollo embrionario del pez cebra[END_REF][START_REF] Mikula | Finite volume schemes for the generalized subjective surface equation in image segmentation[END_REF][START_REF] Mikula | 3D curve evolution algorithm with tangential redistribution for a fully automatic finding of an ideal camera path in virtual colonoscopy[END_REF][START_REF] Mikula | 3D embryogenesis image segmentation by the generalized subjective surface method using the finite volume technique[END_REF] showing results of the tracking procedure. In Fig. 13 one can see the visualization of the cell trajectories in developing embryo. To that goal we built software, running on graphics card, where we can fluently zoom, rotate and animate in time a 3D scene even with a very high number of trajectories. The trajectories are displayed as short lines (in random colors) connecting few subsequent spatio-temporal points with a freely chosen starting time.

In Fig. 14 we use a similar visualization, but we show trajectories only for a subset of all cells and different colors represent different cell populations trajectories. The subset represents seven cell populations forming different presumptive organs in the zebrafish embryo brain, e.g. hypoblast including pre-chordal-plate and notochord (red), presumptive hypothalamus (yellow), ventral telencephalon (brown), right eye (blue) or left eye (green). By the colored ellipsoid we also plot a motion of the mean trajectory of every cell population. One can see the evolution of cells from the 0th still chaotic stage, through 100th, 200th, 300th, 400th time steps where the cells are becoming more compactly localized up to 479th time step where the cells were marked by biologists. In order to have more dynamical impression, in Fig. 14 we show long parts of cell trajectories to illustrate the cell population movements from the central part up to the presumptive organ stage.

The last series of plots, in Fig. 16, presents the tracking of a single cell in several subsequent time steps of embryogenesis. One can see that the cell moves in a complicated manner, however, the visual inspection in 3D+time data shows that its trajectory is correct.

At the end of this section we shortly discuss an important difference between our new approach and previously suggested approaches [START_REF] Bellaïche | New robust algorithm for tracking cells in videos of drosophila morphogenesis based on finding an ideal path in segmented spatio-temporal cellular structures[END_REF][START_REF] Mikula | 4D numerical schemes for cell image segmentation and tracking[END_REF]. In [START_REF] Bellaïche | New robust algorithm for tracking cells in videos of drosophila morphogenesis based on finding an ideal path in segmented spatio-temporal cellular structures[END_REF][START_REF] Mikula | 4D numerical schemes for cell image segmentation and tracking[END_REF], only the distance function D 1 was used instead of combined distance function D. It means that we can extract only cell trajectories laying in simply connected components of 4D segmentation starting in the first time step of the image sequence. After the centering step such trajectories are composed by a subset of all detected cell centers. For real data experiment presented in this section, such extracted trajectories would cover only 52% of all cell centers, which means we would lose many partial cell trajectories which can be extracted from the data by the new approach. In fact, by using the method presented in this paper, the extracted cell trajectories cover all detected cell centers, which is an important improvement obtained by using the combined distance function D.

Conclusions

In this paper we developed a new approach to the cell tracking in 3D+time confocal microscopy data. We applied the tracking method to complex stages of the zebrafish early embryogenesis images. Together with [START_REF] Melani | Algoritmos de procesamiento de imagenes para la reconstruccion del desarrollo embrionario del pez cebra[END_REF] it is a first cell tracking method able to deal with such type of data. We discussed mathematical models in form of nonlinear PDEs the tracking method is based on and we also present their numerical discretizations which leads to parallel computing algorithms.

Fig. 1 .

 1 Fig. 1. Membrane images. Left: volume rendering of data in the first 3D image (4 hours after fertilization), middle: 2D slice of that raw data, right: 2D slice of the raw data in 20th hour of embryo development.

Fig. 2 .

 2 Fig. 2. Nuclei images. Top left: volume rendering of data in the first 3D image (4 hours after fertilization), top right: 2D slice of that raw data, bottom left: volume rendering of the data in the 3D image in 20th hour of embryo development, bottom right: 2D slice of that raw data.

Fig. 3 .

 3 Fig. 3. Image filtering. Left: isosurface visualization of the original 3D data in 100th time step, right: the result after 10 filtering steps by GMCF.

  length m(b pqr ijk ). The faces of the finite volume V ijk are denoted by e pqr ijk with area m(e pqr ijk ) and normal ν pqr ijk . Let x pqr ijk be the point where the line b pqr ijk crosses the face e pqr ijk . Finally, for any (p, q, r) ∈ P ijk , let y pqr ijk denote the midpoints of the voxel edges. The approximate value of u n-1 at x pqr ijk and y pqr ijk , where (p, q, r) belongs to the corresponding index set, is denoted by u pqr ijk , omitting the time index, as only the values from the time level n -1w i l l be needed in these points. Now let us divide Eq. (

Fig. 4 .

 4 Fig. 4. Cell center detection for the lth 3D image, l = 100. We plot a detail of the solution u(x 1 , x 2 , x 3 ), x 3 = 60h 3 in the 0th (left), 10th (middle) and 30th (right) time step of LSCD evolution.

Fig. 5 .

 5 Fig. 5.Examples of segmentation of cell nuclei using 3D GSUBSURF method[START_REF] Mikula | Segmentation of 3D cell membrane images by PDE methods and its applications[END_REF].

Fig. 6 .

 6 Fig. 6. Projection of the 4D segmentation to one time step. For each cell identifier we can see central sphere (inside the red circle) surrounded by the overlaps from the previous time step (inside the green circle) and from the next time step (inside the blue circle). Three further parts of 4D segmentation are shown as well. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) s l m , m = 1, ..., n lC , l = 1, ..., N θ , and create 4D ellipsoids around all these points. Currently we are using the same halfaxes in space equal to 2dx, where dx = dy = dz is voxel size (same in all three spatial dimensions) and with halfaxis equal to dθ

Fig. 7 .Fig. 8 .

 78 Fig. 7. On the left we can see the plot of the constrained distance function D 1 , in the middle we can see the plot of the constrained distance function D l , l > 1a n d on the right we can see the plot of the combined distance function D after time step l. Information about the (most left) simply connected component of the 4D segmentation which was not detected in the 1st time step is added into the function D in time step l. The values of D in the previously detected simply connected components were not changed by values of D l .

Fig. 9 .

 9 Fig. 9. 3D projection of the 4D segmentation of the artificial dataset where we can see cells moving and dividing inside the 4D volume.

Fig. 10 .

 10 Fig. 10. 3D projection of the 4D points of the extracted cell trajectories. (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)

Fig. 11 .

 11 Fig. 11. Separated trajectories after tracking, before post processing centering.

Fig. 12 .

 12 Fig. 12. A detail of the extracted cell lineage tree.

Fig. 13 .

 13 Fig. 13. Visualization of the zebrafish embryo cell trajectories in time steps 0, 100, 200, 300, 400 and 479.

Fig. 14 .

 14 Fig. 14. Visualization of the formation of presumptive organs during embryogenesis obtained by backward tracking of cells selected by biologists in time step 479. The visualization presents seven cell populations in time steps 0, 100, 200, 300, 400 and 479. (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)

Fig. 15 .

 15 Fig. 15. The cell population movement visualized by long parts of trajectories. Top left: trajectories from time step 0 to time step 120, top right: from 120 to 240, bottom left: from 240 to 360, bottom right: from 360 to 479.

Fig. 16 .

 16 Fig. 16. Tracking of a single cell (green cross) displayed together with 2D slice of 3D data where the current cell identifier is localized. The visualization is presented for time steps 0, 100, 200 and 300. On the left: the 2D slice of current cell location, on the right: detail of the 2D slice with tracked cell better visible. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Acknowledgement

We thank to all former members of EC Embryomics and BioEmergences project teams, particularly to Paul Bourgine whose continuous encouragement has helped to obtain presented results.

The authors were supported by the grant APVV-0184-10.