Using a fluorescent marker for cytokinesis (Figure 3—figure supplement 6a), we observed that cellular growth has two phases (Figure 3—figure supplement 6b–c). During G1, the mother’s growth rate peaks; during S/G2/M, which we identify by the cells having buds, the bud dominates growth with its growth rate peaking approximately midway to cytokinesis (Ferrezuelo et al., 2012).

This tight coordination between bud growth rate and cytokinesis suggested that the peak in bud growth rate preceding cytokinesis may mark a regulatory transition. Comparing growth rates over S/G2/M for buds in different carbon sources, we found that the maximal growth rate occurs at similar times relative to cytokinesis despite substantial differences in the duration of the S/G2/M phases (Figure 4a).

Figure 4 with 1 supplement see all

Download asset Open asset

Daughters born in rich media are larger than those born in poor media, and some of this regulation occurs during S/G2/M (Johnston et al., 1977; Jorgensen et al., 2004; Leitao and Kellogg, 2017). Understanding the mechanism, however, is confounded by the longer S/G2/M phases in poorer media (Leitao and Kellogg, 2017; Figure 4a), which counterintuitively allow daughters that should be smaller longer to grow.

Given that the time between maximal growth and anaphase appears approximately constant in different carbon sources (Figure 4a), we hypothesised that the growth rate falls because the bud has reached a critical size. Compared to how their sizes vary immediately after cytokinesis, buds have similar sizes when their growth rates peak — in all carbon sources (Figure 4b): the longer S/G2/M phase in poorer media compensates the slower growth rates. During the subsequent constant time to cytokinesis, the faster growth in richer carbon sources would then generate larger daughters, and we observe that the bud’s average growth rate correlates positively with the volume of the daughter it becomes (Figure 4—figure supplement 1). Cells likely therefore implement some size regulation in S/G2/M as they approach cytokinesis.

Although such regulation in M phase is known (Leitao and Kellogg, 2017; Garmendia-Torres et al., 2018), our data suggest a sequential mechanism to match size to growth rate, with a nutrient-independent sizer followed by a nutrient-dependent timer. To detect the peak in bud growth generated by the sizer, cells may use Gin4-related kinases (Jasani et al., 2020).

An important advantage of the BABY algorithm is that we can estimate single-cell growth rates without fluorescence markers, freeing fluorescence channels for other reporters. Here we focus on Sfp1, a transcription factor that helps coordinate ribosome synthesis with the availability of nutrients (Jorgensen et al., 2004).

Sfp1 promotes the synthesis of ribosomes by activating the ribosomal protein (RP) and ribosome biogenesis (RiBi) genes (Jorgensen et al., 2004; Albert et al., 2019). Upon being phosphorylated directly by TORC1 and likely protein kinase A (Jorgensen et al., 2004; Lempiäinen et al., 2009; Singh and Tyers, 2009) – two conserved nutrient-sensing kinases, Sfp1 enters the nucleus (Figure 5a). In steady-state conditions, levels of ribosomes positively correlate with growth rate (Metzl-Raz et al., 2017), and we therefore assessed whether Sfp1’s nuclear localisation predicts changes in instantaneous single-cell growth rates.

Figure 5 with 1 supplement see all

Download asset Open asset

Shifting cells from glucose to the poorer carbon source palatinose and back again, we observed that Sfp1 responds quickly to both the up- and downshifts and that growth rate responds as quickly to downshifts, but more slowly to upshifts (Figure 5b). As a target of TORC1 and PKA, Sfp1 acts as a fast read-out of the cell’s sensing of a change in nutrients (Granados et al., 2018). In contrast, synthesising more ribosomes is likely to be slower and explains the lag in growth rate after the upshift. The fast drop in growth rate in downshifts is more consistent, however, with cells deactivating ribosomes, rather than regulating their numbers. Measuring the half-times of these responses (Figure 5b boxplots), there is a mean delay of 30 ± 2 minutes (95% confidence; $n=245$) from Sfp1 localising in the nucleus to the rise in growth rate in the upshift. This delay is only 8 ± 1 minutes (95% confidence; $n=336$) in the downshift, and downshift half-times are less variable than those for upshifts, consistent with fast post-translational regulation. Although changes in Sfp1 consistently precede those in growth rate, the higher variability in half-times for the growth rate is not explained by Sfp1’s half-time (Pearson correlation 0.03, $p=0.6$).

By enabling both single-cell fluorescence and growth rates to be measured, BABY permits correlation analyses (Kiviet et al., 2014; Appendix 7). Both Sfp1’s activity and the growth rate vary during the cell cycle. The autocorrelation functions for nuclear Sfp1 and for the growth rate are periodic with periods consistent with cell-division times (Figure 5c): around 90 min in glucose and 140 min in palatinose for Sfp1; and 95 min and 150 min for the growth rate. If Sfp1 acts upstream of growth rate, then its fluctuations in nuclear localisation should precede fluctuations in growth rate. Cross-correlating nuclear Sfp1 with growth rate shows that fluctuations in Sfp1 do lead those in growth rate, by an average of 25 min in glucose and by 50 min in palatinose (Figure 5d). Nevertheless, the weak strength of this correlation suggests substantial control besides Sfp1.

During the downshift, we note that the growth rate transiently drops to zero (Figure 5b), irrespective of a cell’s stage in the cell cycle (Figure 5—figure supplement 1), and there is a coincident rise in the fraction of cells in G1 (Figure 5b bottom), suggesting that cells arrest in that phase.

With BABY, we can use growth rate as a control variable in real time because BABY’s speed and accuracy enables cells to be identified in images and their growth rates estimated during an experiment (Figure 6a). As an example, we switched the medium to a poorer carbon source and used BABY to determine how long to keep cells in this medium if we want 50% to have resumed dividing before switching back to the richer medium (Appendix 8). After 5 hr in glucose, we switched the carbon source to ethanol, or galactose – Figure 6—figure supplement 1. There is a lag in growth as cells adapt. Using BABY, we automatically determined the fraction of cells that have escaped the lag at each time point — those cells that have at least one bud or daughter whose growth rate exceeds a threshold (Figure 6b). The software running the microscopes reads this statistic and triggers the switch back to glucose when 50% of the cells have escaped (Figure 6c). We note that all cells resume dividing in glucose and initially grow synchronously because of the rapid change of media. This synchrony is most obvious in those cells that did not divide in ethanol (Figure 6c).

