AMEBaS: Automatic Midline Extraction and Background Subtraction of Ratiometric Fluorescence Time-Lapses of Polarized Single Cells

Cell polarity is a fundamental biological process in which the concerted action of a collection of spatially concentrated molecules and structures culminates in the establishment of specialized morphological subcellular domains¹. Cell division, growth and migration rely on such polarity sites, while its loss has been associated with cancer in epithelial tissue-related disorders².

Apically growing cells are a dramatic example of polarity, where the polarity site at the tip typically reorients to extracellular cues³. These include developing neurites, fungal hyphae, root hairs, and pollen tubes, where multiple cellular processes show pronounced differences from the tip of the cell toward the shank. In pollen tubes, in particular, actin polymerization, vesicle trafficking, and ionic concentrations are markedly polarized, showing tip-focused gradients⁴. Pollen tubes are the male gametophytes of flowering plants and are responsible for delivering the sperm cells to the ovule by growing exclusively at the apex of the cell at one of the fastest growth rates known for a single cell. The tip-focused gradients of ions such as calcium⁵ (Ca²⁺) and protons⁶ (H⁺) play a major role in sustaining pollen tube growth, which is essential to accomplish its main biological function that culminates in a double fertilization^5,6. Thus, quantitative methods to analyze the spatiotemporal dynamics along the midline of apically growing cells are essential to investigate the cellular and molecular mechanisms underlying polarized growth^7,8,9. Researchers often use kymographs, i.e., a matrix that represents the pixel intensities of the cell's midline (e.g., columns) through time (e.g., rows), which allows visualizing cell growth and migration in the diagonal (Figure 1). Despite their usefulness, kymographs are frequently extracted by manually tracing the midline, being prone to biases and human errors while also being rather laborious. This calls for an automated method of midline extraction that is the first feature of the pipeline introduced herein named AMEBaS: Automatic Midline Extraction and Background Subtraction of ratiometric fluorescence time lapses of polarized single cells.

In terms of experimental procedures, quantitative imaging of ions/molecules/species of interest in single cells can be achieved with genetically encoded fluorescent probes¹⁰. Among the ever-expanding choices, ratiometric probes are one of the most accurate since they emit different fluorescence wavelengths when bound/unbound to the molecules of interest¹¹. This allows for correction of the spatial heterogeneity in the intracellular concentration of the probe by using the ratio of two channels with their channel specific background subtracted. However, estimating the background threshold for each channel and time point can be a complex task since it often varies in space due to effects like shading, where the corners of the image have luminosity variation relative to the center, and in time due to fading of the fluorophore (photobleaching)¹². Although there are multiple possible methods, this manuscript proposes determining the background intensity automatically using the segmentation threshold obtained with the Isodata algorithm¹³, which is then smoothed across frames through polynomial regression as a standard. Spatial components stemming from fluorescence heterogeneity unrelated to the target cell removed in¹², however, were ignored by this method. Automatic thresholding can be performed by several methods, but the Isodata algorithm produced the best results empirically. Thus, automatic background value subtraction and ratiometric calculation are the second main feature of AMEBaS (Figure 1), which, taken together, receives as input a stack of dual-channel fluorescence microscopy images, estimates the cell's midline and the channel-specific background, and outputs kymographs of both channels and their ratio (main output #1) after background subtraction, smoothing, and outlier removal, together with a stack of ratiometric images (main output #2).

AMEBaS was tested with fluorescence time lapses of growing Arabidopsis pollen tubes obtained under a microscope, either with Ca²⁺ (CaMeleon)⁸ or pH (pHluorin)⁶ ratiometric sensors expressed under the pollen-specific LAT52 promoter. Images from each channel were taken every 4 s coupled to an inverted microscope, a front-illuminated camera (2560 pixels × 2160 pixels, pixel size 6.45 μm), a fluorescence illuminator, and a water immersion objective lens 63x, 1.2NA. Filters settings used for CaMeleon were: excitation 426-450 nm (CFP) and 505-515 nm (YFP), emission 458-487 nm (CFP) and 520-550 nm (YFP), while for pHluorin, excitation 318-390 nm (DAPI) and 428-475 nm (FITC), emission 435-448 nm (DAPI) and 523-536 nm (FITC). A complete data set was added for testing at Zenodo (DOI: 10.5281/zenodo.7975350)¹⁴.

