Quantitation of cellular dynamics in growing Arabidopsis roots with light sheet microscopy.

Plant imaging and growing chamber For each experiment a plantlet was grown in a closed and sterile micro-chamber ( Fig. 1A ) designed to constrain the root growth to an optically transparent vertical narrow path, while at the same time providing illumination and gas exchange for the leaves. The main body of the chamber is made up of a glass fluorimetry cuvette (inner size 10×10×40 mm, Starna) with a glass capillary tube (o.d. 1.5 mm, Fisher) placed vertically against one corner of the cuvette to create a channel (diameter approximately 300 µm) in which the root grows. Two teflon collars with a thickness of 80 µm are added at the ends of the tube to distance it from the walls of the cuvette and create enough space for medium flow in the channel, while still preventing root escape. The tube is held in position by filling the cuvette outside the channel with 3 mm glass beads (Fisher), and the whole space is filled with liquid medium (see Perfusion, below). The cuvette is closed with a custom-made transparent top cap and the junction sealed with several rounds of Parafilm. The cap is obtained from a disposable plastic spectrophotometer cuvette (inner size 10×10×40 mm, Fisher) whose square base has been removed and whose three long sides have been cut to create windows. The base is then closed with a glass square coverslip (for illumination) glued with silicone and the three side windows with adhesive breathable sealing tape (for gas exchange, Nunc). The sample is illuminated from above by an array of 12 bright white LEDs (Luxeon Star/O) powered by a current source controlled by the computer. The LED current is gated by a signal indicating camera exposure so that the LEDs are shut off during exposure. The optics and micro-chamber are housed in a thermally insulated light-tight box with a hinged door allowing sample loading and a thermoelectric assembly (AC-046, TE Technologies, Inc.) maintains constant temperature within the box.

Perfusion