Figure 6 with 1 supplement see all

Download asset Open asset

This proof-of-principle shows that BABY is applicable for more complex feedback control, where a desired response is achieved by comparing behaviour with a computational model to predict the necessary inputs, such as changes in media (Harrigan et al., 2018; Milias-Argeitis et al., 2011; Toettcher et al., 2011; Uhlendorf et al., 2012; Lugagne et al., 2017; Menolascina et al., 2014). Unlike previous approaches though, which typically measure fluorescence, BABY not only allows single-cell fluorescence but also growth rates to be control variables, and growth rate correlates strongly with fitness (Orr, 2009).

We inoculated overnight cultures with low cell numbers so that they would reach mid-log phase in 13–16 hr. We diluted cells in fresh medium to OD₆₀₀ of 0.1 and incubated for an additional 3–4 hr before loading them into microfluidic devices at ODs of 0.3–0.4. To expose multiple strains to the same environmental conditions and to optimise data acquisition, we use multi-chamber versions of ALCATRAS (Crane et al., 2014; Granados et al., 2017; Crane et al., 2019), which allow for either three or five different strains to be observed in separate chambers while being exposed to the same extracellular medium. The ALCATRAS chambers were pre-filled with growth medium with added 0.05% bovine serum albumin (BSA) to facilitate cell loading and reduce clumping. We passed all microfluidics media through 0.2 μm filters before use.

We captured images on a Nikon Ti-E microscope using a 60×, 1.4 NA oil immersion objective (Nikon), OptoLED light source (Cairn Research) and sCMOS (Prime95B), or EMCCD (Evolve) cameras (both Photometrics) controlled through custom MATLAB software using Micro-manager (Edelstein et al., 2014). We acquired bright-field and fluorescence images at five Z sections spaced 0.6 μm apart. A custom-made incubation chamber (Okolabs) maintained the microscope and syringe pumps containing media at 30 °C.

For experiments in which the cells experience a change of media, two syringes (BD Emerald, 10 ml) mounted in syringe pumps (Aladdin NE-1002X, National Instruments) connected via PTFE tubing (Smiths Medical) to a sterile metal T-junction delivered media through the T-junction and via PTFE tubing to the microfluidic device. Initially the syringe with the first medium infused at 4 μL/min while the second pump was off. To remove back pressure and achieve a rapid switch, we infused medium at 150 μL/min for 40 s from the second pump while the first withdrew at the same rate. The second pump was then set to infuse at 4 μL/min and the first switched off. We reversed this sequence to achieve a second switch in some experiments. Custom Matlab software, via RS232 serial ports, controlled the flow rates and directions of the pumps.

We evaluated BABY’s segmentation on the training data’s test set and compared with recent algorithms for processing yeast images (Padovani et al., 2022): Cellpose version 2.1.1 (Stringer et al., 2021; Pachitariu and Stringer, 2022), YeaZ (Dietler et al., 2020) from 11 October 2022, and our previous segmentation algorithm DISCO (Bakker et al., 2018). For Cellpose and YeaZ, we also trained new models on the images and annotations from both our training and validation sets, following their suggested methods (Pachitariu and Stringer, 2022; Dietler et al., 2020). Because neither handles overlapping regions, we applied a random order to the cell annotations such that pixels in regions of overlap were assigned the last observed label. We augmented the input data for each model by resampling the images five times, thus avoiding bias by forcing the models to adapt to uncertainty in the regions of overlap.

We assessed performance by calculating the intersection over union (IoU) of all predicted masks with the manually curated ground-truth masks from our test set. We paired predicted masks with the ground truth masks beginning with the highest IoU score; we assigned unmatched predictions an IoU of zero. To calculate the average precision for each annotated image, we used the area under the precision-recall curve for varying thresholds on the IoU score (Manning et al., 2008). Not all of the algorithms we tested give a confidence score, and so we generated precision-recall curves assuming ideal ordering of the predicted masks, by decreasing IoU. For the BABY models, ordering by mask probability produces similar results. We report the mean average precision over all images in the test set. To evaluate segmentation on microcolony images, we performed a similar analysis using the ground-truth annotations of the YeaZ bright-field training data (Dietler et al., 2020), but excluding the 18 images annotated and used to train BABY. We also re-trained the Cellpose and YeaZ models using our training data set supplemented with the microcolony images and evaluated the pre-trained bright-field YeaZ model, which includes this evaluation data in its training set, and the general-purpose pre-trained cyto2 Cellpose model, which segments cells from multiple different organisms.

We evaluated tracking on independent, manually curated data, comprising time series with 180–300 time points for 10 randomly selected traps from two experiments and four different growth conditions, making a total of 128 tracks. We initially generated the annotations using an early version of our segmentation and tracking models, but we manually corrected all tracking and lineage assignment errors and any obviously misshapen, misplaced or missing outlines, including removing false positives and adding outlines to the first visible appearance of buds. Unedited outlines, however, remain and will inevitably impart a bias. By requiring a mask IoU score of 0.5 or higher to match masks for the tracking, we expect to negate this bias. We compared BABY with YeaZ (Dietler et al., 2020) and btrack (Ulicna et al., 2021) because Cellpose cannot track. For YeaZ, we used the model trained on our data; for btrack, we used the Cell-ACDC platform (Padovani et al., 2022) to combine segmentation by Cellpose with tracking by btrack.

The output of each model comprises masks with associated labels. We matched predicted and ground-truth masks at each time point to obtain maps from predicted to ground-truth labels, in descending order of mask IoUs but providing the mask IoU was greater than 0.5. We then calculated a track IoU between all predicted and ground-truth tracks: the number of time points where a predicted label mapped to a ground-truth label divided by the number of time points for which either track had a mask. This approach gave a map between predicted and ground-truth tracks in descending order of track IoUs. Using the mapping, we reported either the fraction of predicted tracks whose duration, the number of masks identified within that track, matched the ground-truth tracks (Figure 3c) or the distribution of track IoUs for all ground-truth tracks (Figure 3—figure supplement 4). For the Multiple Object Tracking Accuracy (MOTA) metric (Bernardin et al., 2006), we used the mask IoU to measure distance and considered correspondences as valid if the mask IoU ≥ 0.5.