In addition, the pipeline was tested with root hair data, where imaging was performed with a light sheet microscope (SPIM) as previously described^15,16 with Arabidopsis root hairs expressing the genetically encoded Ca²⁺ reporter NES-YC3.6 under the control of the UBQ10 promoter¹⁷. The home-made LabView software that controlled the camera acquisition, sample translation and shutter of the light sheet microscope permitted the observation of the two cpVenus and CFP channels, but also the visualization of their ratio in real time. Every ratio image of the time-lapse represented a maximum intensity projection (MIP) between the cpVenus and CFP fluorescent channels images obtained from 15 slices of the sample spaced 3 μm apart. The time-lapse cpVenus/CFP ratio of MIPs was saved and directly used for the AMEBaS analysis.

Although this pipeline can work with multiple types of growing and migrating cells, it was specifically designed to analyze growing cells that grow exclusively at the tip, such as pollen tubes, root hairs, and fungal hyphae, where there is a correspondence of the non-growing cytoplasmic regions between frames. When such a correspondence is not present, the user should choose the complete_skeletonization option in step 1.3.1.1 (see the Discussion section for more details).

Figure 1: An overview of the pipeline workflow. The AMEBaS pipeline analyses and processes microscopic time lapses in three main steps: Single-Cell Segmentation, Midline Tracing, and Kymograph Generation. Please click here to view a larger version of this figure.

1. Interactive notebook protocol

The Jupyter notebook can be used directly on the web using Google Colab at https://colab.research.google.com/github/badain/amebas/blob/main/AMEBAS_Colab.ipynb, where the instructions below were based. Alternatively, the Jupyter notebook is available at https://github.com/badain/amebas, where it can be downloaded and configured to run locally in Jupyter (Anaconda can provide an easy and cross-platform installation process). A complete test data can be found at Zenodo (https://doi.org/10.5281/zenodo.7975350), containing single and dual channel data of Arabidopsis pollen tubes expressing either pH or Ca​²⁺ reporters¹⁴. The pipeline has been divided into parts, where every step can be executed by clicking on the play button after setting the user-specific options. The necessary files for this study are available in AMEBaS- main zip folder (Supplementary Coding File 1).

1. Open Jupyter notebook and read the time-lapse files.
   1. Navigate to the homepage of the interactive notebook in Google Colab referred above or download and open the AMEBaS_Local.ipynb notebook from GitHub.
   2. Prepare the directory setup for the input and output data:
      1. If using the local version, place the fluorescence time-lapse as a TIFF file or DV file inside a folder named data that must be in the root folder of the program. A folder named out must be created to receive the generated data. Then run the Setup code block.
      2. If using the notebook on Google Colab, run the Setup code block to automatically generate the data and out folders.
   3. Run the File Input code block to read the time lapse data by clicking the play button. If using the Google Colab version of the notebook, click on the Choose File button to directly upload the time lapse file to the data folder.
      ​NOTE: The number of channels will be automatically detected based on the dimension of the image.
   4. Choose if additional outputs of each step will be generated by setting the 'verbose' parameter to True or False.
2. Detect the main cell and segment from the background (Figure 2).
   1. Run the Single Cell Segmentation code block to automatically separate the cell of interest from the background by clicking the play button.
      NOTE: Median and Gaussian filters will be applied as a preprocessing step to remove unwanted noise before segmenting the foreground from the background by Isodata thresholding and isolating the region of the largest area to remove undesired artifacts.
      1. Adjust the sigma value used by the Gaussian in the 'sigma' variable to fine-tune the smoothness of the segmentation mask. The default value is 2.0.
      2. Set the variable estimate to False to store the threshold estimated from Isodata directly or to True to smooth it across neighboring frames using local polynomial regression (LOESS). Fine-tune its function by changing the n_points variable. The default value is 40.
