Mathematical imaging methods for mitosis analysis in live-cell phase contrast microscopy

Highlights • A workflow for mitosis analysis in phase contrast microscopy is proposed.• It is derived from mathematical imaging concepts.• Those include the circular Hough transform and variational segmentation methods.• Results for live-cell imaging experiments are presented.


Introduction
Mathematical image analysis techniques have recently become enormously important in biomedical research, which increasingly needs to rely on information obtained from images. Applications range from sparse sampling methods to enhance image acquisition through structure-preserving image reconstruction to automated analysis for objective interpretation of the data [1]. In cancer research, observation of cell cultures in live-cell imaging experiments by means of sophisticated light microscopy is a key technique for quality assessment of anti-cancer drugs [2,3]. In this context, analysis of the mitotic phase plays a crucial role. The balance between mitosis and apoptosis is normally carefully regulated, but many types of cancerous cells have evolved to allow uncontrolled cell division. Hence drugs targeting mitosis are used extensively during cancer chemotherapy. In order to evaluate the effects of a given drug on mitosis, it is desirable to measure average mitosis durations and distribution of possible outcomes such as regular division into two daughter cells, apoptosis, division into an abnormal number of daughter cells (one or P 3) and no division at all [4,5].
Since performance of technical equipment such as microscopes and associated hardware is constantly improving and large amounts of data can be acquired in very short periods of time, automated image processing tools are frequently favoured over manual analysis, which is expensive and prone to error and bias. Generally, experiments might last several days and images are taken in a magnitude of minutes and from different positions. This leads to a sampling frequency of hundreds of images per sequence with an approximate size of 1000 2 pixels.

Image characteristics in phase contrast microscopy
In live-cell imaging experiments for anti-cancer drug assessment, the imaging modality plays a key role. Observation of cell cultures originating from specific cell lines under the microscope requires a particular setting ensuring that the cells do not die during image acquisition and that they behave as naturally as possible [6]. Here, phase contrast is often preferred to fluorescence microscopy because the latter requires labelling or transgenic expression of fluorescent markers, both causing phototoxicity and possibly changes of cell behaviour [7][8][9]. As opposed to this, cells do not need to be stained for phase contrast microscopy. Moreover, phase shifts facilitate visualisation of even transparent specimens as opposed to highlighting of individual specific cellular components in fluorescence microscopy. We believe that one main advantage of our proposed framework is that it can be applied to data acquired with any standard phase-contrast microscope, which are prevalent in many laboratories and more widespread than for instance recently established quantitative phase imaging devices (e.g. Q-Phase by Tescan).
There are two common image characteristics occurring in phase contrast imaging (cf. Fig. 1). Both visual effects highly impede image processing and standard algorithms are not applicable in a straightforward manner. The shade-off effect leads to similar intensities inside the cells and in the background. As a result, edges are only weakly pronounced and imaging methods such as segmentation relying on intensity gradient information (cf. Section 2.2.2) often fail. Moreover, region-based methods assuming that average intensities of object and background differ from one another (cf. Section 2.2.3) are not applicable either. Secondly, the halo effect is characterised by areas of high intensity surrounding cell membranes. The brightness levels increase significantly immediately before cells enter mitosis due to the fact that they round up, form a nearly spherically-shaped volume and therefore the amount of diffracted light increases. In addition, both effects prohibit application of basic image pre-processing tools like for example thresholding or histogram equalisation (cf. [10]).