We evaluated BABY’s lineage assignment using the lineage annotations included in the tracking evaluation data. These assignments pair bud and mother track labels. We used the track IoU to match ground-truth and predicted tracks above a given track IoU threshold and then compared lineage assignments based on this map. We counted true positives for each ground-truth bud-to-mother mapping if the ground-truth bud track had a matching predicted track and this predicted track had a predicted mother track matching the ground-truth mother track. False negatives were any ground-truth mother-bud pairs not counted as true positives; false positives were any predicted mother-bud pairs that were not counted as true positives. We repeated this analysis only for buds assigned to the central trapped cell or its matching predicted track.

We evaluated how well BABY estimates growth rates on independent, manually curated data comprising annotated time series of mother-bud pairs. We did not include this image data, which has growth in glucose, raffinose, pyruvate, and palatinose, in our training data. To select positions, traps, and time points, we randomly selected mother-bud pairs, rejecting samples only if there was no pair with a complete bud-to-bud cycle. We segmented this data with BABY and Cellpose and YeaZ trained on our data. To avoid penalising YeaZ and Cellpose for tracking errors, we found the matching predicted outlines with highest positive IoU for each ground-truth mask. We then used our method to estimate volumes (Appendix 3) to derive volumes for all masks, both ground-truth and predicted. Associating the masks with the ground-truth track, we fit a Gaussian process to each time series of volumes, omitting any time points with no matching mask. From the Gaussian process, we estimated a growth rate for each time point. Finally, we calculated the Root Mean Square Distance (RMSD) between the predicted and ground-truth estimates.

For epifluorescence microscopy, samples are typically prepared to constrain cells in a monolayer. For cells with similar sizes that match the height of this constraint, they will be physically prevented from overlapping. If cells are of different sizes, however, then a small cell can potentially fit in gaps and overlap with others. This phenomenon is especially prevalent for cells that divide asymmetrically, where a small bud grows out of a larger mother.

Few segmentation algorithms identify instances of overlapping cells. Most, including recent methods for budding yeast (Wood and Doncic, 2019; Dietler et al., 2020; Lugagne et al., 2018), assume that cells can be labelled semantically, with each pixel of the image identified with at most one cell. Similarly, most tools for annotating also label semantically, and consequently curated training data does not allow for overlaps (Dietler et al., 2020), even when the segmentation algorithm could (Lu et al., 2019). Our laboratory’s previous segmentation algorithm included limited overlap between neighbouring cells (Bakker et al., 2018), but not the substantial overlap we see between the smaller buds and their neighbours.

We rely on two consequences of the height constraint to segment overlapping instances. First, cells of different sizes show different patterns of overlap; second, the cells’ centres are rarely coincident. Very occasionally, we do observe small buds stacked directly on top of each other, but neglecting these rare cases does not degrade performance. We therefore use morphological erosions to obtain semantic images by shrinking cell masks within a size category and, later, morphological dilations to approximate the original cell outlines from each resulting connected region.

To separate overlapping cells, we define three size categories and treat instances in each category differently. Appendix 1—figure 1 illustrates our approach, where we segment a bud (orange outline) that overlaps a mother cell (green outline). The bud is only visible in the third and fourth Z sections of the bright-field images (Appendix 1—figure 1a). If used for training, we would split the manually curated outlines in this example (Appendix 1—figure 1b) into different size categories (Appendix 1—figure 1c). The bud is assigned to the small category. When we fill the outlines in this category and convert the image to a binary one (Appendix 1—figure 1d), the individual cell masks are distinct. For the large category, however, the masks are not separable when immediately converted, but become so when the filled outlines are morphologically eroded (Appendix 1—figure 1d). The largest size category tolerates more erosions than smaller ones, for which the mask may disappear or lose its shape.

Appendix 1—figure 1

Download asset Open asset

Using the training data – curated masks for each cell present at each trap at each time point, we identify the size categories that best separate overlapping cells. To begin, we calculate the overlap fraction – the intersection over union – between all pairs of cell masks. Its distribution reveals that the most substantial overlaps occur between cells of different sizes (Appendix 1—figure 2a – upper triangle).

Appendix 1—figure 2

Download asset Open asset

We therefore choose the size categories so that most overlaps occur between pairs of cells in different categories and little overlap occurs between pairs of cells within a category. For example, rather than converting the cell masks directly into a single binary image for training, if first we divide cells into two size categories and convert the masks within each category to a separate binary image, giving two images rather than one, then in these two images we have eliminated all overlaps occurring between cells in the smaller category with cells in the larger category (Appendix 1—figure 1 and Appendix 1—figure 2 – lower triangle).

To divide the cell masks into two categories, we define a fuzzy size threshold using a threshold $T$ and padding value $P$. The set of smaller masks is all masks whose area is less than $T+P$; the set of larger masks is all masks whose area is greater than $T-P$. Consequently, the same mask can be in both sets (Appendix 1—figure 1c). This redundancy ensures the CNN produces confident predictions even for cells close to the size threshold, and we eliminate any resulting duplicate predictions in post-processing. BABY prevents a pair of masks overlapping by converted each into distinct binary images if the padded threshold separates their sizes: the smaller cell must have a size ${\displaystyle <T-P}$ and the larger cell must have a size ${\displaystyle >T+P}$. To scale with pixel size, we set $P$ to be 1% of the area of the largest mask in the training set.

To determine an optimal fuzzy threshold, we test $B=100$ values evenly spaced between the minimal and maximal mask sizes and choose the threshold that minimises the summed overlap fraction for all mask pairs not excluded by the threshold. Even with one fuzzy threshold (Appendix 1—figure 2a), we exclude most of the pairs with substantial overlap – typically buds with neighbouring cells (Appendix 1—figure 2c).

