A workflow to process 3D+time microscopy images of developing organisms and reconstruct their cell lineage

Abstract

The quantitative and systematic analysis of embryonic cell dynamics from in vivo 3D+time image data sets is a major challenge at the forefront of developmental biology. Despite recent breakthroughs in the microscopy imaging of living systems, producing an accurate cell lineage tree for any developing organism remains a difficult task. We present here the BioEmergences workflow integrating all reconstruction steps from image acquisition and processing to the interactive visualization of reconstructed data. Original mathematical methods and algorithms underlie image filtering, nucleus centre detection, nucleus and membrane segmentation, and cell tracking. They are demonstrated on zebrafish, ascidian and sea urchin embryos with stained nuclei and membranes. Subsequent validation and annotations are carried out using Mov-IT, a custom-made graphical interface. Compared with eight other software tools, our workflow achieved the best lineage score. Delivered in standalone or web service mode, BioEmergences and Mov-IT offer a unique set of tools for in silico experimental embryology.

ARTICLE

Received 28 Apr 2015 | Accepted 18 Sep 2015 | Published 25 Feb 2016

DOI: 10.1038/ncomms9674 OPEN

A workow to process 3Dtime microscopy

images of developing organisms and reconstruct their cell lineage

Emmanuel Faure1,2,3,w,*, Thierry Savy1,2,3,4,*, Barbara Rizzi1,3,5,w,*, Camilo Melani1,4,5,w,*, Olga Staov6,*,

Dimitri Fabrges1,3,4,*, Rbertpir6, Mark Hammons1,3,4, Rbertnderlk6, Galle Recher1,3,

Benot Lombardot1,2,3, Louise Duloquin1,3, Ingrid Colin3, Jozef Kollr6, Sophie Desnoulez3, Pierre Affaticati7,w, Benot Maury3, Adeline Boyreau3,4, Jean-Yves Nief8, Pascal Calvat8, Philippe Vernier7,w, Monique Frain3,4, Georges Lutfalla9, Yannick Kergosien1,3,4,10, Pierre Suret11, Mariana Remekov6,**, Ren Doursat1,2,3,4,w,**,

Alessandro Sarti5,w,**, Karol Mikula6,**, Nadine Peyriras1,3,4,** & Paul Bourgine1,2,3,4,**

The quantitative and systematic analysis of embryonic cell dynamics from in vivo 3Dtime

image data sets is a major challenge at the forefront of developmental biology. Despite recent breakthroughs in the microscopy imaging of living systems, producing an accurate cell lineage tree for any developing organism remains a difcult task. We present here the BioEmergences workow integrating all reconstruction steps from image acquisition and processing to the interactive visualization of reconstructed data. Original mathematical methods and algorithms underlie image ltering, nucleus centre detection, nucleus and membrane segmentation, and cell tracking. They are demonstrated on zebrash, ascidian and sea urchin embryos with stained nuclei and membranes. Subsequent validation and annotations are carried out using Mov-IT, a custom-made graphical interface. Compared with eight other software tools, our workow achieved the best lineage score. Delivered in standalone or web service mode, BioEmergences and Mov-IT offer a unique set of tools for in silico experimental embryology.

1 Complex Systems Institute Paris Ile-de-France (ISC-PIF, UPS3611), CNRS, 75013 Paris, France. 2 Research Center in Applied Epistemology (CREA, UMR7656), CNRS and Ecole Polytechnique, 75005 Paris, France. 3 Multiscale Dynamics in Animal Morphogenesis (MDAM), Neurobiology & Development (N&D, UPR3294), CNRS, 91198 Gif-sur-Yvette, France. 4 BioEmergences Laboratory (USR3695), CNRS, Universit Paris-Saclay, 91198 Gif-sur-Yvette, France.

5 Department of Electronics, Information and Systems, University of Bologna, 40126, Italy. 6 Department of Mathematics and Descriptive Geometry, Slovak University of Technology, 81005 Bratislava, Slovakia. 7 Neurobiology & Development (N&D, UPR3294), CNRS, 91198 Gif-sur-Yvette, France. 8 Computing Center of the National Institute for Nuclear Physics and Particle Physics (CC-IN2P3, USR6402), CNRS, 69100 Villeurbanne, France. 9 Dynamics of Membrane Interactions in Normal and Pathological Cells (DIMNP, UMR5235), CNRS and Universit Montpellier 2, 34090 Montpellier, France. 10 Medical Informatics and Knowledge Engineering in e-Health (LIMICS, UMR1142), CNRS and Universit Paris 13, 93017 Bobigny, France. 11 Laser, Atomic and Molecular Physics Laboratory (UMR8523), CNRS and Universit Lille 1-Science and Technology, 59650 Villeneuve-dAscq, France. w Present addresses: IRIT (UMR5505), CNRS, Universit de Toulouse, 31400 Toulouse, France (E.F.); TEFOR Core Facility, Paris-Saclay Institute of Neuroscience (UMR9197), CNRS, Universit Paris-Sud, 91190 Gif-sur-Yvette, France (B.R., P.A.); Department of Computer Science, University of Buenos Aires, 1053 Buenos Aires, Argentina (C.M.); Paris-Saclay Institute of Neuroscience (UMR9197), CNRS, Universit Paris-Sud, 91190 Gif-sur-Yvette, France (P.V.); School of Computing, Mathematics and Digital Technology (SCMDT), Manchester Metropolitan University, M1 5GD Manchester, UK (R.D.); Center for Social Analysis and Mathematics (CAMS, UMR8557), CNRS, EHESS, 75013 Paris, France (A.S.). * These authors contributed equally to this work. ** These authors jointly supervised this work. Correspondence and requests for materials should be addressed to R.D. (email: mailto:[email protected]

Web End [email protected] ) or to K.M. (mailto:[email protected]

Web End [email protected]) or to N.P. (mailto:[email protected]

Web End [email protected])

NATURE COMMUNICATIONS | 7:8674 | DOI: 10.1038/ncomms9674 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms9674

Cells being the necessary level of integration of biological processes1, multicellular organization is best described by the cell lineage tree deployed in space and time. Thus, the

quantitative investigation of cell behaviour based on lineage branches annotated with relevant measurements at the individual cell level is the indispensable basis for reconstructing the multilevel dynamics of developing organisms2. Accurate and precise data about cell positions, trajectories, divisions, nucleus and cell shapes can be derived from the automated processing of 3Dtime images. Contributions in the eld point to the

necessary co-optimization of 4D multimodal imaging techniques and algorithmic image processing workows35. Ideally, going from the microscopy data to the interactive visualization of the cell lineage tree and segmented shapes should be automated, easily manageable and fast enough to allow a quantitative comparison of individuals6.

In recent years, decisive breakthroughs were made in the microscopy imaging of living systems, thanks to progress in uorescent protein engineering79 and microscopy imaging techniques, including multiphoton laser scanning microscopy (MLSM) and selective plane illumination microscopy (SPIM)10. Concomitantly, image processing methods for cell segmentation, cell tracking and the analysis of new types of quantitative data have diversied and improved35,1113. The huge data ow produced by 3Dtime imaging of live specimens has also greatly

beneted from faster computer hardware and computing grid architectures able to cope with high-dimensional data sets14. Finally, computer-aided data analysis and visualization software have completed the toolbox of quantitative developmental biology15.

Successful applications are still rare, however, and producing an accurate cell lineage tree for any developing organism remains a difcult challenge. In 2006, the automated reconstruction of the nematode cell lineage from confocal images established the rst standard16, although it did not yield reliable results beyond the 194-cell stage. Later, reconstructions were attempted on more complex organisms, such as the zebrash embryo imaged by digital scanned laser light-sheet uorescence microscopy (DSLM)17 or Drosophila imaged by MLSM during gastrulation18, but they did not provide long-term accurate single-cell tracking either. A concurrent work4 on semiautomated cell lineage reconstruction from harmonic generation imaging of non-labelled zebrash embryos provided six digital specimens with precise nucleus and membrane segmentation, yet was limited to the rst 10 divisions of the egg cell. Most recently, Amat et al.5 delivered a standalone software for the reconstruction of cell lineages from Drosophila embryos with uorescently stained nuclei. Their method, also tested on developing mouse and zebrash embryos, is well suited for the low background and high temporal resolution of SPIM data. Among all the state-of-the-art algorithmic image processing strategies, whether commercial or open-source software, the latter is the only one offering 3Dtime cell tracking with detection