Brief literature review
Over the past few years a lot of cell tracking frameworks have been established (cf. [11]) and some publications also feature mitosis detection. In [12], a two-step cell tracking algorithm for phase contrast images is presented, where the second step involves a level-set-based variational method. However, analysis of the mitotic phase is not included in this framework. Another tracking method based on extended mean-shift processes [13] is able to incorporate cell divisions, but does not provide cell membrane segmentation. In [14] an automated mitosis detection algorithm based on a probabilistic model is presented, but it is not linked to cell tracking. A combined mitosis detection and tracking framework is established in [15], although cell outline segmentation is not included. Li et al. [16] provide a comprehensive framework facilitating both tracking and lineage reconstruction of cells in phase contrast image sequences. Moreover, they are able to distinguish between mitotic and apoptotic events.
In addition, a number of commercial software packages for semi-or fully automated analysis of microscopy images exist, for example Volocity, Columbus (both PerkinElmer), Imaris (Bitplane), ImageJ/Fiji [17] and Icy [18] (also cf. [19]). The last two are open source platforms and the latter supports graphical protocols while the former incorporates a macro language, allowing for individualisation and extension of integrated tools. However, the majority of plugins and software packages are limited to analysis of fluorescence data.
A framework, which significantly influenced development of our methods and served as a basis for our tracking algorithm, was published in 2014 by Möller et al. [20]. It incorporates a MATLAB Ò Graphical User Interface that enables semi-automated tracking of cells in phase contrast microscopy time-series. The user has to manually segment the cells of interest in the first frame of the image sequence and can subsequently execute an automatic tracking procedure consisting of two rough and refined segmentation steps. In the following section, the required theoretical foundations of mathematical imaging methods are discussed, starting with the concept of the circular Hough transform and continuing with a review of segmentation and tracking methods leading to a more detailed description of the above-mentioned framework. For a more detailed discussion, we refer the interested reader to [10] and the references therein.

The circular Hough transform
The Hough transform is a method for automated straight line recognition in images patented by Paul Hough in 1962 [21]. It was further developed and generalised by Duda and Hart in 1972 [22]. More specifically, they extended the Hough transform to different types of parametrised curves and in particular, they applied it to circle detection.
The common strategy is to transform points lying on straight line segments or curves in the underlying image into a parameter space. Its dimension depends on the number of variables required in order to parametrise the sought-after curve. For the parametric representation of a circle, which can be written as the radius r as well as two centre coordinates ðc 1 ; c 2 Þ are required. Hence, the corresponding parameter space is threedimensional. Each point ðx; yÞ in the original image satisfying the above equation for fixed r; c 1 and c 2 coincides with a cone in the parameter space. Then, edge points of circular objects in the original image correspond to intersecting cones and from detecting those intersections in the parameter space one can again gather circles in the image space.
For simplification, we fix the radius and consider the twodimensional case in Fig. 2. On the left, we have the image space, i.e. the x-y-plane, and a circle in light blue with five arbitrary points located on its edge highlighted in dark blue. All points fulfil Eq. (1) for fixed centre coordinates ðc 1 ; c 2 Þ. On the other hand, fixing those specific values for c 1 and c 2 in the parameter space, i.e. c 1 -c 2 -plane, on the right, and keeping x and y in (1) arbitrary, leads to the dashed orange circles, where the corresponding edge points are drawn in grey for orientation. All of the orange circles intersect in one point, which exactly corresponds to the circle centre in the original image. Hence, from intersections in the parameter space one can reference back to circular objects in the image space.
A discussion on how the circular Hough transform is embedded and implemented in MitosisAnalyser can be found in Section 3.1.

Image segmentation and tracking
In the following, we would like to introduce variational methods (cf. e.g. [23,24]) for imaging problems. The main aim is minimi- sation of an energy functional modelling certain assumptions on the given data and being defined as It is dependent on the solution /, which represents the processed image to be obtained, and shall be minimised with respect to /. The given image to be processed is denoted by w. The functions / and w map from the rectangular image domain ues. In the case of 8-bit phase contrast microscopy images, d ¼ 1 and R ¼ f0; . . . ; 255g, where 0 and 255 correspond to black and white, respectively.
The first part D on the right-hand side of (2) ensures data fidelity between / and w, i.e. the solution / should be reasonably close to the original input data w. This can be obtained by minimising a norm measuring the distance between w and /, where the choice of norm naturally depends on the given problem. The regulariser R in (2) incorporates a priori knowledge about the function /. For example, / could be constrained to be sufficiently smooth in a particular sense. The parameter a is weighting the two different terms and thereby defines which one is considered to be more important. Energy functionals can also consist of multiple data terms and regularisers. Eventually, a solution that minimises the energy functional (2) attains a small value of D assuring high fidelity to the original data, of course depending on the weighting. Similarly, a solution which has a small value of R can be interpreted as having a high coincidence with the incorporated prior assumptions.
Here, we focus on image segmentation. The goal is to divide a given image into associated parts, e.g. object(s) and background. This can be done by finding either the objects themselves or the corresponding edges, which is then respectively called regionbased and edge-based segmentation. However, those two tasks are very closely related and even coincide in the majority of cases. Tracking can be viewed as an extension of image segmentation because it describes the process of segmenting a sequence of images or video. The goal of object or edge identification remains the same, but the time-dependence is an additional challenge.
Below, we briefly discuss the level-set method and afterwards present two well-established segmentation models incorporating the former. Furthermore, we recap the methods in [20] building upon the above and laying the foundations for our proposed tracking framework.