After applying the threshold, overlaps between cells within a size category remain, and we reduce such overlaps using morphological erosions (Appendix 1—figure 1). We use the training data to optimise the number of erosions per size category. As the number of iterations increases, there is a trade-off between the number of overlapping mask pairs and the number of masks whose eroded areas become too small to be confidently predicted by the CNN (Appendix 1—figure 2d). Without erosion, the large cells can show overlaps; with too much erosion, the smallest masks distort their shapes or disappear. We therefore optimise the number of iterations separately for each size category, picking the highest number of iterations that do not let any of that category’s training masks either fall below an absolute minimal size, defined as 10 pixels squared, or fall below 50% of the category’s median cell size before any erosions.

Combining categorising by size with eroding reduces the number of pairs of overlapping masks almost to zero (Appendix 1—figure 2e). We arrive at three size categories by first introducing an additional fuzzy threshold for each of the two initial size categories. These thresholds are similarly determined by testing $B=100$ fuzzy threshold values and calculating the overlap fraction for all mask pairs not excluded by either the original or the new threshold. We only keep one of the new thresholds – the one minimising the overlap fraction, giving three size categories in total. This extra category results in a further, although proportionally smaller, decrease in the number of overlapping masks.

After erosion, mask interiors within each size category are easily identified, but with less resolved edges. To help alleviate this loss, we generate edge targets for the CNN from the training data (Appendix 1—figure 1e) – the outlines of all cells within each size category.

The microcolony training images for YeaZ (Dietler et al., 2020) include a larger range of cell sizes than in our training set. We therefore increased $B$ to 200 (Figure 3—figure supplement 2) and determined the thresholds on square-root transformed sizes. We transformed these thresholds back to the original scale when providing targets for the CNN.

We further annotate the curated data with lineage assignments (Appendix 1—figure 1f), which BABY uses to generate ‘bud neck’ targets for the CNN (Appendix 1—figure 1g). The final target is another binary image, which is only true wherever masks of any size overlap.

In total, the eight training targets for the CNN are the mask interiors and edges for three size categories, the bud necks, and the overlap target. We weighted the targets according to their difficulty and importance in post-processing steps: the large and medium edge targets and small interior target with a weight of two and the small edge target with one of four.

We trained fully convolutional neural networks (Goodfellow et al., 2016) to map a stack of bright-field sections to multiple binary target images. We show some example inputs and outputs in Appendix 1—figure 1, but we also trained networks with only one or three bright-field sections. The intensities of the bright-field sections were normalised to the interval $[-1,1]$ by subtracting the median and scaling according to the range of intensities expected between the 2nd and 98 percentiles.

Each output layer of the CNN approximates the probability that a given pixel belongs to the target class, being a convolution with kernel of size 1 × 1 and sigmoidal activation. All other convolutions had kernels of size 3 × 3 with ReLU activation and used padding to ensure consistent dimensions for input and output layers.

To prevent over-fitting and improve generalisation, we augmented the training data (Goodfellow et al., 2016). Each time the CNN sees a training example, it sees too a randomly selected series of image manipulations applied to the input and target. The same training example therefore typically appears differently for each epoch.

Three augmentations were always applied and the others applied with a certain probability. The fixed augmentations were horizontal and vertical translations and if the bright-field input had more Z sections than expected by the network, we selected a random subset, excluding any subsets with selected sections separated by two or more missing sections. Those augmentations applied with a probability $p$ comprised elastic deformation ($p=0.35$), image rotation ($p=0.35$), re-scaling ($p=0.35$), vertical and horizontal flips (each with $p=0.35$), addition of white noise ($p=0.35$), and a step shift of the Z sections ($p=0.35$). The probability of not augmenting was thus $p=0.05$. To show a different region of each image-mask pair at each epoch, translation, rotation, and re-scaling were all applied to images and masks before cropping to a consistent size (128 × 128 pixels for a pixel size of 0.182μm). Using reflection to handle the boundary, translations were over a random distance and rotations over a random angle. To apply elastic deformation, as described for the original U-Net (Ronneberger et al., 2015b), we used the elasticdeform package (van Tulder, 2022) for an evenly spaced grid with target distance between points of 32 pixels and standard deviation of displacement of 2. Augmentation by re-scaling was for a randomly selected scaling factor up to 5%. Augmentation by addition of white noise involved adding random Gaussian noise with a standard deviation picked from an exponential distribution with rate to each pixel of the (normalised) bright-field images.

To reduce aliasing errors when manipulating binary masks during augmentation, we applied all image transformations independently to each filled mask before converting the transformed masks into one binary image. Further, before a transformation, we smoothed each binary filled outline with a 2D Gaussian filter and found the transformed binary outline with the Canny algorithm. To determine the standard deviation of this Gaussian filter, σ, we tested a range of values on the training outlines. For each filled outline and σ, we applied the filter followed by edge detection and filling. We then calculated the intersection over union of the resulting filled outline with the original filled outline. We observed that as a function of edge length, defined as the number of edge pixels, the σ producing the highest intersection over union increased exponentially. We consequently used an exponential fit of this data to estimate an appropriate σ for each outline.

Appendix 1—figure 3

Download asset Open asset

We trained networfks using Keras with TensorFlow 2.8. We used Adam optimisation with the default parameters except for a learning rate of 0.001 and regularised by keeping only the network weights from the epoch with the lowest validation loss (similar in principle to the early stopping method) (Goodfellow et al., 2016). We train for 600 epochs, or complete iterations over the training data set.

The loss function is the sum of the binary cross-entropy and one minus the Dice coefficient across all targets:

where $y$ is the tensor of true values, $\widehat{y}$ is the CNN’s sigmoid tensor output of the CNN, and $i$ is a vectorised index.

Each CNN is trained to a specific pixel size, and we ensured that training images and masks with different pixel sizes were re-scaled appropriately

We trialled two core architectures for the CNN – U-Net (Ronneberger et al., 2015b; Appendix 1—figure 3a) and Mixed-Scale-Dense (MSD) (Pelt and Sethian, 2018) – and optimised hyperparameters to find the smallest loss on the validation data.