of mitotic events to reconstruct the branching dynamics of cell lineage. Altogether, the growing number of solutions available today conrms that the automated reconstruction of cell trajectories and cell shapes, together with their interactive visualization, is at the cutting edge of developmental biology. Obviously, the performance achieved so far in terms of accuracy, scalability and ease of operation leaves plenty of room for improvement. There are still a great number of methods to explore, and more to invent in the elds of image processing and machine vision.

We deliver here an original image processing workow, BioEmergences, in the form of standalone software. Although optimized for MLSM data and fast cell movements in gastrulating

zebrash embryos, it generally performs well on 3Dtime

imaging data without heavy requirements in terms of spatial and temporal resolution, or signal-to-noise ratio. In addition to the reconstruction of the cell lineage branching process, the BioEmergences workow includes segmentation algorithms for cell nucleus and membrane shapes. These are based on the subjective surface method, which can complete cell contours from heterogeneous uorescent membrane staining19. The standalone version of the workow can be operated through a graphical user interface, and its output data are connected to Mov-IT, a custom-made interactive visualization software. Alternatively, our web service offers users customized assistance and fast processing on computer clusters or on the European Grid Infrastructure (EGI), together with the possibility to explore a large parameter space for the optimization of results (see Methods to request access).

We demonstrate the reconstruction and analysis of six digital embryos from three different species. All the data obtained, raw and reconstructed, is made available to the community. The BioEmergences workow is compared with eight other software tools from four different providers on the basis of gold standard data sets obtained by manual validation and correction of cell lineages. It scores best in all three tested categories: nucleus centre detection, linkage and mitosis detection. Thus, the combined BioEmergences/Mov-IT platform can contribute to the denition of standard procedures for the reconstruction of lineage trees from 4D in vivo data. The validation, annotation and analysis tools provided here support detailed, large-scale cell clonal analysis and characterization of cell behaviour along the lineage tree. This leads the way to the creation of benchmarks for a new type of interdisciplinary and quantitative integrative biology.

ResultsOverview of the BioEmergences workow. The phenomenological reconstruction of embryonic cell lineage starts from multimodal 4D data, typically comprising at least one channel for the uorescent signals emitted by cell nuclei, and possibly another one for cell membranes. Although nuclear staining is essential for cell tracking, membrane staining is necessary to assess cell morphology and neighbourhood topology. We deliver the rst public version of the standalone BioEmergences workow, with graphical user interface, able to launch a succession of algorithmic steps on two-channel raw data (Fig. 1 and Methods). Image ltering, nucleus centre detection, membrane shape segmentation and cell tracking methods were all designed and tuned to deal with the inherent noise and incompleteness of 4D imaging data generated by MLSM (Fig. 2). The BioEmergences standalone pipeline produces reconstructed data that can be directly displayed by the interactive visualization software, Mov-IT (Methods). By superimposing reconstructed data on raw data, Mov-IT adds visual inferences to create an augmented phenomenology. This allows the user to control data quality, measure the error rate, easily correct cell detection and tracking errors, and investigate the clonal history of cells and their behaviour.

The automated tracking of cells from 4D image data sets across whole living embryos involves a difcult trade-off between interdependent variables: signal-to-noise ratio, spatial and temporal resolution, imaging depth and cell survival. Microscopy techniques for in vivo and in toto imaging of developing organisms are evolving rapidly. In particular, the combination of two-photon excitation uorescence for deep-tissue imaging over extended periods of time20,21 with parallelized microscopy based on SPIM/DSLM is a promising approach22, although it is not yet widely available for routine time-lapse imaging in developmental biology.

2 NATURE COMMUNICATIONS | 7:8674 | DOI: 10.1038/ncomms9674 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms9674 ARTICLE

We explored here the coupling of two femtosecond lasers with two upright laser-scanning microscopes (Supplementary Fig. 1). Our methods are illustrated on 6-h spatiotemporal sequences

covering the cleavage and gastrulation periods in Danio rerio (Dr1 and Dr2) and Phallusia mammillata (Pm1 and Pm2), and the cleavage stages in Paracentrotus lividus (Pl1 and Pl2; Fig. 2 and Supplementary Movies 1, 2 and 3). The ascidian and echinoderm embryos were imaged in toto, allowing the reconstruction of the complete cell lineage across entire specimens. The subvolume imaged in the developing zebrash contained up to 9,500 cells, adding the difculties of high nucleus density and noise increase along the imaging depth. The user can choose between a standalone version, easily deployed on a personal computer but limited in terms of scalability, and a web service with user support, made to rapidly process large data sets through parallel implementation on our computer clusters and the EGI (Supplementary Note 1 on computational speed). The latter is especially powerful for optimizing the choice of algorithm parameters under expert supervision. The comparison of the BioEmergences workow with other available strategies demonstrates the performance and usefulness of our tools.

Image processing steps. The BioEmergences automated image processing workow (Fig. 1) is delivered as a standalone code that can be executed from the command line and, in part, from a graphical user interface. It starts by ltering the images and detecting the cell centres from local maxima, then performs a tracking of the cells trajectories, with optional segmentation of the shapes of nuclei and membranes, producing in the end a digital specimen in output (Methods, Supplementary Movies 4, 5 and 6).

Filtering and nucleus centre detection algorithms based on multiscale image analysis23 and partial differential equations (PDEs) are particularly useful for data sets that contain a high

Nuclear staining membrane staining

3D+time imaging

Nucleus channel membrane channel (optional)

Embryo preparation & mounting

Image acquisition

Automated image processingIn vivoimaging

Filtering Geodesic mean

curvature flow (GMCF)

Flux-based level set (FBLS)

Subjective surface(SubSurf) Tracking

Difference of gaussians (DoG)

Simulated annealing (SimAnn)

Mov-IT

Center detection

Segmentation

Interactive visualization

Validation, exploration, correction, annotation, analysis

Figure 1 | The BioEmergences reconstruction workow. From top to bottom: successive steps starting from embryo preparation and leading to the reconstructed data, all readily available for display and analysis by the interactive visualization tool Mov-IT. Each processing step is described in greater detail in the Methods. We propose two alternative nuclear centre detection methods. Either output can be used for shape segmentation and/or cell tracking.

Sphere

a b

4.2 h.p.f. 68 m 4.2 h.p.f. 10.9 h.p.f. 76 m

c

d

10.9 h.p.f.

1 s

NK

PCP

YSL

Sphere 1 s

e 5 h.p.f.

f

5 h.p.f.

Gastrula

g

8.1 h.p.f.

h

8.1 h.p.f.

Tailbud

Gastrula

Tailbud

i 5 h.p.f. 8.9 h.p.f.

j

5 h.p.f.

57-cell

k

l

8.9 h.p.f.

Blastula

57-cell

Blastula

Figure 2 | Reconstructing early embryogenesis from time-lapse optical sectioning. (ad) Danio rerio data set Dr1, animal pole (AP) view; (a,b) sphere stage, volume cut at 68 mm of the AP to visualize deep cells; (c,d) 1-somite stage (1 s), volume cut at 76 mm of the AP; NK, neural keel; PCP, prechordal plate; YSL, yolk syncytial layer. (eh) Phallusia mammillata data set Pm1, vegetal view, circum notochord side up; (e,f) gastrula stage; (g,h) tailbud stage.

(il) Paracentrotus lividus data set Pl1, lateral view, AP up; (i,j) 57-cell stage; (k,l) blastula stage. (a,c,e,g,i,k) 3D raw data visualization (with Amira software), hours post fertilization (h.p.f.) indicated top right, developmental stage indicated bottom left, scale bars, 50 mm. (a,c,i,k) Nuclear staining from