The level-set method
In 1988 the level-set method was introduced by Osher and Sethian [25]. The key idea is to describe motion of a front by means of a time-dependent partial differential equation. In variational segmentation methods, energy minimisation corresponds to propagation of such a front towards object boundaries. In two dimensions, a segmentation curve c is modelled as the zero-level of a three-dimensional level-set function /. Two benefits are straightforward numerical implementation without need of parametrisation and implicit modelling of topological changes of the curve. The level-set evolution equation can be written as @/ @t ¼ FÁ j r/ j with curvature-dependent speed of movement F and suitable initial and boundary conditions.
For implementation, the level-set function / is assigned negative values inside and positive values outside of the curve c, commonly chosen to be the signed Euclidean distances (cf. Fig. 3).

Geodesic active contours
Active contours or ''snakes" have been developed and extended for decades [26][27][28][29][30] and belong to the class of edge-based segmentation methods. As the name suggests, the goal is to move segmentation contours towards image edges and stop at boundaries of objects to be segmented (e.g. by using the level-set method described above). Geodesic active contours constitute a specific type of active contours methods and have been introduced by Caselles, Kimmel and Sapiro in 1997 [31]. The level-set formulation reads @/ @t with appropriate initial and boundary conditions and g is an edge-detector function typically depending on the gradient magnitude of a smoothed version of a given image w. A frequently used function is with G r being a Gaussian kernel with standard deviation r. The function g is close to zero at edges, where the gradient magnitude is high, and close or equal to one in homogeneous image regions, where the gradient magnitude is nearly or equal to zero. Hence, the segmentation curve, i.e. the zero-level of /, propagates towards edges defined by g and once the edges are reached, evolution is stopped. In the specific case of g ¼ 1, (4) coincides with mean curvature motion. Geodesic active contours are a well-suited method of choice for segmentation if image edges are strongly pronounced or can otherwise be appropriately identified by a suitable function g.

Active contours without edges
As the name suggests, the renowned model developed by Chan and Vese [32] is a region-based segmentation method and in contrast to the model presented in 2.2.2, edge information is not taken into account. It is rather based on the assumption that the underlying image can be partitioned into two regions of approximately piecewise-constant intensities. In the level-set formulation the variational energy functional reads which is to be minimised with respect to / as well as c 1 and c 2 . Recalling (3), we define the Heaviside function H as indicating the sign of the level-set function and therefore the position relative to the segmentation curve.
In (6) the structure in (2) is resembled. The first two data terms enforce a partition into two regions with intensities c 1 inside and c 2 outside of the segmentation contour described by the zerolevel-set. The third and fourth terms are contour length and area regularisers, respectively.
The optimal c 1 and c 2 can be directly calculated while keeping / fixed: In order to find the optimal / and hence the sought-after segmentation contour, the Euler-Lagrange equation defined as @/ @t ¼ À @E @/ ¼ 0 needs to be calculated, which leads to the evolution equation where d e is the following regularised version of the Dirac delta function: Eq. (8) can be numerically solved with a gradient descent method.
This model is very advantageous for segmenting noisy images with weakly pronounced or blurry edges as well as objects and clustering structures of different intensities in comparison to the background.