The U-Net performed best (see ‘Optimising hyperparameters’ below). The U-Net has two parts: an encoder that reduces the input into multiple small features and a decoder that scales these features up into an output (Ronneberger et al., 2015b). Each step of the encoder comprises a convolutional block, which creates a new, larger set of features from its input. To force the network to keep only small, relevant features, a down-sampling step is applied after three convolutional blocks. This maximum pooling layer reduces the size of the features by half by replacing each two-by-two block of pixels by their maximal value. The decoder also comprises convolutional blocks, but with up-sampling instead of down-sampling. The up-sampling step is the inverse of down-sampling: each pixel is turned into a two-by-two block by repeating its value. Finally, most characteristic of the U-Net is its skip layers. These layers preserve information on the global organisation of the pixels by passing larger-scale information from the encoder to the decoder after each up-sampling step. They act by concatenating the same-size layer of the encoder into the decoder layers, which are then used as inputs for the next step of the decoder. The decoder is therefore able to create an output from both the local features that it up-sampled and from the global features that it obtains from the skip layers.

For the U-Net, we optimised for depth, for a scaling factor for the number of filters output by each convolution, whether or not to include batch normalisation, and for the proportion of neurons to drop out on each batch. For the MSD, we optimised for depth, defined as the total number of convolutions, for the number of dilation rates to loop over with each loop increasing dilation by a factor of two, for an overall dilation-rate scaling factor, and whether or not to include batch normalisation.

We used KerasTuner with TensorFlow 2.4 to optimise hyperparameters, choosing random search with default settings, training for a maximum of 100 epochs, and having 10 training and validation steps per epoch. The U-Net and MSD networks with the lowest final validation loss were then re-trained as described, and the network with the lowest validation loss chosen.

For our data, the best performing model was a U-Net with depth four, and so three contractions, with a scaling factor of 16 for the number of filters output by each convolution, giving 16, 32, 64 and 128 filters for each of the two chained convolution layers of the encoding and decoding blocks, with batch normalisation, and with no drop-out. We show its performance for the 5Z model in Appendix 1—figure 3c–e.

To identify cell instances from the semantic predictions of the CNN, we developed a post-processing pipeline with two parts (Appendix 1—figure 4a): proposing unique cell outlines and then refining their edges.

The pipeline includes multiple parameters that we optimise on validation data by a partial grid search. We favour precision, the fraction of true predicted positives, over recall, the fraction of ground truth positives we predict, by maximising the $F_{\beta }$ score with $\beta =0.5$. Recall that for true positives $\mathrm{TP}$, false negatives $\mathrm{FN}$, and false positives $\mathrm{FP}$,

We measure how well two masks match using the intersection over union (IoU) and consider a match to occur if ${\displaystyle \mathrm{I}\mathrm{o}\mathrm{U}>0.5}$. Nevertheless, multiple predictions may match a single target mask because predicted masks can overlap too. We therefore count true positives as target masks for which there is at least one predicted mask with ${\displaystyle \mathrm{I}\mathrm{o}\mathrm{U}>0.5}$. Any predicted masks in excess of the true positive count are false positives, thus avoiding double counting. Unmatched target masks are false negatives.

The post-processing pipeline starts by identifying candidate outlines independently for each size category. The CNN’s outputs are images $p_{xy}^{(S,C)}\in [0,1]$ approximating the probability that a pixel at position $(x,y)$ belongs to either the small, medium, or large size categories, denoted $S$, and to one of the other classes, denoted $C$: either the interior (Appendix 1—figure 4b), edge (Appendix 1—figure 4e and f), bud neck, or general overlap classes.

In principle, we could find instances for each size category by thresholding the interior probability $p_{xy}^{(S,\mathrm{interior})}$ and identifying connected regions as outlines. To further enhance separability, however, we also re-weight the interior probabilities using the edge probabilities. Specifically, we identify connected regions from semantic bit masks $b_{xy}^{(S,\mathrm{interior})}$ by those pixels that satisfy

where ${\mathrm{Dilate}}^N$ specifies $N$ iterations of a gray scale morphological dilation and $T_{\mathrm{interior}}$ is a threshold. We optimise the thresholds $T_{\mathrm{interior}}\in [0.3,0.95]$, number of dilations $N_{\mathrm{dilate}}\in \{0,1,2\}$, and the order of connectivity (one- or two-connectivity) for each size category.

Appendix 1—figure 4

Download asset Open asset

The connected regions in $b_{xy}^{(S,\mathrm{interior})}$ define masks that are initial estimates of the cells’ interiors (darker shading in Appendix 1—figure 4c). We generate the cell interiors for training the CNN by iterative, binary morphological erosions of the full mask, where the number of iterations $N_{\mathrm{erosion}}$ is pre-determined for each size category. First, we remove small holes and small foreground features by applying up to two binary morphological closings followed by up to two binary morphological openings. Second, we estimate full masks $b_{xy}^{(S,\mathrm{full})}$ from each putative mask by applying $N_{\mathrm{erosion}}$ binary dilations (light shading in Appendix 1—figure 4c) undoing the level of erosion on which the CNN was trained. We optimise both the numbers of closing, $N_{\mathrm{closing}}$, and opening, $N_{\mathrm{opening}}$, operations.

Any masks whose area falls outside the limits for a size category, we discard. For each category, however, we soften the limits, on top of the fuzzy thresholds, by optimising an expansion factor $F_{\exp }\in [0,0.4]$, which extends the limits by a fractional amount of that category’s size range. We also optimise a single hard lower threshold $T_{\min }\in [0,20]$ on mask area.

To prepare for refining edges and to further smooth and constrain outlines, we use a radial spline to match the edge of each of the remaining masks (Appendix 1—figure 4d). As in DISCO (Bakker et al., 2018), we define radial splines as periodic cubic B splines using polar coordinates whose origin is at the mask’s centroid. We generalise this representation to have a variable number $n^{(S,i)}$ of knots per mask specified by $n^{(S,i)}$-dimensional vectors of radii ${\displaystyle {\mathit{r}}^{(S,i)}}$ and angles ${\displaystyle {\mathit{\theta }}^{(S,i)}}$:

A mask’s outline is then fully specified by those pixels that intersect with this spline.