3. Trace midline along the cell extension (Figure 3).
   1. Run the Cell Midline Tracing code block by clicking the play button to automatically skeletonize the cell using Lee's method¹⁸ and extend the tip of the last skeleton through linear extrapolation.
      1. Choose to trace the midline only on the last frame or once per frame by adjusting the complete_skeletonization argument.
         NOTE: When all frames are skeletonized, extrapolation is skipped.
      2. Set the fraction of points in the skeleton to be interpolated during extrapolation by adjusting the interpolation_fraction variable. The default value is 0.25.
      3. Choose the length of the midline extrapolation by changing the variable extrapolation_length. The default is -1, which extends the skeleton up to the nearest edge.
4. Generate kymographs for each channel ( Figure 4).
   1. Run the first Data Visualization code block by clicking the play button to automatically generate kymographs for both channels.
      1. Choose the size of the Gaussian kernel used for smoothing by adjusting the variable kymograph_kernel.
         NOTE: It corresponds to the size of the neighborhood (in pixels) over which the pixel intensities are averaged. The default is 3 pixels x 3 pixels.
      2. Nonextended skeletons generate capped Kymographs that must use a custom colormap to properly display their intensities. Choose the fraction percentage of the intensities that will be assigned to the background color, black, adjusting the shift_fraction variable. The default value is 0.7.
5. Calculate the ratio between channels (Figure 5).
   1. Run the second Data Visualization code block by clicking the play button to automatically generate a ratiometric kymograph and a ratiometric timelapse (Figure 6).
      NOTE: This step is only available when using dual-channel time lapses. The background intensity threshold stored in step 1.2.1.2 is subtracted from each channel.
      1. Adjust the switch_ratio variable to switch the order of channels used as the numerator and denominator during ratio calculations. The default is False.
      2. Choose if the ratio timelapse must be further smoothed with a median filter pass by adjusting the smooth_ratio variable. The default is False.
      3. Choose whether to remove outliers produced by the low signal of the denominator channel by manipulating the reject_outliers variable. Defaults to True and defines outliers as values 1.5 times the interquartile range above the third quartile (where 75% of the values lie).
      4. Choose whether the background in the ratiometric output must be exported by adjusting the variable background_ratio. The default is False, which replaces it with zeros.

Figure 2: Single cell segmentation step. Image processing techniques such as filtering, thresholding, and area labeling are used to isolate the signal of interest (step 1.2). This particular data had the following values for Lowest Intensity: 2556, Median: 3441, and Highest Intensity: 32125. Please click here to view a larger version of this figure.

Figure 3: Midline tracing overview- The midline of the single cell is obtained by calculating its skeleton (white). The tip (magenta) is extrapolated linearly from the last points at the end of the skeleton (step 1.3). In this composition both the midline and its tip are superposed over the original cell. Please click here to view a larger version of this figure.

Figure 4: Timelapse's kymographs- Comparison of the kymographs for each channel generated with 'complete_skeletonization' turned off (step 1.4). The vertical axis describes the progression of time, and the horizontal axis plots the average intensity of the extrapolated midline path followed by a single cell. For this particular data, the colormap represents the following values for Channel 1 Lowest Intensity: 2886, Median: 3167, Highest Intensity: 21021. Channel 2 Lowest Intensity: 3030, Median: 3400, Highest Intensity: 29688. Please click here to view a larger version of this figure.

Figure 5: Background threshold smoothing- The background segmentation threshold is estimated via the isodata algorithm and then smoothed via local polynomial regression (step 1.5). Please click here to view a larger version of this figure.

Figure 6: Ratiometric results- (A) Comparison between the last frame of the ratiometric timelapse and the segmented original first channel. (B) Kymograph generated from the ratiometric timelapse (step 1.5). Please click here to view a larger version of this figure.

2. Batch mode protocol

1. Download and place the file pipeline.py from the AMEBaS GitHub repository in the same directory as the data.
2. Type in the file location on the command line after the program file.
3. Include --v as a positional argument to show the internal steps of the pipeline, if desired.
4. Include --s to choose in the sigma value used in the Gaussian Filter preprocessing step in preparation for cell segmentation. The default is 2.
5. Include --a to trace the midline for each frame of the time lapse. By default, the pipeline uses only the last frame.
6. Include --f to choose the fraction [0,1] of the skeleton to be used in the interpolation. The default value is 0.25.
7. Include -e to choose the length in pixels of the extrapolated skeleton. The default is -1, which extends the skeleton up to the nearest edge.
8. Include --sf to choose the fraction of the color range that will be shifted to the background in non-extrapolated kymographs. The default value is 0.7.
9. Include --k to determine the size of the kernel used in the kymograph Gaussian filtering. The default is 3.
10. Include --eb to estimate global background threshold intensity via LOESS polynomial regression of the frame-specific background threshold intensities.
    1. Customize the number of points used in LOESS smoothing of the background threshold values by modifying the parameter --n. The default value is 40.