Tracking framework by Möller et al.
The cell tracking framework developed in [20] is sub-divided into two steps. First, a rough segmentation based on the model in Section 2.2.3 is performed. The associated energy functional reads In contrast to (6), the area or volume regularisation term weighted by m is altered such that the current volume shall be close to the previous volume V old . Moreover, the data terms weighted by k 1 and k 2 incorporate the normal velocity image j v j instead of the image intensity function w: where the expression in the denominator is a regularisation of the gradient magnitude defined as for small e. The novelty here is that in contrast to only considering the image intensity both spatial and temporal information is used in order to perform the regionbased segmentation. Indeed, cells are expected to move between subsequent frames. In addition, the gradient magnitude shall be increased in comparison to background regions. Therefore the incorporation of both temporal and spatial derivative provides a better indicator of cellular interiors. In a second step, a refinement is performed using the geodesic active contours Eq. (4). The edge-detector function is customised and mainly uses information obtained by the Laplacian of Gaussian of the underlying image. In addition, topology is preserved throughout the segmentation by using the simple points scheme [33][34][35] and in order to reduce computational costs this is combined with a narrow band method [36], which we inherit in our framework as well.

MitosisAnalyser framework
In the following we present our proposed workflow designed in order to facilitate mitosis analysis in live-cell phase contrast imaging experiments. We specifically focused on applicability and usability while providing a comprehensive tool that needs minimal user interaction and parameter tuning. The MATLAB Ò Graphical User Interface MitosisAnalyser (The corresponding code is available at github.com/JoanaGrah/MitosisAnalyser.) provides a userfriendly application, which involves sets of pre-determined parameters for different cell lines and has been designed for non-experts in mathematical imaging.
In Fig. 4 the main application window is displayed on the top left. The entire image sequence at hand can be inspected and after analysis, contours are overlaid for immediate visualisation. Moreover, images can be examined and pre-processed by means of a few basic tools (centre), although the latter did not turn out to be necessary for our types of data. Parameters for both mitosis detection and tracking can be reviewed, adapted and permanently saved for different cell lines in another separate window (bottom left). Mitosis detection can be run separately and produces intermediate results, where all detected cells can be reviewed and parameters can be adjusted as required. Consecutively, running the cell tracking algorithm results in an estimate of average mitosis duration and provides the possibility to survey further statistics (right). Fig. 5 summarises the entire workflow from image acquisition to evaluation of results. First, live-cell imaging experiments are conducted using light microscopy resulting in 2D greyscale image sequences. Next, mitosis detection is performed. For each detected cell, steps 3-5 are repeated. Starting at the point in time where the cell is most circular, the circle-shaped contour serves as an initialisation for the segmentation. The tracking is then performed backwards in time, using slightly extended contours from previous frames as initialisations. As soon as cell morphology changes, i.e. area increases and circularity decreases below a predetermined threshold, the algorithm stops and marks the point in time at hand as start of mitosis. Subsequently, again starting from the detected mitotic cell, tracking is identically performed forwards in time until the cell fate can be determined. As already mentioned in Section 1, different cases need to be distinguished from one another: regular, abnormal and no division as well as apoptosis. The final step comprises derivation of statistics on mitosis duration and cell fate distribution as well as evaluation and interpretation thereof.
The double arrow connecting steps 1 and 5 indicates what is intended to be subject of future research. Ideally, image analysis shall be performed in on-line time during image acquisition and intermediate results shall be passed on to inform and influence microscopy software. Consequently, this may in turn lead to enhancement of image processing. Recently established concepts of bilevel optimisation and parameter learning for variational imaging models (cf. [37,38]) might supplement our framework.