To initially place the knots, we search along rays originating at the centroid of each mask $b_{xy}^{(S,\mathrm{full})}$ and find where these rays intersect with the mask edge. We determine the edge by applying a minimum filter with two-connectivity to the mask and set to true all pixels in the filtered image that are different from the original one. We then smooth the resulting edge image using a Gaussian filter with $\sigma =0.5$. For a given polar angle $\theta$, we find the radius of the corresponding knot by averaging the edge pixels that intersect with the ray, weighted by their values. We use the major axis of the ellipse with the same normalised second central moment as the mask (regionprops function from Scikit-image van der Walt et al., 2014) to determine both the number of rays, and so knots, and their orientations. The length ${\mathrm{\ell }}^{(S,i)}$ of the major axis gives the number of rays: four for ${\displaystyle 0<{\ell }^{(S,i)}<5}$; six for ${\displaystyle 5\le {\ell }^{(S,i)}<20}$; and eight for ${\displaystyle {\ell }^{(S,i)}\ge 20}$. For this initial placement, we choose equiangular ${\mathit{\theta }}^{(S,i)}$, with the first knot on the ellipse’s major axis.

The quality of the outline masks ${\widehat{o}}_{xy}^{(S,i)}$ derived from these initial radial splines are then assessed against the edge probabilities generated by the CNN (Appendix 1—figure 4e) and masks of poor quality discarded. We calculate the edge score for a given outline as

We discard those outlines for which the edge score is less than a threshold, where the thresholds $T_{\mathrm{edge}}\in [0,1)$ are optimised for each size category based on the range of edge scores observed.

With a smoothed and filtered set of outlines, we proceed by detecting and eliminating any outlines duplicated between size categories. We start by filling the outlines to form a set of full masks ${\widehat{m}}_{xy}^{(S,i)}$. We then compare these masks between neighbouring size categories $S_j$ and $S_k$. We consider the pair of masks i₁ and i₂ as duplicates if one of the masks is almost wholly contained within the other:

for some threshold ${\displaystyle T_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{t}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{t}}\in [0,1]}$, optimised on validation data. For pairs that exceed this threshold, we keep only the mask with the highest edge score given by Equation 5.

For each size category, the first part of the post-processing pipeline finishes with the set of outlines that pass these size, edge probability, and containment thresholds. Appendix 1—table 1 gives values for the optimised post-processing parameters.

The outlines ${\widehat{o}}_{xy}^{(S,i)}$, defined by the radial splines, do not directly make use of the CNN’s edge targets for their shape and deviate from $p_{xy}^{(S,\mathrm{edge})}$, particularly for those in the large size category (Appendix 1—figure 4e).

We therefore optimise the radial splines to better match the predicted edge. This optimisation is challenging because $p_{xy}^{(S,\mathrm{edge})}$ provide only a semantic representation of the edge – the association of a given pixel $(x,y)$ with a particular instance $i$ is unknown. Our approach is to use the outlines to generate priors on whether predicted edge pixels associate with a given instance. We then apply standard techniques to optimise the fit of the radii and angles of the knots for each outline’s spline to its instance’s likely edge pixels.

Appendix 1—figure 5

Download asset Open asset

To associate pixels with instances, we first calculate the radial distance of each pixel from the initial radial spline function $s$ proposed for an instance in Equation 4. To increase speed, we consider only pixels where ${\displaystyle p_{xy}^{(S,\mathrm{e}\mathrm{d}\mathrm{g}\mathrm{e})}>0.2}$. Expressing the edge pixels in polar coordinates as $(\rho ,\varphi )$ with the origin at the instance’s centroid, this distance is

which we will refer to as a pixel’s residual. We give two examples of edge pixels (Appendix 1—figure 5a) and of the corresponding residuals (Appendix 1—figure 5b), which highlight the need to associate pixels with a given instance before attempting to optimise the spline.

We use the residuals, Equation 7, to assign prior weights to pixels:

where ${\sigma }_G=5$ and ${\sigma }_E=1$. The function $W$ is a Gaussian function of the residual for pixels exterior to the proposed outline and an exponential function for pixels interior (Appendix 1—figure 5b). This asymmetry should increase tolerance for interior edge pixels, which may belong to neighbouring instances overlapping with the cell of interest. In such cases, we should thus improve instance association, particularly where the edges of each of the cells intersect.

With these prior weights, we find the probability that each edge pixel associates with a particular instance and not with the others via:

where $n$ is the number of detected instances in this and adjacent size categories, with $j$ running over all these instances. We filter the result, keeping only those edge pixels with ${\displaystyle p^{(S,\mathrm{e}\mathrm{d}\mathrm{g}\mathrm{e},i)}>0.1}$. Examples are shown in Appendix 1—figure 5c.

We optimise the knot radii ${\displaystyle {\mathit{r}}^{(S,i)}}$ and angles ${\displaystyle {\mathit{\theta }}^{(S,i)}}$ for each radial spline by minimising the squared radial residual between the spline and the edge pixels, Equation 7. With residuals weighted by $p_{xy}^{(S,\mathrm{edge},i)}$, Equation 9, and initial values taken from each ${\widehat{o}}_{xy}^{(S,i)}$, we constrain radii to a 30% change from their initial values and angles to a change of ±25% of the initial angular separation between knots: ${\theta }_{i+1}-{\theta }_i$. The resulting optimised radial splines provide the outlines output by the BABY algorithm.

In our setup, daughter cells may be washed out of the microfluidic device and so disappear from one time point to the next. These absences undermine other approaches to tracking, such as the Hungarian algorithm (Versari et al., 2017).

To track cells (Appendix 2—figure 2), we use the changes in their masks over time to indicate identity. From each mask, we extract an array of attributes, such as the mask’s area, major axis length, etc., and to compare a mask at one time point to a mask at another time point, we subtract the two corresponding arrays of features. This array of differences is the array of features we use for classification.

Appendix 2—figure 2

Download asset Open asset

Our training data comprises a series of manually labelled time-lapse images from four experiments. For two consecutive time points, we calculated the difference in feature arrays between all pairs of cells and grouped these difference arrays into two classes: one for identical cells – cells with the same label – and one for all other cells.

To generate additional training data, we use multiple time points backwards in time. For example, for time $t$, we generate not only feature vectors by comparing with cells at $t-1$, but also with cells at $t-2$ and $t-3$. We found this additional data increased generalisability, maintaining accuracy across a wider range of imaging intervals and growth rates. For the purpose of training, we treat the additional data as consecutive time points: the algorithm does not know whether the features come from one or more than one time point in the past.