11. Switch the channels used as the numerator and denominator during ratio calculations, including -r or --switch_ratio, if the timelapse has two channels. By default, the second channel is the numerator and the first is the denominator.
12. Choose if the ratio timelapse must be further smoothed with a Median Filter pass with the --sm argument. The default is False.
13. Include -o to reject pixels with abnormal intensities during ratiometric time-lapse generation.
14. Choose if the background in the ratiometric output must be exported using the argument --b. The default is False, which replaces it with zeros.
15. Press enter to run. The output will be generated in the same directory as the program file.

The AMEBaS pipeline automates the extraction of midline dynamics of polarized single cells from fluorescence microscopy image stacks, making it less time consuming and less prone to human errors. The method quantifies these time lapses by generating kymographs and ratiometric image stacks (Figure 1) in growing single cells. It can be adjusted to work on migrating single cells, but further experiments are necessary. AMEBaS is implemented in Python as an interactive Jupyter Notebook (described in the section Interactive notebook protocol), allowing easier usage without requiring programming experience, and as a command line tool (in the section Batch mode protocol), where multiple stacks can be analyzed with the same parameter set. Although either one or two fluorescence channels can be used, ratiometric probes with two emission channels should yield more reliable results since the fluorescence emission of the unbound state of the probe can alleviate the spatial heterogeneity caused by the uneven distribution of the protein in the cytoplasm.

The pipeline first yields a binary mask of the largest cell on the strongest channel, exported as a tif file named Filename_binary_mask.tiff (Figure 2). Threshold estimates obtained with isodata for each channel are optionally smoothed with loess and saved in the table Filename_background_treshold.csv (Figure 5). The cell's midline, extracted from the binary mask, is exported as a tif file named Filename_skeletonized.tiff (Figure 3). Kymographs of each channel are produced from the midline, named Filename_kymograph_c_*.csv, where * corresponds to the channel number (Figure 4). Finally, a ratiometric kymograph is saved as Filename_kymograph_ratio.csv, while the full ratiometric stack is named Filename_ratiometric.tiff (Figure 6). Plots corresponding to Figure 2, Figure 3, Figure 4, Figure 5, and Figure 6 are optionally saved as PNG files when 'verbose == True' or '--v' is specified by the user in the first code chunk (step 1.1).

These results can be coupled to other image analysis pipelines to further investigate the spatiotemporal dynamics of the cell, such as CHUKNORRIS⁸, which takes the kymographs of each channel as an input and performs growth rate analysis with subpixel resolution, together with various time-series analysis methods.

AMEBaS was tested on pollen tube data sets with genetically encoded ratiometric fluorescent reporters for intracellular Ca²⁺ (CaMeleon; Figure 7A,B) and H⁺ (pHluorin; Figure 7C,D) concentrations acquired on an optical fluorescence microscope, as well as maximum intensity projection of cpVenus/CFP ratios of growing root hairs expressing the CaMeleon NES-YC3.6 (Figure 7E,F) acquired on a light-sheet microscope as previously described^15,16. The pipeline worked successfully despite differences in the growth direction, imaging techniques, fluorescent reporters, and cell types. Segmentation, midline tracing, and kymograph extraction are shown for these datasets (Figure 7), demonstrating the potential of applying AMEBaS to a wide range of experimental setups.

Figure 7: Representative results for different fluorescence image datasets of tip-growing cells. (A,B) Pollen tube expressing the Ca²⁺ reporter CaMeleon; (C,D)Pollen tubes expressing the pH indicator pHluorin; (E,F) Root hairs expressing the Ca²⁺ reporter NES-YC3.6. Please click here to view a larger version of this figure.

Supplementary Coding File 1: AMEBaS-main.zip. Please click here to download this File.