Mitosis detection
In order to implement the circular Hough transform (CHT) described in Section 2.1, both image and parameter space need to be discretised. The former is naturally already represented as a pixel grid or matrix of grey values. The latter needs to be artificially discretised by binning values for r; c 1 and c 2 and the resulting representation is called accumulator array. Once the CHT is performed for all image pixels, the goal is to find peaks in the accumulator array referring to circular objects.
There are several options in order to speed up the algorithm, but we will only briefly discuss two of them. First, it is common to perform edge detection on the image before applying the CHT, since pixels lying on a circle very likely correspond to edge pixels. An edge map can for instance be calculated by thresholding the gradient magnitude image in order to obtain a binary image. Then, only edge pixels are considered in the following steps. Furthermore, it is possible to reduce the accumulator array to two dimensions using the so-called phase-coding method. The idea is using complex values in the accumulator array with the radius information encoded in the phase of the array entries. Both enhancements are included in the built-in MATLAB Ò function imfindcircles.
The mitosis detection algorithm implemented into MitosisAnalyser uses this function in order to perform the CHT and search for circular objects in the given image sequences. Fig. 6 visualises the different steps from calculation of the gradient image, to identification of edge pixels, to computation of the accumulator matrix and transformation thereof by filtering and thresholding, to detection of maxima.
This method turned out to be very robust and two main advantages are that circles of different sizes can be found and even not perfectly circularly shaped or overlapping objects can be detected. At the beginning of analysis, the CHT is applied in every image of the given image sequence in order to detect nearly circularly shaped mitotic cells. Afterwards, the circles are sorted by significance, which is related to the value of the detected peak in the corresponding accumulator array. The most significant ones are picked while simultaneously ensuring that identical cells are neither detected multiple times in the same frame nor in consecutive frames. The complete procedure is outlined in Supplementary Algorithm 1.

Cell tracking
We have already introduced variational segmentation methods in general as well as three models our framework is based on in more detail in Section 2.2. Here, we would like to state the cell tracking model we developed starting from the one presented in Section 2.2.4. The energy functional reads: with j v j and H defined as in (10) and (7), respectively.
The two terms weighted by k 1 and k 2 are identical to the ones in (9). Instead of having two separate segmentation steps as in [20], we integrate the edge-based term weighted by m into our energy functional. However, using a common edge-detector function based on the image gradient like the one in (5) was not suitable for our purposes. We noticed that the gradient magnitude image contains rather weakly pronounced image edges, which motivated us to search for a better indicator of the cells' interiors. We realised that the cells are very inhomogeneous in contrast to the background and consequently, we decided to base the edge-detector function on the local standard deviation of grey values in a 3Â3neighbourhood around each pixel. Additionally smoothing the underlying image with a standard Gaussian filter and rescaling intensity values leads to an edge-detector function, which is able to indicate main edges and attract the segmentation contour towards them.
Furthermore, we add a standard length regularisation term weighted by l. We complement our energy functional with an area regularisation term that incorporates a priori information about the approximate cell area and prevents contours from becoming too small or too large. This penalty method facilitates incorporation of a constraint in the energy functional and in this case the area shall not fall below the threshold t area .   Optimal parameters c 1 and c 2 can be calculated directly. We numerically minimise (11) with respect to the level-set function / by using a gradient descent method (cf. 2.2.3). The third term weighted by l is discretised using a combination of forwards, backwards and central finite differences as proposed in [32]. We obtain the most stable numerical results by applying central finite differences to all operators contained in the fourth term weighted by m. In Fig. 7 we visualise level-set evolution throughout the optimisation procedure.
In order to give an overview of the backwards and forwards tracking algorithms incorporated in the mitosis analysis framework, we state the procedures in Supplementary Algorithm 2 and 3. Together with the mitosis detection step they form the foundation of the routines included in MitosisAnalyser.

Material and methods
The MitosisAnalyser framework is tested in three experimental settings with MIA PaCa-2 cells, HeLa Aur A cells and T24 cells. Below, a description of cell lines and chemicals is followed by details on image acquisition and standard pre-processing.