As part of testing if all features contribute to the learning, we divided the features into two overlapping sets. One set had no features that explicitly depend on distance, comprising area, lengths of the minor and major axes, convex area, and area of the bounding box; the second set did include distance-dependent features, comprising area, lengths of the minor and major axes, and convex area again, but additionally including the mask’s centroid, and the distance obtained from the x- and y-axis locations.

We compared three standard algorithms for classification (Bishop, 2006): the Support Vector Classifier (SVC), Random Forest, and Gradient Boosting, specifically Xtreme Gradient Boosting (Chen and Guestrin, 2016). We used scikit-learn (Pedregosa, 2011) to optimise over a grid of hyperparameters.

For the SVC, we considered a regularisation parameter $C$ of 0,1, 10, or 100; a $\mathrm{\Gamma }$ kernel coefficient of 1, 10−3, or 10−4; no shrinking heuristic to speed up training; and either a radial basis function or sigmoid kernel.

For the Random Forest, we explored a range between 10 and 80 estimators and a depth between 2 and 10 levels.

For Gradient Boosting, we used a maximal depth of either 2, 4, or 8 levels; a minimal child weight of 1, 5, or 10; gamma, the minimal reduction in loss to partition a leaf node, of 0.5, 1, 1.5, 2, or 5; and a sub-sampling ratio of 0.6, 0.8, or 1.

Within the training data, the number of time points for each experiment is different. To prevent biases toward long experiments, we define the accuracy as the fraction of true positives – cells correctly linked between images – and compare the precision and recall of this time-averaged accuracy.

After training, we evaluated which features were important using the Random Forest. The distribution of the feature weights depends on whether we include distance (Appendix 2—figure 3), and excluding distance distributes the weights more evenly.

Appendix 2—figure 3

Download asset Open asset

The precision-recall curve indicates that using the distance-explicit features is best, although both sets of features have high accuracy (Appendix 2—figure 3c). Despite performing better on our test data, we expect that using the distance-explicit features may perform worse if the cells pivot or become displaced. Therefore, we use the non-explicit features as our main model, but also use the distance-explicit features to resolve any ambiguous predictions. The ensemble model performs similarly to the distance-implicit classifier, but for more stringent thresholds behaves like the distance-explicit one.

To predict with the classifier, we use data from the current time point and the two most recent previous time points. We generate feature arrays between $t$ and independently $t-1$ and $t-2$ and feed both arrays to the classifier. If the probability returned is greater than 0.9, we accept the result; if the probability lies between 0.1 and 0.9, we use instead the probability returned by the backup classifier, which uses the distance-explicit features.

Using multiple time points to track cells has two advantages: first, it reduces noise generated by artefacts, either in image acquisition, such as a loss of focus, or in segmentation; second, it ensures that cells are more consistently identified if their position or shape transiently changes. Including data further back in time is neither computationally efficient nor more accurate, and greater than three time points is long, over 15 minutes in our experiments and about a sixth of a cell cycle.

We apply the linear sum assignment algorithm, via SciPy, on the probability matrix of predictions to assign labels (Appendix 2 Algorithm 1). This approach guarantees at most one outline is assigned to each cell by choosing the set of probabilities whose total sum is highest. To match a cell with its previous self, we pick the cell in the recent past that generates the highest probability when paired with the cell of interest, providing this probability is greater than 0.5. We label a cell as new if the probabilities returned from pairing with all cells in the recent past is below 0.5.

We wish to identify which cells are buds of mothers and which mothers have buds. This problem is analogous to tracking, but, rather than identifying pairs of cells that are the same cell at different time points, we must identify pairs of cells that are a mother-bud pair at one time point. We therefore seek to determine the probability that a pair of cells is a mother-bud pair (Appendix 2—figure 4). Unlike tracking, however, we anticipated that the cell outlines alone would be at best a weak indicator of this probability.

We observed that cytokinesis is sometimes visible in bright-field images as a darkening of the bud neck and designed features to exploit this characteristic of mother-bud pairs.

Such features often rely on the CNN’s prediction of bud necks. For generalisability and to avoid ambiguity, we chose to define the corresponding training target using manually annotated outlines and lineage relationships, rather than relying on a fluorescent bud-neck marker. Specifically, we define a binary semantic ‘bud-neck’ training target that is true only at pixels where a mother mask, dilated twice by morphological dilation, intersects with its assigned bud (Appendix 1—figure 1). Assigning a time of cytokinesis by eye is challenging, and so we included two constraints to identify a bud. First, the bud must be current – as soon as BABY finds another bud associated to the mother, we drop the current one. Second, we exclude buds if their area is larger than and has always been larger than 10 μm² (300 pixels for our standard training target with a pixel size of 0.182 μm and corresponding to a sphere of ∼ 24 μm³).

Appendix 2—figure 4

Download asset Open asset

We used multiple image features to characterise a mother-bud relationship. For an ordered pairing of all cells in an image, we consider a putative mother-bud pair and define a mask, $M_{\mathrm{joining}}$, as the joining rectangle between the centres of the mother and bud with a width equal to one quarter the length of the bud’s minor axis. Given $M_{\mathrm{joining}}$, we consider five features:

1. $F_{\mathrm{size}}$, which is the ratio of the mother’s to bud’s area. Mothers generally have a greater size than their bud so that ${\displaystyle F_{\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e}}>1}$.

2. $F_{\mathrm{adjacency}}$, which is the fraction of $M_{\mathrm{joining}}$ intersecting with the union of the mother’s and bud’s masks. Mothers should be proximal to their buds so that $F_{\mathrm{adjacency}}$ is close to one: only a small fraction of $M_{\mathrm{joining}}$ should lie outside of the mother and bud outlines.

3. $F_{\mathrm{bud}}$, which is the mean over the union of the CNN’s output for a small, interior targets and all pixels contained in the bud. $F_{\mathrm{bud}}$ approximates the probability that a cell is a bud and should be close to one for mother-bud pairs.