H2B-mCherry mRNA injection (in orange), membrane staining from eGFP-HRAS mRNA injection (in green). (e,g) Nuclear staining from H2B-eGFP mRNA injection (in blue-green). (b,d,f,h,j,l) Reconstructed embryo visualized with the Mov-IT tool, corresponding to (a,c,e,g,i,k), respectively. Each cell is represented by a dot with a vector showing its path over the next time steps (15 for the zebrash data, 9 for the ascidian and sea urchin data). Colour code indicates cell displacement orientation.

NATURE COMMUNICATIONS | 7:8674 | DOI: 10.1038/ncomms9674 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms9674

density of nuclei, as observed in the zebrash by the end of gastrulation and beyond. Starting from raw images (Fig. 3a), an initial image ltering step removes noise through geodesic mean curvature ow (GMCF)24 relying on a nonlinear geometrical PDE that can simultaneously smooth and sharpen the image25. Next, the position of each nucleus is found by a ux-based level set (FBLS) centre detection method (Fig. 3b)26. It identies objects with humps in the image intensity function, and places nuclear centres at local maxima. This is done through a nonlinear advection-diffusion equation that moves each level set of image intensity by a constant normal velocity and curvature.

Alternatively, we also provide a module performing well with a low density of nuclei and highly contrasted images. Based on a difference of Gaussians (DoG) convolution lter, it is able to simultaneously smooth the image and keep the most salient features. In the web service implementation of the workow, parameters of the Gaussians are interactively selected upon visual inspection of the resulting centres.

In an optional step, both nucleus and membrane geometries are automatically found by shape segmentation (Fig. 3ce) using

the subjective surface (SubSurf) technique initialized with the previously detected cell centres2729. The numerical discretization of GMCF, FBLS and SubSurf is based on the co-volume method30 and its parallel implementation19. The SubSurf method was also used to segment the global volume of imaged embryonic tissue, useful to estimate the average cell density and its evolution through time (Fig. 3fh). Finally, our original cell tracking algorithm inputs the list of approximate centres of cell nuclei and outputs their lineage in space and time following a three-step strategy. The rst step initializes the lineage links by a nearest-neighbour heuristic method. The second step uses simulated annealing31, a variant of the Metropolis algorithm, to progressively enforce a set of constraints reecting a priori biological knowledge. This is achieved by repeated random trials of link modication, and validation of possible changes according to a cost function, based on a weighted sum of contributing terms and a temperature parameter. The last step uses simulated annealing to link childless and motherless cells to the tree in an acceptable way, leaving open the possibility to make a posteriori changes in the detected nuclei.

a

10.1 h.p.f. z67 m

b c

10.1 h.p.f. z67 m

Raw

d

e f

4.2 h.p.f.

g

10.1 h.p.f. z67 m

h

10,000

0.80.60.40.2

8,000

4,000 2,000

Total cell number

6,000

Cell number per (10 m)3

5 h.p.f.

6

7

8

9

10

11

Figure 3 | Chain of PDE-based algorithms for the 3D segmentation of embryonic cells. Results from the zebrash data set Dr1. Scale bars in mm (100 in a,b,f,g, 10 in c,d,e). (a) Raw data image; single z section at 67 mm depth (orientation as in Fig. 2d), magnication in inset. (b) Nucleus centre detection by the FBLS method, after ltering by GMCF; approximate centres (cyan cubes) superimposed on raw data (grey levels, same orthoslice as a). Magnication in insets: left inset at depth 67 mm; right inset located one section deeper, showing that several centres not displayed at 67 mm were in fact correctly detected and visible below the chosen plane (white arrowheads); remaining centres can be found on other planes using the Mov-IT interactive visualization tool.

(c) Nucleus segmentation by SubSurf method; left panels: an interphase nucleus; right panels: a metaphase nucleus; top panels: initial segmentation (pink contour) superimposed on raw data (grey levels); bottom panels: nal segmentation (orange contours) superimposed on the same raw data. (d) Nucleus segmentation by SubSurf method; top panel: three nuclear contours superimposed on two raw-data orthoslices; bottom panel: 3D rendering of the segmented nuclei on a single orthoslice. (e) Membrane segmentation by SubSurf method; top panel: one cell membrane contour superimposed on two-raw data orthoslices; bottom panel: 3D rendering of the segmented membrane shape on a single orthoslice. (f) 3D rendering of segmented cell shapes.(g) Embryo shape segmentation; local cell density increasing from blue to red (same orthoslice as a). (h) Total cell number (blue curve) and cell density (brown curve) as a function of time; arrow indicates the end of gastrulation correlating with a plateau in cell density.

4 NATURE COMMUNICATIONS | 7:8674 | DOI: 10.1038/ncomms9674 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms9674 ARTICLE

Visualization and validation. After processing, each step of the reconstruction can be interactively visualized and analysed with the Mov-IT software. Developed to address the needs of biologists investigating the potential of in silico embryology based on 3Dtime microscopy imaging, Mov-IT allows easy validation

and correction of nucleus detection and cell lineage by superimposing reconstructed data on raw data. To this aim, it provides a complete set of tools including 3D volume rendering, 2D orthoslice views, cell lineage and segmentation display.

The Mov-IT software was used to validate the reconstruction of our six specimens cell lineage trees. The sea urchin (Pl1,2) and ascidian (Pm1,2) data sets were extensively checked and curated, producing quasi error-free digital specimens (gold standards) over an average of 30,000 temporal links. The sh data sets (Dr1,2) were only partly checked because of their large number of cells. More than 80,000 tracking links and 1,200 cell divisions were validated and, when necessary, corrected in the data set Dr1 chosen to establish our standard validation protocol (Fig. 4ad, Methods, Supplementary Movies 7, 8 and 9).

The performance of our automated workow on data set Dr1 was evaluated on three processing outputs: nuclear centre detection, linkage (by tracking cells one time step in the past) and mitosis detection. For each output, three success or error percentages were calculated: a sensitivity (representing the rate of true positives (TPs)), a false detection rate (representing false positives (FPs)) and a false negative (FN) rate. Finally, the product of centre detection sensitivity and linkage sensitivity produced a global lineage score (Methods, Supplementary Table 1 and Supplementary Fig. 2). On average over four windows of 21 time steps distributed in the rst 360 time steps of data set Dr1 (1h30 out of 6h40 of imaging from prior gastrulation to the 1 somite stage), the average linkage sensitivity of BioEmergences was 97.89%, meaning that only 2.11% of the cells were missed or not correctly tracked between time steps. FP nuclei were a negligible source of tracking errors (0.21%). False trajectories, in which tracking jumped from one cell to another between two time steps were observed in a very small number of cases, too(1.03%). The performance on small organisms, sea urchin and ascidian, was close to the gold standards for the best part of the images. Cell detection and tracking were occasionally poorer,

depending on the image quality, essentially the signal-to-noise ratio, which degraded with imaging depth.

The relevance and usefulness of the BioEmergences reconstruction workow are demonstrated by comparing its outputs with those obtained on the same dataset Dr1 using state-of-the-art commercial and open-source software. Eight image processing tools from four different providers were deemed suitable to handle gigabytes of time-lapse data. Their performance was tested on the detection of nuclear centres, links and cell divisions in several time intervals (Fig. 5, Methods and Supplementary Table 1). Based on measures of TP rates, called sensitivity, and FP rates, our methods produced the best results in every category. In particular, BioEmergences had the lowest rate of false linkage at 1% (Fig. 5h, column e). It also obtained an average mitosis sensitivity of 67% against 13% for Amat et al. and 0% for the rest (Fig. 5h, column c), meaning that the other tools were actually not designed to join lineage branches through the detection of divisions. Our software and Amat et al.s software were the only ones allowing the linkage of one mother cell at time step t to two daughter cells at time step t 1. In essence, mitosis