Cell lines and chemicals
The FUCCI (Fluorescent Ubiquitination-based Cell Cycle Indicator [39])-expressing MIA PaCa-2 cell line was generated using the FastFUCCI reporter system and has previously been characterised and described [40,41]. Cells were cultured in phenol red-free Dulbecco's modified Eagle's medium (DMEM) supplemented with 10% foetal calf serum (FBS).
HeLa Aur A cells, HeLa cells modified to over-express aurora kinase A, were generated by Dr. Jennifer Harrington with Dr. David Perera at the Medical Research Council Cancer Unit, Cambridge, using the Flp-In T-REx system from Invitrogen as described before [42]. The parental HeLa LacZeo/TO line, and pOG44 and pcDNA5/ FRT/TO plasmids were kindly provided by Professor Stephen Taylor, University of Manchester. The parental line grows under selection with 50 lg/ml Zeocin TM (InvivoGen) and 4 lg/ml Blasticidin (Invitrogen). HeLa Aur A cells were cultured in DMEM supplemented with 10% FBS and 4 lg/ml blasticidin (Invitrogen) and 200 lg/ml hygromycin (Sigma Aldrich). Transgene expression was achieved by treatment with 1 lg/ml doxycycline (Sigma Aldrich).
In all experiments, all cells were grown at 37°C and 5% CO 2 up to a maximum of 20 passages and for fewer than 6 months following resuscitation. They were also verified to be mycoplasma-free using the Mycoprobe Ò Mycoplasma Detection Kit (R&D Systems). Paclitaxel (Tocris Bioscience), MLN8237 (Stratech Scientific) and Docetaxel (Sigma Aldrich) were dissolved in dimethylsulphoxide (DMSO, Sigma) in aliquots of 30 mM, kept at À20°C and used within 3 months. Final DMSO concentrations were kept constant in each experiment ð6 0:2%Þ.

Acquisition and processing of live-cell time-lapse sequences
Cells were seeded in l-Slide glass bottom dish (ibidi) and were kept in a humidified chamber under cell culture conditions (37°C, 5% CO 2 ). For experiments with T24 and HeLa Aur A cells they were cultured for 24 h before being treated with drugs or DMSO control. They were then imaged for up to 72 h. Images were taken from three to five fields of view per condition, every 5 min, using a Nikon Eclipse TE2000-E microscope with a 20X (NA 0.45) long-working distance air objective, equipped with a sCMOS Andor Neo camera acquiring 2048 Â 2048 images, which have been binned by a factor of two. Red and green fluorescence of the FUCCI-expressing cells were captured using a pE-300white CoolLED source of light filtered by Nikon FITC B-2E/C and TRITC G-2E/C filter cubes, respectively. For processing, an equalisation of intensities over time was applied to each channel, followed by a shading correction and a background subtraction, using the NIS-Elements software (Nikon).

Results and discussion
In this section we present and discuss results obtained by applying MitosisAnalyser to the aforementioned experimental live-cell imaging data. A list of parameters we chose can be found in Supplementary Table 1. For each cell line, we established a unique set of parameters. Nevertheless, the individual values are in reasonable ranges and do not differ significantly from one another. We did not follow a specific parameter choice rule, but rather tested various combinations and manually picked the best performing ones.

MIA PaCa-2 cells
In a multi-modal experiment with FUCCI-expressing MIA PaCa-2 cells, both phase contrast images and fluorescence data were acquired. The latter consist of two channels with red and green intensities corresponding to CDT1 and Geminin signals, respectively. In this case we do use fluorescence microscopy imaging data as well, but we would like to stress that this analysis would not have been possible without the mitosis detection and tracking performed on the phase contrast data. As before, mitotic cells are detected using the circular Hough transform applied to the phase contrast images. Cell tracking is performed on the phase contrast images as well, but in addition, information provided by the green fluorescent data channel is used. More specifically, stopping criteria for both backwards and forwards tracking are based on green fluorescent intensity distributions indicating different stages of the cell cycle, which can be observed and is described in more detail in Supplementary Fig. 1.
The whole data set consists of nine imaging positions, where three at a time correspond to DMSO control, treatment with 3nM paclitaxel and treatment with 30nM paclitaxel. Fig. 8 visualises exemplary courses of the mitotic phase, which could be measured by means of our proposed workflow. Table 1 presents estimated average mitosis durations for the three different classes of data. Indeed, the average duration of 51 min for the control is consistent with that obtained from manual scoring (cf. [41], Figure S3D). Moreover, we can observe a dose-dependent increase in mitotic duration for the two treatments, which was anticipated, since paclitaxel leads to mitotic arrest.