4. $F_{\mathrm{p}}$, which is the mean over the union of the pixels contained in $M_{\mathrm{joining}}$ with the CNN’s output for bud-necks, only including those pixels whose probability is greater than 0.2. $F_{\mathrm{p}}$ approximates the probability that a bud neck joins a mother and bud.

5. $F_{\mathrm{c}}$, which is the number of the CNN’s bud-neck target pixels with a probability greater than 0.2 that are in $M_{\mathrm{joining}}$ normalised by the square root of the bud’s area, or effectively the bud’s perimeter. We interpret $F_{\mathrm{c}}$ as a confidence score on $F_{\mathrm{p}}$ because a single spurious pixel with high bud-neck probability could produce high $F_{\mathrm{p}}$.

With these features, we train a random forest classifier to predict the probability that a pair of cells is a mother and bud. We train on all pairs of cells in the validation data. We optimised the hyperparameters, including the number of estimators and tree depth, using a grid search with five-fold cross-validation. We optimise for precision because true mother-bud pairs are in the minority and because our strategy for assigning lineages aggregates over multiple time points (as detailed below).

For our standard 5Z CNN trained to a pixel size of 0.182, the random forest classifier had a precision of 0.83 and recall of 0.54 on the test data. This precision and recall of the classifier to assign mother-bud pairs within a single time point is distinct from the precision and recall of mother-bud pairs after accumulating across time (reported in Figure 3d). This data has 211 true mother-bud pairs out of 1678 total pairs. The classifier assigned feature weights of 0.46 for $F_{\mathrm{size}}$, 0.11 for $F_{\mathrm{adjacency}}$, 0.24 for $F_{\mathrm{bud}}$, 0.06 for $F_{\mathrm{p}}$, and 0.13 for $F_{\mathrm{c}}$.

To establish lineages, we need to assign at most one tracked mother cell to each tracked cell object. We use the classifier to assign a mother-bud probability $p_{mb}^{(t)}$ for each time point at which a pair of tracked objects are together. We then estimate the probability that a tracked object $\widehat{c}$ has ever been a mother using

as well as the probability that it has ever been a bud with

Finally, we calculate a cumulative score for a putative mother-bud pair and reduce this score if the candidate bud has previously shown a high probability of being a mother:

At each time point, we then propose lineages by assigning each putative cell object ${\widehat{c}}_b$ with a bud probability ${\displaystyle p_b^{(t)}\left({\hat{c}}_b\right)>0.5}$ and a mother ${\widehat{c}}_m={\mathrm{argmax}}_{\widehat{c}}\left(S_{mb}^{(t)}(\widehat{c},{\widehat{c}}_b)\right)$. We treat the mother-bud assignments proposed at the final time point as definitive because they have integrated information over the entire time series. To avoid spurious assignments, we require all buds to be present for at least three time points.

Though rare, we do have to mitigate occasional detection, tracking and assignment errors. For example, debris can occasionally be mistakenly identified as a cell and tracked.

We discard tracks that have both small volumes and show limited growth over the experiment. Specifically, we discard a given cell track $i$ with duration $T_i$, minimal volume $V_i^{(\min )}$ at time $T_i^{(\min )}$, and maximal volume $V_i^{(\max )}$ at time $T_i^{(\max )}$, if both ${\displaystyle V_i^{(\mathrm{m}\mathrm{a}\mathrm{x})}<7}$ μm³ and the estimated average growth rate ${\displaystyle G_i<10}$ μm³/hour, where

Our tracking algorithm usually identifies correctly instances where a mother and bud pivot with the flow of the medium, but exceptions do arise. For a given mother, we therefore join contiguous bud tracks – pairs of bud tracks where one ends with the other starting on the next time point – if the extrapolated volume of the old track falls within a threshold difference of the volume of the new track. Specifically, for the pair of contiguous tracks $i$ and $j$, with track $i$ ending at time point $t$ and track $j$ beginning at time point $t+1$, we calculate

where $V_i^{(t)}$ is the volume of track $i$ at time point $t$ and $\mathrm{\Delta }T$ is the time step between time points $t$ and $t+1$. We join these tracks if ${\displaystyle V_{ij}^{(\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f})}<7}$ μm³.

Finally, we discard any tracks with fewer than five time points.

To estimate cell volumes, we model a 3D cell from our 2D outline. We use a conical method (Gordon et al., 2007), which is robust to common cell shapes, to estimate cell volume from an outline. This method makes two assumptions: that the outline obtained cuts through the widest part of the cell and that the cell is ellipsoidal. We build a cone with a base shape that is the filled outline of the cell by iteratively eroding the segmented mask of the cell and stacking these masks in the $Z$ dimension. We find the volume of the cone by summing the voxels in the corresponding 3D mask. Finally, we multiply this sum by four to obtain the volume of the cell: a cone whose base is the equatorial plane of an ellipsoid will have a volume that is a quarter of the corresponding ellipsoid’s volume (Gordon et al., 2007).

Depending on the need for computational speed, we use one of two methods for estimating instantaneous growth rates.

For long-term, and stored, analysis, we estimate growth rates by fitting a Gaussian process with a Matern covariance function to the time series of each cell’s volume (Swain et al., 2016). We set the bounds on the hyperparameters to prevent over-fitting. Maximising the likelihood of the hyperparameters, we are able to obtain the mean and first and second time derivatives of the volume, as well as estimates of their errors. The volume’s first derivative is the single-cell growth rate.

During real-time processing where we may wish to use to the growth rate to control the microscope, fitting a Gaussian process is too slow. Instead we estimate growth rates from the smoothed first derivative obtained by Savitzky-Golay filtering of each cell’s volume time series. Though faster, this method is less reliable and does not estimate errors. For time series of mothers, we use a third-order polynomial with a smoothing window of seven time points; for time series of buds, we use a fourth-order polynomial also with a smoothing window of seven time points.

We estimate growth rates separately for mothers and their buds because both are informative. We find that the summed results are qualitatively similar to previous estimates of growth rate, which fit the time series of the combined volume of the mother and its bud (Ferrezuelo et al., 2012; Cookson et al., 2010).

• Peter S Swain
• Julian MJ Pietsch

• Peter S Swain
• Ivan BN Clark
• Iseabail Farquhar

• Alán F Muñoz
• Diane-Yayra A Adjavon