detection is a special task that only two packages among the nine that we tested were able to solve. Typically, it is at the time of anaphase and telophase that an algorithm can detect two separate nuclei, and as soon as it does it must also link them back to their mother. If this is not accomplished, the lineage tree will lack a junction. Even when a mitotic event is missed, most methods are still able to correctly track the resulting two cells, but these cells will not be properly annotated as daughters (Supplementary Movie 10).

Overall, our workow reached an average lineage score of 96% (the product of centre sensitivity and linkage sensitivity), whereas Amat et al. had 83% and all other tools remained below 50% (Fig. 5h, column g). Behind these average values, there was also a noticeable degradation of performance at later developmental stages in all systems except BioEmergences (Fig. 5g).

Analysis of cell fate and proliferation in digital specimens. The cell lineage tree is the cornerstone of a detailed understanding of morphogenetic processes (Supplementary Note 2 on in silico fate

5.80 h.p.f.

5.73 h.p.f. 5.95 h.p.f.

5.88 h.p.f.

a b

c

5.88 h.p.f. 5.90 h.p.f. 5.94 h.p.f.

11.61 h.p.f. 11.61 h.p.f.

d

Figure 4 | Visualization and validation of the lineage tree reconstruction (zebrash dataset Dr1). Results from the zebrash data set Dr1. All screenshots taken from the Mov-IT visualization interface, then tagged. (a) Cell division illustrated by three snapshots; time in hours post fertilization (h.p.f.) indicated top right; cell centres (cyan cubes) and cell paths (cyan lines) superimposed on two raw-data orthoslices showing the membranes (grey levels). (b) Flat representation of the cell lineage tree for three cell clones over 17 consecutive time steps: each cell is represented by a series of cyan squares, linked according to the cells clonal history; the cell dividing in Fig. 4a is circled. (c,d) Nucleus centre detection in a subpopulation of cells chosen at11.61 h.p.f. Mov-IT visualization in checking mode adding short vertical white lines to the detected centres, and displaying correct nuclei in green and false positives in red; (c) all detected nuclei; (d) validated nuclei only.

NATURE COMMUNICATIONS | 7:8674 | DOI: 10.1038/ncomms9674 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms9674

a b c

Sensitivity of centre detection

Sensitivity of linkage

Sensitivity of mitosis detection

100

80 60 40 20

90

50 30 10

%

100

80 60 40 20

90

50 30 10

%

100

80 60 40 20

90

50 30 10

%

70

70

70

0 4.36 6.22 8.08 9.95

h.p.f.

0 4.36 6.22 8.08 9.95

h.p.f.

0 4.36 6.22 8.08 9.95

h.p.f.

d

Rate of false centre detection

%

e f

Rate of false linkage

%

100

80 60 40 20

70

90

50 30 10

100

80 60 40 20

70

90

50 30 10

100

80 60 40 20

70

90

50 30 10

0 4.36 h.p.f.

6.22 8.08 9.95

0 4.36 6.22 8.08 9.95

Rate of false mitosis detection

%

h.p.f.

0 4.36 6.22 8.08 9.95

h.p.f.

g h

Lineage score

%

Bio emergences workflow

Imaris autoregressive motion expert

Imaris brownian motion

Imaris autoregressive motion

Imaris connected component

Icy spot tracking

Volocity shortest path

Volocity trajectory variation

Amat et al.5

100

80 60 40 20

70

90

50 30 10

100

80 60 40 20

70

90

50 30 10

0 4.36 6.22 8.08 9.95

Average rates

%

h.p.f.

0 a b c d e f g

Figure 5 | Comparative performance of nine software tools on zebrash data set Dr1. The BioEmergences workow presented here (red) was evaluated alongside Imaris (four shades of blue), Icy (yellow), Volocity (two shades of green) and Amat et al.5 (pink). Measurements were made by comparing the outcome of three reconstruction methods: nucleus detection, cell tracking (linkage) and mitosis detection, to a set of validated events (positions, trajectories, divisions) registered in a gold standard. This was done inside four time intervals centred in 4.36, 6.22, 8.08 and 9.95 h.p.f. (a) Sensitivity of nuclear centre detection, representing the rate of true positive (TP) centres. (b) Sensitivity of linkage, representing the rate of TP links (restricted to the subset of TP centers that possessed a validated link). (c) Sensitivity of mitosis detection, representing the rate of TP divisions (restricted to the subset of TP centres). (df) Rates of false positive (FP) centres, links and divisions (only two software tools applicable to the latter). (g) Global lineage score, equal to the linkage sensitivity times the centre detection sensitivity: all methods except BioEmergences deteriorate noticeably at later developmental stages. (h) Average rates calculated over the four time intervals: each column displays the mean heights of one of the previous seven charts. BioEmergences obtained the best results in every category: highest values in ac,g (success rates), lowest values in df (failure rates). Formulas can be seen in Methods, detailed values in Supplementary Table 1.

mapping). Digital specimens, such as the examples presented here, constitute a unique source of in-depth knowledge into embryonic patterning and the individuation of morphogenetic elds (Fig. 6). In the zebrash embryo, cells from the epiblast, hypoblast or epithelial enveloping layer were distinguished according to their position and behaviour32, and selected in the digital specimen with the Mov-IT tool. The selected cell populations were tracked, either backward or forward, revealing their relative movements with resolution at the cellular level during gastrulation stages (Fig. 6a,b and Supplementary Movie11). The more stereotyped cell lineage of the tunicate Phallusia mammillata allowed us to transpose the state-of-the-art fate map established at the 110-cell stage33. Cells were marked in the digital specimen according to their presumptive fate (Fig. 6c,d and Supplementary Movie 12). The propagation of the 110-cell stage

fate map along the cell lineage revealed the patterning of all the morphogenetic elds with a temporal resolution in the minute range. In the sea urchin embryo at cleavage stages, the cell membrane channel was used to mark cells according to their volume. This led to the identication of the three cell populations: micromeres, macromeres and mesomeres, organized along the vegetal-animal axis, with further segregation of macromeres into Vg1 and Vg2 subtypes, and micromeres into small and large subtypes34 (Fig. 6e,f and Supplementary Movie 13). The propagation of colours along the cell lineage showed that, despite limited cell dispersion, the cellular organization in the sea urchin embryo varied from one embryo to another.

The Mov-IT interface was also designed for fast import of processed data les in order to perform a systematic analysis of cell and tissue properties (Mov-IT tutorial). This is illustrated by

6 NATURE COMMUNICATIONS | 7:8674 | DOI: 10.1038/ncomms9674 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms9674 ARTICLE

a b

4.2 h.p.f. 10.9 h.p.f.

Sphere 1 s

32-cell Blastula

c d

5 h.p.f. 8.6 h.p.f.

e 4.3 h.p.f. 9.5 h.p.f.

f

Gastrula TB

Figure 6 | Automated fate map propagation. Cells were manually selected according to their identity or fate, using the interactive visualization tool Mov-IT. Developmental stages indicated bottom left, developmental times in h.p.f., top right. All scale bars, 50 mm. (a,b) Three cell populations in

Danio rerio specimen Dr1, animal pole view; (a) sphere stage, ventral (anterior) up, enveloping layer cells in cyan, epiblast cells in blue;(b) 1-somite stage, same colour code plus hypoblast cells in yellow.(c,d) Fate map in Phallusia mammillata specimen Pm1, vegetal view, circum notochord up, colour code as in ref. 33; (c) gastrula stage; (d) tailbud (TB) stage with automated propagation along the cell lineage of the fate map shown in c. (e,f) Cell populations in Paracentrotus lividus specimen Pl1, lateral view, animal pole up; small micromeres in dark purple, large micromeres in light purple; mesomeres in cyan. (e) 32-cell stage, macromeres in red; (f) blastula stage, Veg1 cell population in red,Veg2 cell population in yellow.

