Efficient microbial colony growth dynamics quantification with ColTapp, an automated image analysis application

Populations of genetically identical bacteria are phenotypically heterogeneous, giving rise to population functionalities that would not be possible in homogeneous populations. For instance, a proportion of non-dividing bacteria could persist through antibiotic challenges and secure population survival. This heterogeneity can be studied in complex environmental or clinical samples by spreading the bacteria on agar plates and monitoring time to growth resumption in order to infer their metabolic state distribution. We present ColTapp, the Colony Time-lapse application for bacterial colony growth quantification. Its intuitive graphical user interface allows users to analyze time-lapse images of agar plates to monitor size, color and morphology of colonies. Additionally, images at isolated timepoints can be used to estimate lag time. Using ColTapp, we analyze a dataset of Staphylococcus aureus time-lapse images including populations with heterogeneous lag time. Colonies on dense plates reach saturation early, leading to overestimation of lag time from isolated images. We show that this bias can be corrected by taking into account the area available to each colony on the plate. We envision that in clinical settings, improved analysis of colony growth dynamics may help treatment decisions oriented towards personalized antibiotic therapies.

Extended description of the implementation 1

Graphical user interface
Global operating buttons to load images and save data are in the top left part of the graphical interface ( Fig. S1-a) and information about currently operating functions and text feedback in the top right part (Fig. S1-b). Upon loading a folder containing images to analyze, the user needs to define the operating mode: Time-Lapse (TL) or Endpoint (EP) mode, which cannot be changed for further analysis. Images in a folder are sorted with a natural sorting algorithm by name and assigned a frame number (from 1 to number of images in the folder). The currently active frame is visible at the center of the graphical user interface (Fig. S1-c). The user can switch between frames by moving the slider, typing the number of the desired frame ( Fig. S1-

Spatial calibration
The program proposes to define a pixel to µm spatial calibration factor. The user may either directly enter the factor if known, or define the diameter of a circle, typically the agar plate (automatically detected), or any reference line of known length on the image.

Area of interest
Automatic colony detection can be narrowed to the defined boundaries of an area of interest (AOI). Users can select the plate as AOI or draw a custom polygon on the image. Using an AOI is suggested to reduce computational time and reduce false positive detection outside the boundaries of e.g. an agar plate. In time lapse mode, both the spatial calibration factor and AOI are propagated from the active frame to the other frames. In EP mode, these variables are only applied to the current image. This can then be manually propagated if all images in a given folder were acquired with the same setup in Options: Apply calibration factor/AOI to all frames.

Color image conversion
ColTapp offers 16 methods to transform a color image into grayscale. Using one of the three RGB channels is computationally quick and generally precise enough, but in some special cases the user may want to use a different conversion as for example a transformation to other classical color spaces (e.g. CIELAB) and using one of the generated channels. The inbuilt MATLAB function rgb2gray which uses a weighted image conversion, retaining information of all three color channels, is also available.

Algorithm parameters
A selection of parameters involved in circle detection and quality control can be tuned to potentially improve performance, although we tested the default parameters values on a selection of images and reported high accuracy. The range of expected radius is an exception: it greatly affects computational speed and accuracy of the algorithm. We thus suggest to always define it, setting it as narrow as possible. The range boundaries (expected minimal and maximal circle radius in pixels) can be either derived automatically from the smallest and biggest circle drawn on the currently displayed image with the function Define radius range or directly specified within the Options dialog. The following other parameters can be modified: the relative size of imgcrop, minimal distance from imgcrop boundaries, a bias towards foreground classification as well as the minimal proportion classified as foreground within a circle. Maximal proportion of a circle overlapping with one other circle as well as minimal distance and radius difference are modifiable.
Additionally, the user can set maximal total area of a circle overlapping with any other circle and define the final minimal distance between any circle.

Manual corrections
Addition of non-detected colonies and false positive removal is possible through simple mouseguided operations. Additionally, polygons can be drawn on the image to clear out entire zones from found circles. In TL mode, it is important for downstream analysis to check if none of the detected circles is composed from two (or more) colonies which merged during the time-lapse imaging. These should be removed and replaced with the appropriate number of circles with corresponding centers. This control step is needed to be done manually.

Algorithm performance
For the example shown in Fig. 6, a 1956x1936 pixels sized grayscale image displaying 42 Staphylococcus aureus colonies after 70.5 h of growth on Columbia sheep blood with an average radius of 53 pixels, the computation from raw grayscale image to isolated objects took  Table S2. Generally, smaller image size and a narrower range of expected radii improve computational speed but may negatively impact accuracy.
Generally, ColTapp accurately detected colonies when the area of interest, expected radius range and the method for conversion of a color image into a grayscale image were set for each image separately (SI text 2.1). Note, that although colonies are usually lighter than the background, phage plaques are darker than the bacterial lawn: an option to find darker-than-background circles is available.
Overall, false positive rate varied substantially, ranging from 0% to 55.6% (median = 3.54%) ( Supplementary Fig. S5A). Almost all images with false positive rates exceeding 5% had most of their wrongly detected circles in clusters (marked with blue symbols in Supplementary Fig. S5A), usually in areas with lighting artifacts (for example the triangular-shaped light reflection on image 30). These clusters can be cleared with a few clicks and do not pose a problem in our opinion. The one other case with a high false positive rate was an image displaying colonies of the marine bacterial species Alteromonas macleodii and a lot of chitinous debris, forming particles which were only marginally different in size as compared to the bacterial colonies. Therefore, false positive rate should not exceed 5% if ideal imaging conditions are maintained. Additionally, we did not observe any correlation between false positive rate and total number of colonies on a plate ( Supplementary Fig. S5A).
The false negative rate varied between 0% and 21.7% (median = 4.06%) and slightly increased with total number of colonies on the plate ( Supplementary Fig. S5B). Some of the images with high false negative rates had strong lighting gradients, high size heterogeneity or a high number of overlapping colonies all of which can decrease performance of ColTapp's colony detection algorithm.

Morphology descriptors
ColTapp can extract basic colony morphology descriptors (including color, outline and texture) from each frame for downstream analyses, such as species identification or observation of mutation-induced phenotypes linked with altered growth on agar plates. Color values can be extracted from the entire colony or from the center (within a 5 pixel radius) and texture is calculated either as standard deviation of the color or image entropy 3 (Supplementary Figure 4).
The length of the perimeter of colonies as well as the standard deviation on the average radius are metrics useful to categorize colonies based on the shape of the border. We also propose to export the mean color of the halo around colonies as this could be useful for colorimetric assays 4 or quantification of hemolysis capacities of the colony.

Spatial metrics
For post-processing steps such as density correction or the study of colonies interactions, the user can export typical spatial metrics calculations. Interactions between neighboring colonies are mostly occurring through the diffusion of small molecules through the agar plate matrix. The diffusion in such 3D gels can sometimes be approximated to be two-dimensional, and basic solutions of the diffusion equations will give interactions based on the distance D between colonies either as 1 or distance from the focal colony of each neighboring colony. This metric, unexplored in the context of colony growth, is biologically plausible since it assumes that distant large colonies (that have consumed large amount of resources) have a similar effect as small but closer colonies, consuming resources closer to the focal colony. A user may define distance cutoffs to limit the interaction terms to colonies within a maximal distance to the focal colony for further calculations.
Finally, Chacón, et al. 5 proposed that the amount of nutrients available to each colony on the agar plate can be approximated with Voronoi cell areas, obtained by tracing perpendicular bisector lines between each pair of neighboring colonies 6 . They showed that the final size colonies attain when agar plate carrying capacity is reached can be predicted by the colonies Voronoi cell areas.
ColTapp computes these areas by performing a Voronoi tessellation within user-specified boundaries 7 .

Drift correction
Slight drift can occur in time-lapse setups during imaging. Therefore, ColTapp incorporates an image registration algorithm which performs a 2-D rigid translation at subpixels resolution (1/κ), where the factor κ (adjustable as RegistrationFactor in Options) can be specified by the user. The algorithm computes the up-sampled cross-correlation with a discrete Fourier transform 8 . The drift correction is done on a user-defined small area of the image by generating a translation vector in respect to the reference frame. Concretely, the translation is not applied to the images themselves, but to the colonies' centers position.

Colony centers correction
The centering of colony circles is crucial to track radii over time, especially at early timepoints when colonies are small. ColTapp offers an automatic center correction, by tentatively detecting circles on sub-images cropped from an early frame based on the location and radii of the colonies detected on the reference frame. These newly detected circles' center coordinates are kept as colony center, unless they are farther away than a user-defined threshold from the original coordinates or if they are set to the center coordinates of another colony in close proximity. In these cases, the correction is skipped, and the colonies are added to a list for subsequent manual center correction. If no circle is detected at all because the colony is not yet visible at that timepoint, the same process is repeated 10 frames later. The user is able to monitor this process visually.
Alternatively, the user may choose to execute this task manually, by clicking the correct center on the sub-images sequentially displayed.

Overlapping colonies
Neighboring colonies can be visible as blurry regions in the top of kymographs, hindering creation of clean binary images (Fig. 5B). To reduce this phenomenon, before the kymograph creation process, overlapping colonies are automatically detected and the ranges of angles corresponding to adjacent colonies are discarded from the polar transformed intensity data. If more than 90% of all angles are discarded because of overlap, the exclusion of angles is omitted completely, to avoid reducing the available data too much. The overlap detection functionality can be deactivated or tuned with Scale radius for overlap (accessible in Options). This scaling factor is multiplied to the radius of the focal colony from which center neighboring colonies are tested for overlap: by increasing it, a user may choose to discard ranges of angles corresponding not only to overlapping colonies but also very close colonies. Note that this increase might lead to high proportions of angles to be discarded. Decreasing the scaling factor leads to reduced ranges of excluded angles.
This might be useful in densely populated plates to still achieve some overlap exclusion to increase quality of kymographs at earlier timepoints.

Radial growth curve correction
ColTapp has an inbuilt function to automatically detect radial growth curves with probable errors based on the number of local maxima, size of radius differences from frame to frame, monotonicity, and number of frames without a successful radius determination. These (or any) radial growth curves can be manually corrected with a dedicated tool (Fig. S3) which allows the user to switch for each colony individually between Global thresholding and Edge detection method, and adjust any parameter of the two methods to derive best parameter combinations to derive the radial growth curve from the kymograph.

Algorithm performance
Computational

Reference growth parameters definition
The user can define the reference radial growth rate and appearance time either by providing known values manually, or by directing ColTapp to a folder containing a control experiment 7 Sequence of endpoint images

Delayed colony detection
Some very late appearing colonies may not be observed at the optimal time for colony size observation (e.g. 24h). In that case, we suggest capturing image(s) at later timepoint(s) (e.g. 48h).
ColTapp allows the user to overlay images taken from the same plate at different times to ease detection of very late growing colonies which did not yet appear at an earlier timepoint. ColTapp assumes that a folder under analysis contains multiple images from the same timepoint and requires the user to navigate to a folder containing the same set of images with matching order, acquired at another timepoint.
To achieve similar image orientation at both timepoints, the user may place a mark on the edge of the plate as reference. Users can manually align images which have not been captured in the exact same orientation by clicking at two matching positions on the images of the two timepoints.
The inbuilt MATLAB function fitgeotrans is used to fit a geometric transformation to the pairs of points with a nonreflective similarity to estimate rotation and translation of the two images to align. If delayed colonies are revealed at this point, the user has the option to add not yet visible colonies to the earlier timepoint, which are defined as colonies with a radius equal to zero.

Sequence of endpoint images
The user can link corresponding folders of multiple image sets taken at different timepoints to create timeseries data without the need for a time-lapse imaging setup. Colonies on the images in different folders are matched based on the minimal distance of centers between colonies.
Mismatches can be corrected by aligning the images at different timepoints by fitting a geometric transformation to user-defined pairs of points with a nonreflective similarity to estimate rotation and translation of the two images to align. It is necessary that the number of detected colonies is the same in all folders.

Options
Many parameters of implemented functions can be tuned by the user within the Options window ( Fig S2). The Options window is divided in four tabs (Global, Detect, Main (TL or EP) and Visualize).
We describe some of the most important parameters and functions here.
The Global tab allows to set the method for converting a colored image into a grayscale version, activate direct data visualization options on the image (e.g. colony circles and Voronoi edges), set averaging strategy for reference growth data extraction and to define the unit (pixel or micrometer) of the data export function.
The The Main-EP tab contains small functions to assign the AOI and spatial calibration factor of the current frame to all frames of the loaded folder and remove the linked folders.
Finally, the Visualize tab allows to tune details of the data previewing functions. Smoothing factor for radial growth curves, units to use for plotting (pixel or micrometer and frame or hour), number of bins used for distributions and titles can be set here.

Data export
The generated data is automatically saved in the same folder as the images are stored in the native MATLAB format (.mat files) when running the application. For export to another software, ColTapp provides a data export functionality. Through the export menu, user can save comma separated values files of position, radius, fitted appearance time, as well as several calculated density measures and simple morphology descriptors of each colony.
Most variables are exported as a three-column table, containing the frame number, the colony number and the exported data of interest, grouped with similar variables. In TL mode, the radius of colonies is exported as a matrix of time by colony to facilitate extraction of growth curves.
The user can choose to export relevant variables either as pixel values or as µm values, provided the spatial calibration factor was set on all frames. Additional data as for example the spatial calibration factor, used color to grayscale transformation method as well as plate center and radius (if specified) can be exported in a separate metadata csv file. can be brought below 50µm, and a short exponential phase is observed when colonies grow, before they reach ~100 µm. In this case, the appearance time cannot be estimated using a simple linear regression and calculations will need to be adapted.   positive of these and number of missed (false negative). We computed false negative and false positive rates with these numbers. Many of the images with high false positive rate had the false negative circles detected in large clusters, which is indicated on the table. Additionally, we measured computational time of the detection algorithm and report computational time per circle as well. Finally, average size of colonies in pixels and the range of expected radius is shown. here the number of colonies and frames of each time-lapse as well as the number of radial growth curves automatically identified as failed after applying the Global thresholding method, the number of growth curves falsely categorized as failed (including rate), and the number of the identified growth curves which needed more manual correction than a simple switch to the Edge detection method. Corresponding rates. Next, we manually assessed all growth curves not classified as failed and record number of falsely categorized as correct (including false detection rate). Of these manually identified, we assessed if more manual correction than a switch to Edge detection method was necessary. The combined number of automatically and manually identified growth curves is reported here and the rate of failure of the Global thresholding method is given. Overall, only few growth curves need manual correction after switching to the Edge detection method. Additionally, we report the total computational time and the time per colony and frame as these two factors are proportional to required time.

Supplementary
Finally, Average colony size at the last frame of the time-lapse is indicated.

Movies
Supplementary Movie 1: Representative images of the SI movie 1 (one image every 5h, starting at 10h) that show a time-lapse of colony growth. The time-lapse movie is generated from the images of plate 11 (Fig.   S9, Fig. S10). It includes 423 frames acquired in 10 min intervals.

Supplementary Movie 2:
This movie proposes to merge Fig 7BD (full dataset) and show the evolution of the correlation of radius of colonies and their appearance time. The red line represents the radius that could be obtained assuming the maximal growth rate, and deviations from this radius means that colonies grow slower than this maximal growth rate.

Supplementary Movie 3:
This movie proposes to visualize as in Fig. 8A the correlation and fitted model of the systematic error obtained after the use of a linear regression with GRmax to obtain tapp. The 0 line thus represents colonies growing at maximal growth rate, and the fit is described in main text.