The novel method presented here is a potent tool to streamline and automate the analysis of fluorescence microscopy image stacks of polarized cells. Current methods described in the literature, such as ImageJ Kymograph plugins, require manual tracing of the midline of the polarized cell of interest, a task that is not only time consuming but also prone to human errors. Since the definition of the midline in this pipeline is supported by a numerical method^18,19 that performs skeletonization, subjective evaluation is removed, introducing a quantitative standard to the procedure. This is particularly useful for researchers with large quantities of data, with the possibility of customizing the pipeline to different demands, including extracting the midline for every frame. Furthermore, the background value subtraction is often arbitrary, with a single background value chosen by hand for all frames of a given channel. Here, isodata segmentation is used to objectively determine the background threshold for each frame, smoothing it (optionally) with a local polynomial regression to capture long-term changes in fluorescence caused by fading (photobleaching). While the possible artifacts and secondary elements are ignored in the background by selecting the larger object by pixel area, other methods such as FRET-IBRA¹², can be used to remove effects such as shading. Spatial artifacts (such as shading) can influence the shape of the cell segmented by thresholding, which may account for the bias of the gradient toward one side of the tube seen in the ratiometric movie featured on the AMEBaS logo (see GitHub page or Colab notebook).

Nonetheless, there are still a few situations where the presented pipeline may fail and should not be considered as a silver bullet solution. Special attention must be paid during image capture because large disturbances can impair the ability of the algorithm to segment the target cell. Preprocessing the image stack should be considered for optimal results, removing unwanted elements and overall faulty frames.

Parameters were chosen aiming to analyze single apically growing cells considering restricted data set of Arabidopsis pollen tubes assaying ionic concentrations with fluorescent reporters. Thus, other data may require sensible choices of parameters. Preprocessing filters smooth the data with the objective to produce a clean binary mask after segmentation (step 1.2), so the sigma value used for the Gaussian and median filters can be increased for noisier data or decreased if the results are overly smoothed. The location of the binary mask obtained in the strongest channel is assumed to be the same as the weakest channel, which can be a problem if the cell location is not the same in both channels. If this is the case, either different masks must be used for each channel or image registration must be performed, as done in FRET-IBRA¹².

Skeletonization is done exclusively in the last frame (step 1.3) by default, which assumes a tip-growing cell that maintains the location of its cytoplasm over time (other than at the apex). By extending the skeleton, it is possible to generate kymographs that allow analysis of the growth rate, even with subpixel resolution, using methods like CHUKNORRIS⁸. This extension is done by linear extrapolation considering the first 25% points in the growing tip of the skeleton (step 1.3) and can be tuned by adjusting the interpolation_fraction from 0%to 100% points of the entire skeleton. However, if the cell drifts in location between frames or is a migrating cell, growth rate analysis becomes more complex. In such a scenario, independent skeletons may be generated for every frame with the choice of parameters complete_skeletonization=TRUE, which will produce skeletons without extracellular extension. It would still be possible to analyze the growth rate, but the resolution would be limited by the binary mask produced with isodata thresholding. Furthermore, the resulting kymograph assumes that consecutive skeletons can be aligned by their default coordinates, which, if not true, renders it unsuitable for analyzing intracellular dynamics.

When generating kymographs (step 1.4), AMEBaS uses averaging with a Gaussian kernel instead of using the traditional pixel margin around the midline. The default value of a 3x3 is the smallest possible size, which may be increased if the data is too noisy and/or if the cell is large. However, this step may cause over-smoothing, in which case the filter should be turned off altogether by setting kymograph_kernel = 0. Finally, the background estimated for each frame of each channel can be smoothed when the user expects fading (photobleaching) or using the raw estimates (step 1.5). Smoothing of the background values with LOESS polynomial regression (step 1.5) can be tuned by setting the number of points used in the window n_points to a minimum of 3 and a maximum of the number of frames in the stack (capped automatically). A simpler function describing the background change over time can be attained with a larger window, yielding a coarser fit.

Since this method consolidates multiple manual steps typically conducted by researchers, the AMEBaS pipeline is a more efficient, unbiased, and precise approach tool to analyze the spatiotemporal behavior of single polarized cell time-lapses. In the future, this method has the potential to be extended to support the analysis of migrating single cells. Moreover, further analysis is necessary to assess the performance of this method across a wider range of cell types.