the analysis of the evolution of cell proliferation and cell density in our sea urchin, ascidian and zebrash specimens (Supplementary Fig. 3). Further statistical analysis and measurements are expected to contribute in a feedback loop to theoretical models and numerical simulations (Supplementary Fig. 4).

DiscussionThe BioEmergences reconstruction pipeline, accessible through a simple graphical user interface, was designed to run as quickly and efciently as possible from the acquisition of microscopy images to the display of cell lineage and cell segmentation aligned with raw data. Although the standalone version was crafted for convenient processing of small data sets on a laptop computer, there is no size limit on the data that can be uploaded to the BioEmergences web service through iRODS, an open-source data management software and the OpenMOLE engine14 for data processing on the EGI.

The concept of augmented phenomenology, coined to describe the overlay of raw and reconstructed data required for validation, correction and analysis, is fully exploited using the custom-made Mov-IT software. This tool serves in particular to

demonstrate that the reconstructed cell lineages meet the best quality of precision and accuracy. So far, cell lineage data can only be validated by eye inspection or by comparison with available gold standards, which are established manually and cross-checked by at least two experts. Corrected data sets are considered error-free as long as experts do not dispute this conclusion. Although gold standards themselves depend on what the eye can achieve, it is generally accepted that automated image processing software is still on average less effective than human vision. Our cross-software comparison method based on a set of validated events (cell positions, temporal links, divisions) is a step in the direction of standardized validation and comparison protocols. All the results led to the conclusion that the BioEmergences workow achieved the best performance on standard data produced by MLSM imaging with a temporal resolution chosen to explore either a small specimen (sea urchin or tunicate embryo) or up to one-third of the zebrash gastrula. The software provided by Amat et al. was the next best, but its performance dropped faster than BioEmergences with the increase in cell density and tissue thickness.

The cell tracking pipeline delivered here should be useful for a large variety of model organisms that have adequate optical properties to allow the acquisition of 3Dtime data sets,

including a channel for stained nuclei. The requirements in terms of spatial and temporal resolution and signal-to-noise ratio are easily fullled with commercial microscopy setups, either MLSM or SPIM. However, the accuracy of cell lineage reconstruction depends to a great extent on the quality of the data. Biologists who try for the rst time the automated processing of their time-lapse images might have to adjust their staining and imaging scheme to increase the quantity of information collected, especially across tissue depth. Image processing outcome also depends on the selection of algorithm parameters. This might require a few trials in the standalone version, if the default parameters do not already lead to satisfactory results. The web service version, on the other hand, offers the advantage of fast processing on a grid infrastructure, allowing the concurrent execution of hundreds of parameter combinations that can be easily explored by the user through Mov-IT, as done for the DoG centre detection method.

The detection of cell nuclei is a critical step, as it is used not only for cell tracking but also for nucleus and membrane shape segmentation, which must be performed on validated nuclei to avoid major errors but is computationally expensive. Although the standalone software performs well on a few selected cells, full-scale shape segmentation can be achieved in a reasonable amount of time only on a computing cluster such as the one accessible through the BioEmergences web service. Segmentation output can also be veried with Mov-IT, but its quantitative validation remains an issue35.

The public availability of the BioEmergences platform and its application to the embryogenesis of model organisms is intended to open the path to in silico embryology based on digital specimens. We illustrate this ambition with fate-map studies. The possibility of constructing complete fate maps in digital specimens, as we achieved for the tunicate embryos, is a revolution in the eld. Finally, although delivering here a standalone software, we also expect to initiate through the web service option a synergistic effort of the scientic community towards further validation, correction and annotation of digital specimens.

Methods

Embryo staining and mounting. Wild-type Danio rerio (zebrash) embryos were injected at the one-cell stage with 200 pg H2B-mCherry and 200 pg enhanced green uorescent protein (eGFP)-HRAS mRNA prepared from PCS2 constructs36,37.

NATURE COMMUNICATIONS | 7:8674 | DOI: 10.1038/ncomms9674 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 7

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms9674

Although mCherry was not as bright as eGFP and bleached more through imaging, this colour combination was the best compromise for a proper staining of cell membranes and further segmentation. Injected embryos were raised at 28.5 C for the next 3 h. Embryos were mounted in a 3-cm Petri dish with a glass coverslip bottom, sealing a hole of 0.5 mm at the dish centre, where a Teon tore (ALPHAnov) with a hole of 780 mm received the dechorionated embryo. The embryo was maintained and properly oriented by inltrating around it 0.5% low-melting-point agarose (Sigma) in embryo medium38. Temperature control in the room resulted in B26 C under the objective, slightly slowing down development with respect to the standard 28.5 C developmental table32 (Supplementary Table 2). After the imaging procedure, the embryo morphology was checked under the dissecting binocular and the animal was raised for at least 24 h to assess morphological defects. Embryo survival depended on total imaging duration, average laser power and image acquisition frequency (time step Dt). The zebrash data sets Dr1 and Dr2 were imaged through a standardized procedure with Dto150 s allowing up to 120 sections per time point with an average laser power of 80 mW delivered to the sample for more than 10 h without detectable photodamage. Lowering the laser power to less than 60 mW and increasing Dt up to 3.5 min allowed imaging embryos for more than 20 h, then raising them to adulthood. These conditions were used to obtain up to 320 sections at a 400 Hz line scan rate (bidirectional scanning) or 200 Hz to improve signal-to-noise ratio.

Oocytes from Phallusia mammillata (ascidian) were dechorionated and injected as described39 with 1 mg ml 1 H2B-eGFP mRNA prior fertilization. Membrane staining was obtained through continuous bathing in articial sea-water containing

FM 464 (Life Technologies) at a concentration of 1.6 mg ml 1. Embryos were deposited in a hole made in 1% agarose in ltered sea water at the centre of a 3-cm

Petri dish.

Oocytes from Paracentrotus lividus (sea urchin) were prepared and injected as described40 with 150 mg ml 1 H2B-mCherry and 150 mg ml 1 eGFP-HRAS synthetic mRNA. Embryos were either maintained between slide and coverslip covered with protamin41 or embedded in 0.25% low-melting-point agarose sea water at the centre of a 3-cm Petri dish.

Image acquisition. Imaging was performed with Leica DM5000 and DM6000 upright microscopes SP5 MLSM, equipped with an Olympus 20/0.95NA W dipping lens objective or a Leica 20/1NA W dipping lens objective. Axial resolution at the sample surface (1.5 mm) was estimated by recording 3D images of 0.1 or1 mm uorescent polystyrene beads (Invitrogen) embedded in 1% agarose.

For the zebrash specimen Dr1, the eld size was 700 700 mm2 in x, y and

142 mm in z, with a voxel size of 1.37 1.37 1.37 mm3 and a time step of 67 s.

For the zebrash Dr2, the eld size was 775 775 mm2 in x, y and 164 mm in z, a

cubic voxel of edge 1.51 mm and a time step 153 s. For the ascidian Pm1, these

dimensions were 384 384 mm2 in x, y, 165 mm in z, voxel 0.75 0.75 1.5 mm3

and time step 180 s. For the ascidian Pm2: 353 353 mm2 in x, y, 188 mm in z, voxel

0.69 0.69 1.39 mm3 and time step 180 s. For the sea urchin Pl1: 280 280 mm2

in x, y, 86 mm in z, voxel 0.55 0.55 1.09 mm3 and time step 207 s. For the sea

urchin Pl2: 266 266 mm2 in x, y, 82 mm in z, voxel 0.52 0.52 1.05 mm3 and

time step 180 s. For two-colour acquisition, simultaneous two-photon excitation42 at two different wavelengths (1,030 and 980 nm) was performed with pulsed laser beams (T-pulse 20, Amplitude Systmes and Ti-Sapphire femtosecond oscillator Mai Tai HP, Newport Spectra physics, respectively). Details of the optical bench are provided in Supplementary Fig. 1. Raw-data movies were made with the Amira software (Mercury Computer Systems).