HeLa cells
In the following we discuss results achieved by applying Mito-sisAnalyser to sequences of phase contrast microscopy images showing HeLa Aur A cells. In addition to DMSO control data, cells have been treated with 25 nM MLN8237 (MLN), 0.75 nM paclitaxel (P), 30 nM paclitaxel (P) and with a combination of 25 nM MLN8237 and 0.75 nM paclitaxel (combined). Fig. 9 shows exemplary results for detected and tracked mitotic events, where DMSO control cells divide regularly into two daughter cells. Particular treatments are expected to enhance multipolar mitosis and indeed our framework was able to depict the three daughter cells in each of the three examples (bottom rows) presented. In addition, mitosis duration is extended, as anticipated, for treated cells and specifically for the combined treatment. The segmentation of the cell membranes seems to work well by visual inspection, even in the case of touching neighbouring cells. Table 2 summarises average mitosis durations that have been estimated for the different treatments. Again, the results are according to our expectations, i.e. mitosis durations for treated cells are extended in comparison to DMSO control.

T24 cells
For this data set we wanted to focus on cell fate determination and in order to distinguish between different fates in the T24 cell data set we combine the MitosisAnalyser framework with basic classification techniques. In particular, we manually segmented three different classes of cells: mitotic and apoptotic ones as well as cells in their normal state outside of the mitotic cell cycle phase (see Fig. 10).
In Fig. 11 we show boxplots of nine features based on morphology as well as intensity values we use for classification. Those include area, perimeter and circularity. Furthermore, we calculate both mean and standard deviation of the histogram. In addition, we consider the maximum of the gradient magnitude, the mean as well as the total variation of the local standard deviation and the total variation of the grey values. One can clearly observe that cells in mitosis have much higher circularity than in any other state. Flat cells differ significantly from the other two classes with respect to features based on intensity values.
In order to train a classifier solely based on those few features we used the MATLAB Ò Machine Learning Toolbox and its accompanying Classification Learner App. We chose a nearest-neighbour classifier with the number of neighbours set to 1 using Euclidean distances and equal distance weights, which yielded a classification accuracy of 93.3% (cf. Supplementary Fig. 2).
Pie charts for T24 cell fate distributions for different drug treatments as preliminary results can be found in Supplementary Fig. 3, although integration of classification techniques will be subject of more extensive future research.

Validation
In order to validate performance of the segmentation, we compare results obtained with MitosisAnalyser with blind manual segmentation. For that purpose, we choose two different error measures: The Jaccard Similarity Coefficient (JSC) [43] and the Modified Hausdorff Distance (MHD) [44], which we are going to define in the following.
Let A and M be the sets of pixels included in the automated and manual segmentation mask, respectively. The JSC is defined as where where dða; MÞ ¼ min m2M ka À mk with Euclidean distance k Á k.
The JSC assumes values between 0 and 1 and the closer it is to 1 the better is the segmentation quality. The MHD on the other hand is equal to 0 if two shapes coincide and the larger the number, the farther they differ from each other. In Fig. 12 and Supplementary Table 2 we can observe that on average, MitosisAnalyser performs better than the standard Chan-Vese method (cf. Section 2.2.3) and Geodesic Active Contours based on the gradient magnitude (cf. Section 2.2.2) (both performed using the MATLAB imageSegmenter application) compared to manual segmentation of ten apoptotic T24 cell images (cf. Fig. 10,, top row). Moreover, Fig. 13 shows successful segmentation of flat T24 cells affected by the shade-off effect in phase contrast microscopy images using Mito- Fig. 8. Three examples of mitotic events detected for FUCCI MIA PaCa-2 ''DMSO control", ''treatment with 3 nM paclitaxel" and ''treatment with 30 nM paclitaxel" data (from top to bottom).  sisAnalyser, where both the method by Chan and Vese and geodesic active contours failed.

Conclusions
We have used concepts of mathematical imaging including the circular Hough transform and variational tracking methods in order to develop a framework that aims at detecting mitotic events and segmenting cells in phase contrast microscopy images, whilst overcoming the difficulties associated with those images. Originating from the models presented in Section 2, we developed a customised workflow for mitosis analysis in live-cell imaging experiments performed in cancer research and discussed results we obtained by applying our methods to different cell line data.