To sustain root growth for days, continuous perfusion of liquid medium is maintained through the micro-chamber. One peristaltic pump (Dynamax RP-1, Rainin, run at 30 rpm with Pharmed BPT tubing 0.51 mm i.d., Cole Parmer) flows liquid medium from a reservoir to the bottom of the micro-chamber, and a second pump (Dynamax RP-1, Rainin, run at 30 rpm with Pharmed BPT tubing 0.89 mm i.d., Cole Parmer) flows medium out of the top of the glass cuvette (right below plantlet's hypocotyl) to a waste bottle. The outflow rate is kept higher than the inflow rate (different i.d.), to prevent flooding. The liquid medium is composed by 1/4× Murashige and Skoog (MS) basal medium enriched with Gamborg's vitamins (Sigma), 3.0% sucrose, 0.05% MES buffer (Sigma) and adjusted to pH 5.7 (KOH). After autoclave sterilization, antibiotic antimycotic solution (A5955, 1/200×, Sigma) and Amphotericin B (A9528, 1.0 mg/l, Sigma) are added the medium. The medium reservoir is maintained saturated with oxygen with filtered air from a standard aquarium air pump.

Plant material and growth conditions

Arabidopsis seeds are sterilized as describe previously [6] and germinated on 0.8% agar plates containing 1× MS, 0.5% sucrose, 0.05% MES buffer, pH 5.7, under constant illumination and 23°C. Between 3 and 5 days post-germination a single plantlet is transferred to the micro-chamber, first making sure that the root is completely aligned along the corner of the cuvette, and then reconstructing the micro-chamber by adding the capillary tube, filling the cuvette with the glass beads, adding some liquid medium and sealing the top cap with Parafilm. The perfusion is started, and the automated tracking and scanning sequence initiated. Throughout the experiment, the micro-chamber is illuminated with constant light at 20,000 lux and the temperature maintained constant at 23°C. Optical apparatus Our imaging system is based on the light sheet microscopy (LSM) concept [9] . Fluorescence excitation is provided by either a 473 nm, 80 mW diode pumped solid state laser (Lasever) or a multi-line 60 mW Argon ion laser (Spectra Physics). The laser beam passes an excitation filter (Chroma) and optional neutral density filter (Thorlabs) and is modulated by a Pockels cell (350–50, Conoptics). Two flat mirrors direct the beam into a 10× Galilean beam expander (Thorlabs) illuminating a circular 9 mm iris (Thorlabs). The circular beam is focused in one dimension by a plano-convex cylindrical lens (f = 30 mm, Thorlabs), creating a “sheet” of 4.2±0.2 µm across 180 µm along the excitation axis. Emission is imaged with a 20× super long working distance objective (NA = 0.35, f = 19.9 mm, Nikon) coupled by a lens tube to an emission filter (Chroma) and cooled CCD camera (1392×1024, 4.65 µm square pixels, QICam, QImaging). Acquisitions are typically performed with 2×2 pixel binning. The x axis is defined as coincident with the optical axis of the cylindrical lens, the y axis is parallel to gravity and z is the optical axis of the objective ( Fig. 1B ). The sample is maintained in the coincident foci of the cylindrical lens and the objective by 3 translation stages coupled to the sample holder ( x,y and z axes, Fig. 1B ) and 2 translation stages coupled to the CCD camera and objective ( x and z axes). All translation stages are positioned by micrometers driven by servo motors (25 mm travel, 40 nm resolution).

Resolution

The axial resolution of the optical system is determined by the sheet thickness near its focus and the depth of field of the microscope objective. Since the former is inversely proportional to the numerical aperture of the cylindrical lens (0.15) and the latter is inversely proportional to the square of the objective's NA (0.35 2 = 0.1225), the sheet thickness dominates in determining the overall axial resolution. The mean of the sheet full width at half maximum (FWHM) over the region of interest is 4.2±0.2 µm, as measured by scanning fluorescent beads (0.3 µm diameter, Duke Sci.) embedded in Phytagel (0.5%). This sheet thickness is maintained within 5% over 180 µm along the excitation beam axis, which is sufficient to section the entire diameter of the root. The subsequent deconvolution step (see below) improves the axial resolution, as shown by the intensity profiles of neighbor nuclei ( Fig. S1 ).

Image acquisition

At the beginning of each experiment, the root tip is positioned in the focus of the excitation beam, whose profile was previously measured by scanning fluorescent beads as described above. Imaging typically occurred every 10 minutes. At the beginning of each scan, the focal plane of the objective is determined by exciting fluorescence from a central region of the root while imaging at several positions along the optical axis of the objective. At each position, the spatial gradient is estimated at each pixel and summed over the image, and the position corresponding to the maximum sum is chosen as the best focus for that excitation plane. During the scan, the sample is repeatedly moved through the laser sheet by stepping its stage in the z direction and acquiring the emission signal from each optical section (typical exposure time was 25–50 ms). At each step of the scan, the stage holding the camera and objective is repositioned in z to maintain the excitation plane in the objective focus, accounting for the defocus caused by variable optical path length through air and medium. The sectioning interval is chosen to satisfy the Nyquist criterion for the measured point spread function (typically 1.25–2.5 µm). The total scan depth (typically 200 µm) is larger than the root thickness (typically 100 µm) to ensure that the entire tip is sectioned. The camera acquisition triggers via TTL a pulse from a custom circuit to modulate the Pockels cell. Typically multiple acquisitions at each z step are averaged to reduce amplifier noise from the CCD. The duration of the scan (typically 90 s) is primarily limited by the repositioning of the motorized stages. Photoenergy Due to the fundamental difference between the way the optical sectioning is performed in LSM and in confocal laser scanning microscopy (CLSM), a scan consisting of n slices (for typical Arabidopsis root, n≈60 with 2.5 µm steps) will produce a total incident photoenergy at least n times reduced in LSM compared to CLSM. Moreover, the calculated irradiance in the focal plane is 57 W/cm 2 with LSM compared to 950 kW/cm 2 calculated for typical confocal imaging of a root expressing the same marker. Such 16,000-fold reduction in irradiance and 60-fold reduction of incident photoenergy for the LSM compared to CLSM allow acquisition of live plant roots at high temporal resolution without significant photo-bleaching or photo-toxicity, as seen in Figure 1D of the main text.

Root tracking

Because the root grows continuously throughout the experiment (approximately 1–3 µm/min) the sample must be repositioned between scans to maintain the root tip in the field of view. This is accomplished via optical registration of the most recent scan's images to those from the previously optically registered scan. The algorithm used for registration is a version of the Lucas-Kanade algorithm [23] . The algorithm calculates by iterated gradient descent the 3D rigid motion which minimizes the squared difference between two images, assuming that the motion is small compared to the length scale over which the image intensity varies. After the optimal rigid motion has been calculated, the stages are translated to reposition a point of interest in the most recently scanned image to a fixed location in image space, chosen to maintain the entire root tip within the field of view. The stage position is recorded at each scan. The root typically rotates about all three axes as it grows. The calculated rigid motion includes rotations, while the repositioning of the sample stage only involves translations - the former are necessary to maintain registration of the root tip. Tracking was accurate for the largest observed translations, which were 15% in the linear dimension of the field of view. Defocus along the excitation beam caused by movement along the x axis of the sample was compensated taking into account the refractive indices of air and medium (water).

Acquisition software

The acquisition and tracking is facilitated by custom C++ software with a graphical interface controlling the stages, camera, LEDs, perfusion and all scan parameters, and displaying relevant images and information from the acquisition as it proceeds. The most computationally demanding task during an acquisition is the sample tracking, which takes about one minute on a desktop computer with a dual core 3.0 GHz processor (Core 2 Duo E8400, Intel) and 2GB RAM.

Deconvolution

The stack at each time point is deconvolved off-line via expectation maximization (EM) assuming Poisson noise and using the point spread function measured by imaging 0.3 µm diameter fluorescent beads (Duke Scientific Corp.) [24] . Typically 30 iterations of the EM algorithm are used. 3D fast Fourier transforms were performed using FFTW [25] . The deconvolution significantly improves the separability of nuclei along the optical axis ( Fig. S1 ), without introducing artifacts affecting the analysis.

Segmentation

The stack of deconvolved images at each time point is automatically segmented by a routine adapted from a previously described method [12] . The main difficulties for image segmentation in our data are heterogeneity of nuclear size and intensity, and variation of optical quality with position within the root. In addition, the spatial extent of each nucleus appears comparable to the internuclear distance, occasionally making the separation of neighboring nuclei even more challenging. Most of the spatial variation can be ascribed to the depth of tissue through which the light passes along the two optical axes. These depths are estimated from images by constructing the three dimensional convex hull of all voxels above an intensity threshold chosen to separate nuclear signal from background. This convex hull approximates the boundary of root tissue and for each imaged voxel the distances along the optical axes to the hull can be calculated. The logarithm of maximum pixel intensity was fit to a second order polynomial in the two optical depths by least squares and the fit was subsequently used to normalize intensities within the root. A segmentation routine was developed consisting of an increasing series of global segmentation thresholds and a set of criteria for classifying a contiguous blob as a nucleus. Once a blob is classified as a nucleus, its voxels are no longer considered as the segmentation threshold increases. The criteria for classification are ranges for number of voxels, mean intensity, lateral size, lateral eccentricity and axial size. To ensure that nuclei are separated along the optical axis, the intensity profile along the optical axis is computed for each blob and tested for the absence of multiple maxima as a further criterion. A blob must fall within the specified range for each parameter to be classified as a nucleus. The ranges for classification resulted from training with 850 manually identified nuclei, located throughout the root tip. The output of the segmentation routine is a time series of centroid location, maximum and mean intensities, voxel count, lateral size, lateral eccentricity and axial size for each detected nucleus. The segmentation validation was performed on a set of 1303 manually identified nuclei. The routine failed to identify 90 of the nuclei in the test set (false negative rate, 6.9%), while labeling as nuclei 95 blobs not present in the test set (false positive rate, 7.3%) ( Table S1 ). Registration of segmented images and frame of reference The global alignment of nuclear positions over time is important, as the nuclear displacement across time points is the essential factor used in tracking nuclei. The sample tracking is predictive, in that it compares the previous two stacks in order to calculate a translation for the next scan. An improved alignment with finer translations and rotations can be calculated after the image has been acquired. For simplicity, we calculate this alignment based on nuclear positions extracted by segmentation, minimizing by iterated gradient descent the mean squared nuclear displacements over adjacent time points with respect to rotations and translations. This frame of reference is used for nuclear tracking. Finally, a new frame of reference is introduced in which the origin of the primary root axis m effectively coincides with the root tip, defined at a constant distance from a single nucleus identifiable in all time points, and the cylindrical angle θ is measured from the objective axis ( Fig. 3A ). All the data presented are in such reference frame, which facilitates the analysis of performance as a function of penetration depth.

Nuclear tracking

To follow nuclei through time and identify cell divisions, the output of the segmentation routine is analyzed by a custom algorithm for nuclear tracking. Our routine exploits similarity in nuclear position within the root from one time point to the next, since other features such as intensity and nuclear size fluctuate more than position due to biological and optical effects. Nuclear identity assignments are represented as an array of trees, with each tree corresponding to a unique identity, each node of a tree representing a segmented blob associated with that tree's nuclear identity, and an edge between two nodes indicating that the corresponding blobs represent the same nucleus at adjacent times. In the ideal case, each level of the tree corresponds to one scan of the experiment, but in practice false negatives in the image analysis necessitate skipping isolated time points for some nuclei, so that a nucleus at time can be directly associated with a nearby nucleus not just at time but also , etc. The energy to be minimized consists of a sum over all edges of all trees of squared spatial displacement, a penalty for temporal gaps, a penalty for terminal nodes and a penalty for branching. Each identified nucleus is numbered and for nucleus the position , time and maximum intensity are considered. An indicator function is defined as equal to when nucleus is directly equated with nucleus and otherwise. The energy is expressed in , since the squared displacement is the most important term and is the only term with physical meaning, so this scale eases interpretation of the various terms. The energy can be written as (1) where is a coefficient controlling the contributions from non-physical quantities, and are the number of nuclei identified with nucleus at immediately previous and later times along its trajectory, respectively: and and are functions defined as follows for penalizing nuclei with numbers of parents and children not equal to one, respectively. Although nuclear fusions ( ) do not occur biologically, they are allowed in the analysis with a finite penalty to prevent false positives from disrupting the continuity of nearby trajectories so that these events can be considered manually: The data contained many events in which two nearby nuclei were properly resolved by the segmentation routine for some time, followed by a period in which they were not resolved. When the nuclei are again resolved, the event resembles a cell division. In order to distinguish these events from true cell divisions, several properties unique to true divisions were incorporated in the branching penalty. In a manually generated set of over 400 true divisions, the displacement of centroids from the parent cell to each daughter cell was on average 2.4-fold greater than the average non-dividing displacement (2.6±0.6 µm for dividing cells compared to 1.1±0.8 µm for non-dividing cells over 10 minute scan intervals). In addition the displacements from the parent were typically in opposite directions. The branching term g(0) was designed to apply negligible penalty to true division events and significant penalty to most events which otherwise resembled divisions, based on these empirical observations. For a nucleus with two daughter nuclei and , define the displacements and and the deviations of these displacements from the average measured displacement upon division , where , along with a threshold deviation: . The penalty for node having children is defined as where accounts for the displacement of the daughter nucleus from the parent: with and accounts for the angle between the two displacements: The fundamental term in the energy is the squared displacement, since nuclear position is the strongest indicator of identity in our data. The penalties f (0) and g (0) for zero children and parents determine the length scales over which identifications can be made, as described below. The remaining terms are chosen to cope with specific events - specifically, segmentation errors and cell divisions. The values for the penalties are empirically constrained. As an example, we depicted in Figure S4 a typical situation: two nuclei, one at time t and one at time t +1, are separated by a distance Δ r . Consider two configurations: (A), in which the nuclei are connected, i.e. assigned to the same trajectory, and (B), in which the nucleus at time t has no child and the nucleus at time t +1 has no parent, i.e. not assigned to the same trajectory. From (1), the energy for (A) is (Δ r ) 2 , while the energy for (B) is f (0)+ g (0). Configuration (A) should be chosen to describe a situation where the same nucleus is observed at both time points with a true displacement Δ r , while configuration (B) describes events in which the first nucleus disappears after time t (by leaving the field of view, for example) and the nucleus observed at time t +1 is its neighbor. Since the algorithm favors lower energy configurations, to achieve the correct assignment f (0)+ g (0) must be large compared to the square of the largest observed nuclear displacement over the time-step ( d max = 5 µm), so that the energy of (B) is always higher than the one of (A) in the case of true displacement of a single nucleus. At the same time, f (0)+ g (0) cannot be much larger than the square of the typically observed internuclear distance ( d nn = 10 µm), otherwise distinct neighbor nuclei will tend to be wrongly assigned to the same trajectory as in configuration (A). Other parameters are determined by similar empirical considerations. The tracking parameters have not been systematically optimized; nevertheless their performance allows for useful analysis, as validated below. The scheme for optimizing identity assignments was as follows for all experiments. The initial configuration is chosen by equating each nucleus at time with its nearest neighbor at time . The move at each step is chosen randomly but with a bias for high energy nodes as targets, so that the probability of choosing a given nucleus for a move is proportional to the energy it contributes to the total. Once node has been chosen, the possible moves include cutting a connection by flipping from 1 to 0 for some , adding a connection by flipping from 0 to 1, or changing a connection by flipping to 0 and to 1 for some and . Moves are accepted according to the Metropolis criterion [26] . We found that allowing two serial moves, where one of the nodes involved in the first move is chosen as the target for the second move, improved the rate of convergence significantly, while allowing three or more such serial moves showed no improvement. We found that the critical temperature for annealing the energy defined in equation (1) was approximately . This was typically chosen as the initial temperature, the system was equilibrated for several thousand moves, then exponentially cooled. Nuclear tracking was manually validated by choosing 100 trajectories at random and identifying errors by considering high energy nodes and comparing to images when necessary. Since some ambiguity remains in the manual tracking, we calculated global lower and upper bounds on the error rate of 2.4% and 3.9%, respectively, over 32,067 total assignments ( Table S1 ). Nuclear tracking was performed on a dual core 2.66 GHz processor (Core 2 Duo E7300, Intel). Identification of cell divisions (or cell-lineage tracking) Each branching event in the set of identified nuclear trajectories was labeled as a potential cell division. The majority of these branching events arose from failure of the segmentation to properly separate two nearby nuclei at one time point followed by proper separation at subsequent time points. Additional dynamic features of true cell divisions were considered after the nucleus tracking was performed. These features did not enter the definition of the energy for tracking due to the expense of calculation arising from their temporally extended nature. This filtering step was implemented by constructing a support vector machine [27] (SVM) for classification based on nuclear size and intensity at time points within a temporal neighborhood of the identified division and on the durations of the trajectories for the parent nucleus and both daughter nuclei. A complete set of 461 true cell divisions was manually built, and the SVM was trained on a subset of 424 true cell divisions and 3805 spurious branchings, yielding an out of sample accuracy of 92.3% in five fold cross validation. The SVM was then used to classify all branching events and compared to the complete manual set, resulting in 340 true positives (73.8% of true divisions) and 58 false positives (12.6% of true divisions).

Source code

Code for deconvolution, segmentation and nucleus tracking is available and maintained as a SourceForge project (qanneal-celltrk.sourceforge.net).

Plant material and growth conditions

Arabidopsis seeds are sterilized as describe previously [6] and germinated on 0.8% agar plates containing 1× MS, 0.5% sucrose, 0.05% MES buffer, pH 5.7, under constant illumination and 23°C. Between 3 and 5 days post-germination a single plantlet is transferred to the micro-chamber, first making sure that the root is completely aligned along the corner of the cuvette, and then reconstructing the micro-chamber by adding the capillary tube, filling the cuvette with the glass beads, adding some liquid medium and sealing the top cap with Parafilm. The perfusion is started, and the automated tracking and scanning sequence initiated. Throughout the experiment, the micro-chamber is illuminated with constant light at 20,000 lux and the temperature maintained constant at 23°C.

Supporting Information Figure S1 Effects of deconvolution in the axial ( z ) dimension. ( A ) Raw (left column) and deconvolved (right column) images of three close nuclei, at different z -positions; bar, 7.5 µm. ( B ) Maximum intensity collected in the region marked by the red ellipses in (A), both for the raw (black) and the deconvolved (blue) images. The increased resolution in z after deconvolution is shown by the improved separation between peaks. (TIF) Click here for additional data file. Figure S2 Contraction of the root radius. Number of nuclei with radial coordinate ρ greater than 55 µm, as a function of time. The red dotted line indicates the time point (21 hours) chosen to distinguish “early” and “late” temporal windows. (TIF) Click here for additional data file. Figure S3 False positives. Positions of the false positive cell divisions found in the same dataset containing the true positives shown in Fig. 5 of the main text, superimposed on the smoothed density estimate of all the nuclei segmented at all time points throughout the same 29 hour interval. (TIF) Click here for additional data file. Figure S4 Tracking nuclei: example. Diagram showing the correct assignments in two typical scenarios (A and B) for two consecutive time-points ( t and t +1). ( A ) The same nucleus a is observed at both times t and t +1, with a spatial displacement (blue arrow): the two nuclei should be assigned to the same trajectory (link between them, red arrow). ( B ) Nucleus a is only observed at time t and nucleus b is only observed at time t +1: the two nuclei should not be assigned to the same trajectory (no link between them, red truncated links). Δr , distance between nuclei observed in t and t +1; d max = 5 µm, maximum nucleus displacement observed between t and t+1 (time-step = 10 min); d nn = 10 µm, average distance between nearest neighbor nuclei. (TIF) Click here for additional data file. Table S1 Performance of segmentation and tracking routines. Segmentation: error rates in nuclei recognition in the regions I–V defined in figure 2C ; N, number of manually identified nuclei; FP, false positive rate; FN, false negative rate. Tracking: error rates in nucleus-trajectory assignment in the same regions I–V; N, number of manually identified trajectories starting in the given region; error range, lower and upper bounds for errors in assigning a nucleus to the correct trajectory. (PDF) Click here for additional data file. Video S1 Growing root. Time-lapse images of the same 29 hour interval used to extract the presented data. Maximum intensity projection along the optical axis of the deconvolved images. The temporal step is 10 minutes. Longitudinal tissue stretching and several cell divisions are detectable by eye. (MOV) Click here for additional data file.