Image processing algorithms. Explanation of the parameters and their useful range, along with the specic values used to process data sets Dr1-2, Pm1-2 and Pl1-2, are all provided in Supplementary Table 3. Filtering by GMCF relied on the following PDE25,43,44:

@tu ru

j jr g rGs u

j j

r

u

ru j j

ru

j j2

q r r

u

@

e2

0

1

q

A; 3

where g gjrGs ufNj for nucleus shapes, g gjrGs ufMj for membrane

shapes, and e, wc, wd are parameters. A detailed description of the role of these parameters, with explicit and semi-implicit numerical schemes, along with a discussion of the computational results is given in refs 28,29. The same equation (3) is also used for the overall embryo shape segmentation29.

Our cell tracking algorithm uses simulated annealing (SimAnn). It takes in input the list of approximate centres of cell nuclei detected at each time step, and produces in output a lineage forest, a graph equal to the union of several disjoint trees, in which each cell present at the rst time step was the root of a lineage tree. A graph is composed of a set of edges, or links, each of them connecting a nucleus centre at time t to a nucleus centre at time t 1. The algorithm was implemented

by an automated three-stage process:

First, edges between centres at consecutive time steps were initialized using a nearest-neighbour heuristic method. This created a set of links that were not necessarily biologically plausible.

Second, simulated annealing31, a variant of the Metropolis algorithm45, was used to progressively enforce a set of predened constraints summarizing together a certain number of biological requirements on the lineage forest, most notably: each cell beyond the rst time step should come from a single cell at a previous time step; each cell should have a single mother (a corresponding centre at t 1);

no cell should have more than two daughters (corresponding centres at t 1); no

cell should disappear (there is no cell death at these developmental stages in the chosen speciesan assumption invalid in other cases, such as mammalian preimplantation embryos); divisions should not occur too frequently; cell displacements should be bounded and so on. At the start, a nite set of allowed link modications and a cost function F measuring the departure from the constraints

were dened. Then, the algorithm relied on random link visits, and acceptance or rejection of a potential link change based on its cost. More precisely, F was the

weighted sum of local contributions, each addressing one of the a priori biological requirements. The cost function weights were selected after some initial trials of the algorithm and visual inspection of the resulting lineages, then recorded in a conguration le and kept constant during computation. This second stage consisted of repeatedly selecting and tentatively modifying a link at random from the list of permitted moves, then evaluating the resulting change through the cost function. In case of cost decrease, the candidate change was systematically accepted, whereas in case of cost increase, it was accepted with a probability proportional to exp( DF/T), where T is a temperature parameter that was

progressively lowered (annealed) over time so that fewer and fewer breaches to the constraints were accepted. The temperature decrease schedule was linear with time. The cost function also included terms to penalize deformations, to favour symmetry in the behaviour of sister cells, to add noise-countering inertia, to bound speed, and to account for division times.

Third, the nal goal was to minimize false-positive and false-negative errors in the nuclear detection steps. This was achieved by looking at the whole biological coherence of the lineage tree and identifying lineage gaps or cells that lived only a few time steps. The algorithm could then delete the centres of short branches or introduce virtual centres to regain continuity. To nd a best solution, it used simulated annealing again. A parallel implementation was written, which partitioned space-time into different cylinders that were run on different processors, then merged the results.

Validation protocol. Obtaining a complete gold standard annotated reference, even for small animals such as the ascidian and sea urchin embryos, was possible only with considerable human effort. The Phallusia and Paracentrotus data sets, Pm1-2 and Pl1-2, were almost completely checked and curated. Cell tracking accuracy depended on the imaging depth, along which the signal-to-noise ratio degraded. Owing to the decrease in image quality with depth, there remained a few errors for which there was no solution, even through manual expertise. We curated around 25,000 temporal links in Pl1, 22,000 in Pl2, 30,000 in Pm1 and 40,000 in Pm2. All divisions were also fully annotated. The four digital embryos produced were validated by two independent experts.

ru

j j2

: 1

It was accompanied by the initial condition u 0; x

u0N x

, or u 0; x

u0M x

,

where u0N x

and u0M x

are the image intensities of the nuclei and the membranes,

respectively, depending on which channel was ltered. In the GMCF model, the mean curvature motion of the level sets of image intensity is inuenced by the edge indicator function g(s) 1/(1 Ks2), KZ0, applied to the image intensity gradient

pre-smoothed by convolution with a Gaussian kernel Gs of small variance s. Details of the numerical method for solving equation (1) and its (parallel)

computer implementation are given in refs 25,43.

Centre detection by FBLS dened the nuclei as the local maxima of a smoothed version of the original image. Our algorithm was based on the following PDE26:

@tu d ru

j j m ru

j jr r

u

ru

j j

; 2

where the initial condition was given by the intensity function of the ltered nucleus image ufN. Equation (2) represents the level-set formulation for the motion of isosurfaces of solution u by a normal velocity V d mk, where d and m are

constants and k is the mean curvature. Owing to the shrinking and smoothing of all level sets, the function ufN was simplied, and we observed a decrease in the

number of spatial positions of local maxima, which could be used as approximate nuclear centre positions. Details of the numerical solution to equation (2) can be found in refs 26,43.

Alternatively, smoothing and centre detection could be achieved by DoG,a convolution of the image with two Gaussians of different standard deviations, here 1.52.5 mm and 1216 mm respectively. Their difference was calculated and the gray values above a threshold between 1 and 10% were selected. This allowed to simultaneously smooth the image and keep the most signicant objects. We ran multiple simulations combining different possible values of standard deviations and thresholds. Optimal values were visually chosen with Mov-IT by interactively checking the detection results.

Nucleus and membrane segmentation extracted the shapes of cell nuclei and/or membranes by evolving an initial segmentation function based on the subjective surface (SubSurf) equation27:

@tu wcrg ru wdg

e2

8 NATURE COMMUNICATIONS | 7:8674 | DOI: 10.1038/ncomms9674 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms9674 ARTICLE

In the case of the zebrash, however, a complete gold standard was almost impossible due to a number of cell samples of the order of several million. Our gold standard for the Danio data set Dr1 was obtained by a large-scale curation strategy:(i) all detected nuclear centres were completely annotated inside 54 spatial boxes of average size 340 70 140 mm3 in x, y, z, over ten time steps starting at various

sample times t0A{1, 20, 100, 200, 300, 350}. By this method, 25,123 links were manually checked, comprising 3,262 links for t0 1 (from t 1 to t 10), 2,735 for

t0 20 (from t 20 to t 29), 3,005 for t0 100, 8,736 for t0 200, 5,474 for

t0 300 and 1,911 for t0 350. This strategy insured even sampling throughout the

image data both in time and space (Fig. 4). (ii) As part of the investigation of the zebrash fate map, several cell clones were manually validated and when required corrected along with the whole data set, corresponding to 50,503 additional curated temporal links and the corresponding curated nuclei. (iii) An entire layer of epithelial cells was completely curated at one specic time, t 4. (iv) In addition,

a large number of links were randomly checked, bringing the total to 80,428.(v) Finally, more than 600 mitoses were checked.

Accuracy estimation protocol. Taking the consistent and representative gold standard of data set Dr1 as a basis, we were able to automatically identify and count various error types in the reconstructed embryos (Supplementary Table 1). This produced TP, FP and FN numbers, which were used as a performance metric for BioEmergences and other software tools. The accuracy of a given processing workow was evaluated on three outputs: nuclear centre detection, linkage (tracking cells one time step in the past) and mitosis detection.

By denition, the sum of TP and FN centres was equal to the total number of centres in the gold standard. This relationship did not hold, however, for TP and FN links or divisions because counting those events was restricted to the subset of TP centres. For links, FN was the sum of the numbers of wrong links (WL) and missing links (ML). This is because WL, which connect a cell to a wrong target, contributed both to FP (by creating new links that do not exist) and to FN (by missing the correct links), whereas ML, which correspond to a cell without any link, contributed only to FN.

Then, for each output, three types of success and error percentages were measured: a sensitivity, equal to TP/(TP FN), a false detection rate, equal to

FP/(TP FP) and a FN rate, equal to FN/(TP FN). Finally, a global lineage

score was calculated as the product of centre detection sensitivity and linkage sensitivity. The rationale for this formula is that linkage alone could appear successful even if many centres were missing, therefore it should be weighted by the actual proportion of detected centres.

Software performance and comparison. We identied eight state-of-the-art tools most relevant for our benchmark comparison: Icy (Spot Tracking, http://icy.bioimageanalysis.org

Web End =http:// http://icy.bioimageanalysis.org

Web End =icy.bioimageanalysis.org ), Imaris (Autoregressive Motion Expert, Autoregressive Motion, Brownian Motion and Connected Component, http://www.bitplane.com/imaris/imaris

Web End =http://www.bitplane.com/ http://www.bitplane.com/imaris/imaris

Web End =imaris/imaris ), Volocity (Shortest Path and Trajectory Variation, http://cellularimaging.perkinelmer.com/downloads

Web End =http:// http://cellularimaging.perkinelmer.com/downloads

Web End =cellularimaging.perkinelmer.com/downloads ) and the last method publishedby Amat et al.5,46. For each software tool, if a command line mode was available, we ran several reconstructions using different sets of parameters and selected the most advantageous one. For example, in the case of Amat et al.,we found that an optimal conguration was backgroundThreshold 16 and

persistanceSegmentationTau 0. If the software provided only a

visualization interface, we produced the best reconstruction by visual inspection based on a few parameter variations.

We used data set Dr1 in input to all methods, preprocessed through our ltering algorithm (GMCF, 5 iterations). The nucleus detection and tracking outputs were compared with our gold standard data. To assess nucleus detection, we explored the neighbourhood of validated centres by looking for other detected centres at a distance of 0.20.6 times the average internuclear distance. The number of correct links was estimated by inspecting the subset of detected nuclei labelled TP that also possessed a link in the gold standard, and counting the correct links from t to t 1. Mitosis detection was also assessed within the set of TP nuclei.

All measurements were made over four separate windows of 21 time steps each: from t 0 to t 20 (corresponding to 4.360.18 h.p.f.), from t 100 to t 120

(6.22 h.p.f.), from t 200 to t 220 (8.08 h.p.f.) and from t 300 to t 320

(9.95 h.p.f.). The nal sensitivity, false detection and FN rates were averaged over these four intervals. Detailed scores are shown in Supplementary Table 1.

Data sets and software. The standalone BioEmergences workow, the visualization tool Mov-IT and the six in vivo 3Dtime image data sets are

provided online at http://bioemergences.iscpif.fr/bioemergences/openworkflow-datasets.php

Web End =http://bioemergences.iscpif.fr/bioemergences/openworkow http://bioemergences.iscpif.fr/bioemergences/openworkflow-datasets.php

Web End =datasets.php . For each specimen, we provide the 4D raw-data images and the corresponding reconstructed embryo. The parameters used to process these data sets are provided in Supplementary Table 3.

Our BioEmergences platform is also available as a web service at http://bioemergences.iscpif.fr/workflow/

Web End =http:// http://bioemergences.iscpif.fr/workflow/

Web End =bioemergences.iscpif.fr/workow/ , which offers users customized assistance in addition to fast processing. The online workow architecture relies on iRODS, an open-source data management software, and the OpenMOLE engine14,a middleware platform facilitating the experimental exploration of complex systems models on a computing cluster, which leverages the power of the EGI.

Interested users are invited to request access to these resources by sending an email to http://[email protected]

Web End [email protected] .

References

1. Jacob, F. & Spillmann, B. E. The Logic of Life: a History of Heredity (Princeton Univ., 1993).

2. Oates, A. C., Gornkiel, N., Gonzalez-Gaitan, M. & Heisenberg, C.-P. Quantitative approaches in developmental biology. Nat. Rev. Genet. 10, 517530 (2009).

3. Fernandez, R. et al. Imaging plant growth in 4D: robust tissue reconstruction and lineaging at cell resolution. Nat. Methods 7, 547553 (2010).

4. Olivier, N. et al. Cell lineage reconstruction of early zebrash embryos using label-free nonlinear microscopy. Science 329, 967971 (2010).

5. Amat, F. et al. Fast, accurate reconstruction of cell lineages from large-scale uorescence microscopy data. Nat. Methods 11, 951958 (2014).

6. Luengo-Oroz, M., Ledesma-Carbayo, M., Peyrieras, N. & Santos, A. Image analysis for understanding embryo development: a bridge from microscopy to biological insights. Curr. Opin. Genet. Dev. 21, 630637 (2011).

7. Mavrakis, M., Rikhy, R., Lilly, M. & Lippincott-Schwartz, J. Fluorescence imaging techniques for studying Drosophila embryo development. Curr. Protoc. Cell Biol. 39, 4.18.14.18.43 (2008).

8. Chudakov, D. M., Matz, M. V., Lukyanov, S. & Lukyanov, K. A. Fluorescent proteins and their applications in imaging living cells and tissues. Physiol. Rev. 90, 11031163 (2010).

9. Buckingham, M. E. & Meilhac, S. M. Tracing cells for tracking cell lineage and clonal behavior. Dev. Cell 21, 394409 (2011).

10. Fischer, R. S., Wu, Y., Kanchanawong, P., Shroff, H. & Waterman, C. M. Microscopy in 3D: a biologists toolbox. Trends Cell Biol. 21, 682691 (2011).

11. Meijering, E., Dzyubachyk, O., Smal, I. & van Cappellen, W. A. Tracking in cell and developmental biology. Semin. Cell Dev. Biol. 20, 894902 (2009).

12. Meijering, E., Dzyubachyk, O. & Smal, I. Methods for cell and particle tracking. Method. Enzymol. 504, 183200 (2012).

13. Khairy, K. & Keller, P. J. Reconstructing embryonic development. Genesis 49, 488513 (2011).

14. Reuillon, R., Leclaire, M. & Rey-Coyrehourcq, S. OpenMOLE, a workow engine specically tailored for the distributed exploration of simulation models. Future Gener. Comp. Sy. 29, 19811990 (2013).

15. Eliceiri, K. W. et al. Biological imaging software tools. Nat. Methods 9, 697710 (2012).

16. Bao, Z. et al. Automated cell lineage tracing in Caenorhabditis elegans. Proc. Natl Acad. Sci. USA 103, 27072712 (2006).

17. Keller, P. J., Schmidt, A. D., Wittbrodt, J. & Stelzer, E. H. Reconstruction of zebrash early embryonic development by scanned light sheet microscopy. Science 322, 10651069 (2008).

18. Supatto, W., McMahon, A., Fraser, S. E. & Stathopoulos, A. Quantitative imaging of collective cell migration during Drosophila gastrulation: multiphoton microscopy and computational analysis. Nat. Protoc. 4, 13971412 (2009).

19. Suri, J. S. & Farag, A. Deformable Models: Theory and Biomaterial Applications vol. 2 (Springer, 2007).

20. Squirrell, J. M., Wokosin, D. L., White, J. G. & Bavister, B. D. Long-term two-photon uorescence imaging of mammalian embryos without compromising viability. Nat. Biotechnol. 17, 763767 (1999).

21. Helmchen, F. & Denk, W. Deep tissue two-photon microscopy. Nat. Methods 2, 932940 (2005).

22. Planchon, T. A. et al. Rapid three-dimensional isotropic imaging of living cells using Bessel beam plane illumination. Nat. Methods 8, 417423 (2011).

23. Alvarez, L., Guichard, F., Lions, P.-L. & Morel, J.-M. Axioms and fundamental equations of image processing. Arch. Ration. Mech. An. 123, 199257 (1993).

24. Caselles, V., Kimmel, R. & Sapiro, G. Geodesic active contours. Int. J. Comput. Vision 22, 6179 (1997).

25. Kriva, Z. et al. 3D early embryogenesis image ltering by nonlinear partial differential equations. Med. Image Anal. 14, 510526 (2010).

26. Frolkovic, P., Mikula, K., Peyrieras, N. & Sarti, A. Counting number of cells and cell segmentation using advection-diffusion equations. Kybernetika 43, 817829 (2007).

27. Sarti, A., Malladi, R. & Sethian, J. A. Subjective surfaces: a method for completing missing boundaries. Proc. Natl Acad. Sci. USA 97, 62586263 (2000).

28. Zanella, C. et al. Cells segmentation from 3D confocal images of early zebrash embryogenesis. IEEE T. Image Process 19, 770781 (2010).

29. Mikula, K., Peyrieras, N., Remesikova, M. & Stasova, O. Segmentation of 3D cell membrane images by PDE methods and its applications. Comput. Biol. Med. 41, 326339 (2011).

30. Corsaro, S., Mikula, K., Sarti, A. & Sgallari, F. Semi-implicit covolume method in 3D image segmentation. SIAM J. Sci. Comput. 28, 22482265 (2006).

31. Kirkpatrick, S., Gelatt, C. D. & Vecchi, M. P. Optimization by simulated annealing. Science 220, 671680 (1983).

NATURE COMMUNICATIONS | 7:8674 | DOI: 10.1038/ncomms9674 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 9

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms9674

32. Kimmel, C. B., Ballard, W. W., Kimmel, S. R., Ullmann, B. & Schilling, T. F. Stages of embryonic development of the zebrash. Dev. Dynam 203, 253310 (1995).

33. Lemaire, P. Unfolding a chordate developmental program, one cell at a time: invariant cell lineages, short-range inductions and evolutionary plasticity in ascidians. Dev. Biol. 332, 4860 (2009).

34. Angerer, L. M. & Angerer, R. C. Animal-vegetal axis patterning mechanisms in the early sea urchin embryo. Dev. Biol. 218, 112 (2000).

35. Rizzi, B. & Sarti, A. Region-based PDEs for cells counting and segmentation in 3Dtime images of vertebrate early embryogenesis. Int. J. Biomed. Imaging

2009, 968986 (2009).36. Moriyoshi, K., Richards, L. J., Akazawa, C., OLeary, D. D. & Nakanishi, S. Labeling neural cells using adenoviral gene transfer of membrane-targeted GFP. Neuron 16, 255260 (1996).

37. Shaner, N. C., Steinbach, P. A. & Tsien, R. Y. A guide to choosing uorescent proteins. Nat. Methods 2, 905909 (2005).

38. Westereld, M. The Zebrash Book: a Guide for the Laboratory Use of Zebrash (Danio rerio) (Univ. of Oregon Press, 2000).

39. Sardet, C. et al. Embryological methods in ascidians: the Villefranche-sur-Mer protocols. Methods Mol. Biol. 365400 (2011).

40. McMahon, A. P. et al. Introduction of cloned DNA into sea urchin egg cytoplasm: replication and persistence during embryogenesis. Dev. Biol. 108, 420430 (1985).

41. Summers, R. G., Piston, D. W., Harris, K. M. & Morrill, J. B. The orientation of rst cleavage in the sea urchin embryo, Lytechinus variegatus, does not specify the axes of bilateral symmetry. Dev. Biol. 175, 177183 (1996).

42. So, P. T., Dong, C. Y., Masters, B. R. & Berland, K. M. Two-photon excitation uorescence microscopy. Annu. Rev. Biomed. Eng. 2, 399429 (2000).

43. Bourgine, P. et al. 4D embryogenesis image analysis using PDE methods of image processing. Kybernetika 46, 226259 (2010).

44. Chen, Y., Vemuri, B. C. & Wang, L. Image denoising and segmentation via nonlinear diffusion. Comput. Math. Appl. 39, 131149 (2000).

45. Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H. & Teller, E. Equation of state calculations by fast computing machines. J. Chem. Phys. 21, 10871092 (1953).

46. Amat, F., Myers, E. W. & Keller, P. J. Fast and robust optical ow for time-lapse microscopy using super-voxels. Bioinformatics 29, 373380 (2013).

Acknowledgements

We thank Patrick Parra and Jean-Yves Tiercelin for custom-made high-precision

mechanics; Amparo Ruiz Vera and Odile Lecquyer for the management of the

Embryomics and BioEmergences EC projects; Audrey Colin, Valrie Lavalle and

Simon Agostini for sh care; Roger Tsien and Richard Harland for DNA constructs

(mCherry and eGFP-HRAS, respectively); Pascal Abbas for technical assistance; the N&D

and ISC-PIF staff for administrative support; Dominique De Waleffe from DENALI SA for

his contribution to the BioEmergences EC projects workow interface. This work was

supported by grants EC NEST Adventure 012916, EC NEST Measuring the Impossible

28892, ZF-Health EC project HEALTH-F4-2010-242048, ACI STIC sant, ANR BioSys

Morphoscale, ANR-11-BSV5-021-03 REGENR, ANR-10-INSB-04 and ANR-11-EQPX-

0029 for N.P. and P.B., ARC 3628 for N.P. and APVV-0184-10 for K.M.

Author contributions

E.F. coordinated the workow development, the manuscript preparation and performed

the software comparison; T.S. developed the visualization interface Mov-IT and the

BioEmergences web service; B.R. developed the Gaussian-based nucleus centre detection

method; C.M. developed the cell tracking method and contributed to the workow

development; O.S. developed ltering and nucleus centre detection methods; D.F. con

tributed to the specication of the Mov-IT functionalities and organization of the

BioEmergences database; R.S. developed the embryo shape segmentation method and

contributed to the workow development; M.H. contributed to the workow develop

ment, in both its standalone and web service versions; R.C. developed algorithm paral

lelization tools; G.R. coordinated the specication of the Mov-IT functionalities and

organization of the BioEmergences database; B.L. worked on the reconstruction valida

tion and standardization of error assessment; L.D. conducted imaging experiments, data

validation and correction; I.C. and S.D. conducted imaging experiments, data validation

and correction; J.K. provided engineering support to the computation; P.A., B.M. and

A.B. conducted data validation and correction; J.-Y.N. and P.C. provided support for data

management, transfer and storage, and cluster management; P.V. contributed to dening

the BioEmergences workow concept; M.F. conducted data validation and correction;

G.L. led the molecular biology of embryo staining; Y.K. co-supervised the design

and implementation of the tracking method, and co-wrote the manuscript;

P.S. custom-designed the MLSM coupling and designed the imaging strategies;

M.R. co-coordinated the project, co-developed nucleus and membrane segmentation

methods and co-wrote the manuscript; R.D. co-coordinated the project, co-designed cell

tracking methods and co-wrote the manuscript; A.S. co-coordinated the project, and led

several algorithm designs; K.M. co-coordinated the project, led several algorithm designs,

and co-wrote the manuscript; N.P. co-founded and co-coordinated the project, led the

biological relevance and co-wrote the manuscript; P.B. co-founded and co-coordinated

the project, designed the formal and epistemological framework and co-wrote the

manuscript.

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/naturecommunications

Web End =http://www.nature.com/

http://www.nature.com/naturecommunications

Web End =naturecommunications

Competing nancial interests: The authors declare no competing nancial

interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions/

Web End =http://npg.nature.com/

http://npg.nature.com/reprintsandpermissions/

Web End =reprintsandpermissions/

How to cite this article: Faure, E. et al. A workow to process 3Dtime microscopy

images of developing organisms and reconstruct their cell lineage. Nat. Commun. 7:8674

doi: 10.1038/ncomms9674 (2016).

This work is licensed under a Creative Commons Attribution 4.0

International License. The images or other third party material in this

article are included in the articles Creative Commons license, unless indicated otherwise

in the credit line; if the material is not included under the Creative Commons license,

users will need to obtain permission from the license holder to reproduce the material.

To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

Web End =http://creativecommons.org/licenses/by/4.0/

10 NATURE COMMUNICATIONS | 7:8674 | DOI: 10.1038/ncomms9674 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications