# Mathematical Imaging Methods for Mitosis Analysis in Live-Cell Phase Contrast Microscopy

###### Abstract

In this paper we propose a workflow to detect and track mitotic cells in time-lapse microscopy image sequences. In order to avoid the requirement for cell lines expressing fluorescent markers and the associated phototoxicity, phase contrast microscopy is often preferred over fluorescence microscopy in live-cell imaging. However, common specific image characteristics complicate image processing and impede use of standard methods. Nevertheless, automated analysis is desirable due to manual analysis being subjective, biased and extremely time-consuming for large data sets. Here, we present the following workflow based on mathematical imaging methods. In the first step, mitosis detection is performed by means of the circular Hough transform. The obtained circular contour subsequently serves as an initialisation for the tracking algorithm based on variational methods. It is sub-divided into two parts: in order to determine the beginning of the whole mitosis cycle, a backwards tracking procedure is performed. After that, the cell is tracked forwards in time until the end of mitosis. As a result, the average of mitosis duration and ratios of different cell fates (cell death, no division, division into two or more daughter cells) can be measured and statistics on cell morphologies can be obtained. All of the tools are featured in the user-friendly MATLAB Graphical User Interface MitosisAnalyser.

Keywords: Phase Contrast Microscopy, Mitosis Analysis, Circular Hough Transform, Cell Tracking, Variational Methods, Level-Set Methods

## 1 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 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 pixels.

### 1.1 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. Figure 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]).

### 1.2 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.

## 2 Mathematical Background

### 2.1 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

(1) |

the radius as well as two centre coordinates are required. Hence, the corresponding parameter space is three-dimensional. Each point in the original image satisfying the above equation for fixed , and 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 two-dimensional case in Figure 2. On the left, we have the image space, i.e. the --plane, and a circle in light blue with five arbitrary points located on its edge highlighted in dark blue. All points fulfil equation (1) for fixed centre coordinates . On the other hand, fixing those specific values for and in the parameter space, i.e. --plane, on the right, and keeping and 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.

### 2.2 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 minimisation of an energy functional modelling certain assumptions on the given data and being defined as

(2) |

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 . The functions and map from the rectangular image domain to containing colour () or greyscale () intensity values. In the case of 8-bit phase contrast microscopy images, and , where 0 and 255 correspond to black and white, respectively.

The first part on the right-hand side of (2) ensures data fidelity between and , i.e. the solution should be reasonably close to the original input data . This can be obtained by minimising a norm measuring the distance between and , where the choice of norm naturally depends on the given problem. The regulariser in (2) incorporates a-priori knowledge about the function . For example, could be constrained to be sufficiently smooth in a particular sense. The parameter 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 assuring high fidelity to the original data, of course depending on the weighting. Similarly, a solution which has a small value of 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 region-based 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.

#### 2.2.1 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 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

with curvature-dependent speed of movement and suitable initial and boundary conditions.

For implementation, the level-set function is assigned negative values inside and positive values outside of the curve ,

(3) |

commonly chosen to be the signed Euclidean distances (cf. Figure 3).

#### 2.2.2 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

(4) |

with appropriate initial and boundary conditions and is an edge-detector function typically depending on the gradient magnitude of a smoothed version of a given image . A frequently used function is

(5) |

with being a Gaussian kernel with standard deviation . The function 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 and once the edges are reached, evolution is stopped. In the specific case of , (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 .

#### 2.2.3 Active Contours without Edges

As the name suggests, the renowned model developed by Chan and Vese in 2001 [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

(6) |

which is to be minimised with respect to as well as and . Recalling (3), we define the Heaviside function as

(7) |

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 inside and outside of the segmentation contour described by the zero-level-set. The third and fourth terms are contour length and area regularisers, respectively.

The optimal and 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 needs to be calculated, which leads to the evolution equation

(8) |

where is the following regularised version of the Dirac delta function:

Equation (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.

#### 2.2.4 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

(9) |

In contrast to (6), the area or volume regularisation term weighted by is altered such that the current volume shall be close to the previous volume . Moreover, the data terms weighted by and incorporate the normal velocity image instead of the image intensity function :

(10) |

where the expression in the denominator is a regularisation of the gradient magnitude defined as for small . 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 region-based 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 equation (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.

## 3 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 user-friendly application, which involves sets of pre-determined parameters for different cell lines and has been designed for non-experts in mathematical imaging.

In Figure 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).

Figure 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.

### 3.1 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 , and 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. Figure 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 A1.

### 3.2 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:

(11) |

The two terms weighted by and 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 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 33-neighbourhood 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 . 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 .

Optimal parameters and 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 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 . In Figure 7 we visualise level-set evolution throughout the optimisation procedure.

## 4 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.

### 4.1 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).

T24 cells were acquired from CLS. The T24 cells were cultured in DMEM/F12 (1:1) medium supplemented with 5% 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 g/ml Zeocin™(InvivoGen) and 4 g/ml Blasticidin (Invitrogen). HeLa Aur A cells were cultured in DMEM supplemented with 10% FBS and 4 g/ml blasticidin (Invitrogen) and 200 g/ml hygromycin (Sigma Aldrich). Transgene expression was achieved by treatment with 1 g/ml doxycycline (Sigma Aldrich).

In all experiments, all cells were grown at 37C and 5% CO 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 30mM, kept at -20C and used within 3 months. Final DMSO concentrations were kept constant in each experiment (0.2%).

### 4.2 Acquisition and Processing of Live-Cell Time-Lapse Sequences

Cells were seeded in -Slide glass bottom dish (ibidi) and were kept in a humidified chamber under cell culture conditions (37C, 5% CO). For experiments with T24 and HeLa Aur A cells they were cultured for 24 hours before being treated with drugs or DMSO control. They were then imaged for up to 72 hours. Images were taken from three to five fields of view per condition, every 5 minutes, 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 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).

## 5 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 A1. 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.

### 5.1 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 Figure A1.

DMSO control | 3nM paclitaxel | 30nM paclitaxel | |||||||

Pos 1 | Pos 2 | Pos 3 | Pos 4 | Pos 5 | Pos 6 | Pos 7 | Pos 8 | Pos 9 | |

Events | 14 | 11 | 13 | 12 | 8 | 19 | 10 | 13 | 35 |

AMD | 51 | 41 | 60 | 52 | 88 | 94 | 146 | 104 | 112 |

Total AMD | 51 | 78 | 121 |

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. Figure 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 minutes 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.

### 5.2 HeLa Cells

In the following we discuss results achieved by applying MitosisAnalyser to sequences of phase contrast microscopy images showing HeLa Aur A cells. In addition to DMSO control data, cells have been treated with 25nM MLN8237 (MLN), 0.75nM paclitaxel (P), 30nM paclitaxel (P) and with a combination of 25nM MLN8237 and 0.75nM paclitaxel (combined).

Figure 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.

DMSO | 25nM MLN | 0.75nM P | 30nM P | Combined | |
---|---|---|---|---|---|

Events | 44 | 75 | 10 | 35 | 43 |

AMD | 58 | 73 | 68 | 116 | 105 |

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.

### 5.3 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 Figure 10).

In Figure 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 Figure A2).

Pie charts for T24 cell fate distributions for different drug treatments as preliminary results can be found in Supplementary Figure A3, although integration of classification techniques will be subject of more extensive future research.

### 5.4 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 denotes the intersection of sets A and M, which contains pixels that are elements of both A and M. The union of sets A and M, denoted by , contains pixels that are elements of A or M, i.e. elements either only of A or only of M or of . The MHD is a generalisation of the Hausdorff distance, which is commonly used to measure distance between shapes. It is defined as

where with Euclidean distance .

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 Figure 12 and Supplementary Table A2 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. Figure 10, top row). Moreover, Figure 13 shows successful segmentation of flat T24 cells affected by the shade-off effect in phase contrast microscopy images using MitosisAnalyser, where both the method by Chan and Vese and geodesic active contours failed.

### 5.5 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.

## 6 Acknowledgements

JSG acknowledges support by the NIHR Cambridge Biomedical Research Centre and would like to thank Hendrik Dirks, Fjedor Gaede [45] and Jonas Geiping [46] for fruitful discussions in the context of a practical course at WWU Münster in 2014 and significant speed-up and GPU implementation of earlier versions of the code. JSG and MB would like to thank Michael Möller for providing the basic tracking code and acknowledge support by ERC via Grant EU FP 7 - ERC Consolidator Grant 615216 LifeInverse. MB acknowledges further support by the German Science Foundation DFG via Cells-in-Motion Cluster of Excellence. CBS acknowledges support from the EPSRC grant Nr. EP/M00483X/1, from the Leverhulme grant “Breaking the non-convexity barrier”, from the EPSRC Centre for Mathematical And Statistical Analysis Of Multimodal Clinical Imaging grant Nr. EP/N014588/1, and the Cantab Capital Institute for the Mathematics of Information. JAH, SBK, JAP, AS and SR were funded by Cancer Research UK, The University of Cambridge and Hutchison Whampoa Ltd. SBK also received funding from Pancreatic Cancer UK.

## References

- [1] Jens Rittscher. Characterization of biological processes through automated image analysis. Annual Review of Biomedical Engineering, 12:315–344, 2010.
- [2] Caroline H Topham and Stephen S Taylor. Mitosis and apoptosis: how is the balance set? Current opinion in cell biology, 25(6):780–785, 2013.
- [3] Karen E Gascoigne and Stephen S Taylor. Cancer cells display profound intra-and interline variation following prolonged exposure to antimitotic drugs. Cancer cell, 14(2):111–122, 2008.
- [4] Conly L Rieder and Helder Maiato. Stuck in division or passing through: what happens when cells cannot satisfy the spindle assembly checkpoint. Developmental cell, 7(5):637–651, 2004.
- [5] Beth AA Weaver and Don W Cleveland. Decoding the links between mitosis, cancer, and chemotherapy: The mitotic checkpoint, adaptation, and cell death. Cancer cell, 8(1):7–12, 2005.
- [6] David J Stephens and Victoria J Allan. Light microscopy techniques for live cell imaging. Science, 300(5616):82–86, 2003.
- [7] Ram Dixit and Richard Cyr. Cell damage and reactive oxygen species production induced by fluorescence microscopy: effect on mitosis and guidelines for non-invasive fluorescence microscopy. The Plant Journal, 36(2):280–290, 2003.
- [8] Jurek W Dobrucki, Dorota Feret, and Anna Noatynska. Scattering of exciting light by live cells in fluorescence confocal imaging: phototoxic effects and relevance for frap studies. Biophysical journal, 93(5):1778–1786, 2007.
- [9] Richard M Lasarow, R Rivkah Isseroff, and Edward C Gomez. Quantitative in vitro assessment of phototoxicity by a fibroblast-neutral red assay. Journal of investigative dermatology, 98(5):725–729, 1992.
- [10] Joana Sarah Grah. Methods for automatic mitosis detection and tracking in phase contrast images. Master’s thesis, WWU - University of Münster, 2014.
- [11] Erik Meijering, Oleh Dzyubachyk, Ihor Smal, and Wiggert A van Cappellen. Tracking in cell and developmental biology. In Seminars in cell & developmental biology, volume 20, pages 894–902. Elsevier, 2009.
- [12] M. E. Ambühl, C. Brepsant, J.-J. Meister, A. B. Verkhovsky, and I. F. Sbalzarini. High-resolution cell outline segmentation and tracking from phase-contrast microscopy images. Journal of microscopy, 245(2):161–170, 2012.
- [13] O. Debeir, P. Van Ham, R. Kiss, and C. Decaestecker. Tracking of migrating cells under phase-contrast video microscopy with combined mean-shift processes. Medical Imaging, IEEE Transactions on, 24(6):697–711, 2005.
- [14] S. Huh, D. F. E. Ker, R. Bise, M. Chen, and T. Kanade. Automated mitosis detection of stem cell populations in phase-contrast microscopy images. Medical Imaging, IEEE Transactions on, 30(3):586–596, 2011.
- [15] K. Thirusittampalam, J. Hossain, and P. F. Whelan. A novel framework for cellular tracking and mitosis detection in dense phase contrast microscopy images. IEEE Transactions on Biomedical Engineering, 17(3):642–653, 2013.
- [16] Kang Li, Eric D Miller, Mei Chen, Takeo Kanade, Lee E Weiss, and Phil G Campbell. Cell population tracking and lineage construction with spatiotemporal context. Medical image analysis, 12(5):546–566, 2008.
- [17] Johannes Schindelin, Ignacio Arganda-Carreras, Erwin Frise, Verena Kaynig, Mark Longair, Tobias Pietzsch, Stephan Preibisch, Curtis Rueden, Stephan Saalfeld, Benjamin Schmid, et al. Fiji: an open-source platform for biological-image analysis. Nature methods, 9(7):676–682, 2012.
- [18] Fabrice De Chaumont, Stéphane Dallongeville, Nicolas Chenouard, Nicolas Hervé, Sorin Pop, Thomas Provoost, Vannary Meas-Yedid, Praveen Pankajakshan, Timothée Lecomte, Yoann Le Montagner, et al. Icy: an open bioimage informatics platform for extended reproducible research. Nature methods, 9(7):690–696, 2012.
- [19] Kevin W Eliceiri, Michael R Berthold, Ilya G Goldberg, Luis Ibáñez, Bangalore S Manjunath, Maryann E Martone, Robert F Murphy, Hanchuan Peng, Anne L Plant, Badrinath Roysam, et al. Biological imaging software tools. Nature methods, 9(7):697–710, 2012.
- [20] M. Möller, M. Burger, P. Dieterich, and A. Schwab. A framework for automated cell tracking in phase contrast microscopic videos based on normal velocities. Journal of Visual Communication and Image Representation, 25(2):396–409, 2014.
- [21] P. V. C. Hough. Method and means for recognizing complex patterns, December 18 1962. US Patent 3,069,654.
- [22] R. O. Duda and P. E. Hart. Use of the hough transformation to detect lines and curves in pictures. Communications of the ACM, 15(1):11–15, 1972.
- [23] G. Aubert and P. Kornprobst. Mathematical problems in image processing: partial differential equations and the calculus of variations, volume 147. Springer, 2006.
- [24] T. F. Chan and J. Shen. Image processing and analysis: variational, PDE, wavelet, and stochastic methods. Siam, 2005.
- [25] S. Osher and J. A. Sethian. Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations. Journal of computational physics, 79(1):12–49, 1988.
- [26] V. Caselles, F. Catté, T. Coll, and F. Dibos. A geometric model for active contours in image processing. Numerische mathematik, 66(1):1–31, 1993.
- [27] Laurent D Cohen. On active contour models and balloons. CVGIP: Image understanding, 53(2):211–218, 1991.
- [28] S. Kichenassamy, A. Kumar, P. Olver, A. Tannenbaum, and A. Yezzi. Gradient flows and geometric active contour models. In Computer Vision, 1995. Proceedings., Fifth International Conference on, pages 810–815. IEEE, 1995.
- [29] N. Paragios and R. Deriche. Geodesic active contours and level sets for the detection and tracking of moving objects. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 22(3):266–280, 2000.
- [30] M. Kass, A. Witkin, and D. Terzopoulos. Snakes: Active contour models. International journal of computer vision, 1(4):321–331, 1988.
- [31] V. Caselles, R. Kimmel, and G. Sapiro. Geodesic active contours. International journal of computer vision, 22(1):61–79, 1997.
- [32] T. F. Chan and L. A. Vese. Active contours without edges. Image processing, IEEE transactions on, 10(2):266–277, 2001.
- [33] G. Bertrand. Simple points, topological numbers and geodesic neighborhoods in cubic grids. Pattern recognition letters, 15(10):1003–1011, 1994.
- [34] X. Han, C. Xu, and J. L. Prince. A topology preserving level set method for geometric deformable models. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 25(6):755–768, 2003.
- [35] T. Y. Kong and A. Rosenfeld. Digital topology: introduction and survey. Computer Vision, Graphics, and Image Processing, 48(3):357–393, 1989.
- [36] D. Adalsteinsson and J. A. Sethian. A fast level set method for propagating interfaces. Journal of computational physics, 118(2):269–277, 1995.
- [37] Luca Calatroni, Cao Chung, Juan Carlos De Los Reyes, Carola-Bibiane Schönlieb, and Tuomo Valkonen. Bilevel approaches for learning of variational imaging models. arXiv preprint arXiv:1505.02120, 2015.
- [38] Karl Kunisch and Thomas Pock. A bilevel optimization approach for parameter learning in variational models. SIAM Journal on Imaging Sciences, 6(2):938–983, 2013.
- [39] Asako Sakaue-Sawano, Hiroshi Kurokawa, Toshifumi Morimura, Aki Hanyu, Hiroshi Hama, Hatsuki Osawa, Saori Kashiwagi, Kiyoko Fukami, Takaki Miyata, Hiroyuki Miyoshi, et al. Visualizing spatiotemporal dynamics of multicellular cell-cycle progression. Cell, 132(3):487–498, 2008.
- [40] Siang-Boon Koh, Patrice Mascalchi, Esther Rodriguez, Yao Lin, Duncan I. Jodrell, Frances M. Richards, and Scott K. Lyons. A quantitative fastfucci assay defines cell cycle dynamics at a single-cell level. Journal of Cell Science, 130(2):512–520, 2017.
- [41] Siang-Boon Koh, Aurélie Courtin, Richard J Boyce, Robert G Boyle, Frances M Richards, and Duncan I Jodrell. Chk1 inhibition synergizes with gemcitabine initially by destabilizing the dna replication apparatus. Cancer research, 75(17):3583–3595, 2015.
- [42] Anthony Tighe, Victoria L Johnson, and Stephen S Taylor. Truncating apc mutations have dominant effects on proliferation, spindle checkpoint control, survival and chromosome stability. Journal of cell science, 117(26):6339–6353, 2004.
- [43] Paul Jaccard. La distribution de la flore dans la zone alpine. 1907.
- [44] M-P Dubuisson and Anil K Jain. A modified hausdorff distance for object matching. In Pattern Recognition, 1994. Vol. 1-Conference A: Computer Vision & Image Processing., Proceedings of the 12th IAPR International Conference on, volume 1, pages 566–568. IEEE, 1994.
- [45] Fjedor Gaede. Segmentation and tracking of cells in complete image sequences. Bachelor’s thesis, WWU - University of Münster, 2015.
- [46] Jonas Geiping. Comparison of topology-preserving segmentation methods and application to mitotic cell tracking. Bachelor’s thesis, WWU - University of Münster, 2014.

## Appendix A Supplementary Data

Parameter | Description | MIA PaCa-2 | HeLa Aur A | T24 |

radiusMin | Minimum radius of mitotic cell | 10 | 10 | 10 |

radiusMax | Maximum radius of mitotic cells | 20 | 25 | 20 |

sensitivity | The higher, the more circular objects are detected | 0.8 | 0.7 | 0.7 |

mitosisThreshold | Maximum mitosis duration | 50 | 25 | 25 |

Weight for normal velocity term inside of cell | 1 | 0.5 | 5 | |

Weight for normal velocity term in the background | 1 | 0.1 | 5 | |

Weight for length regularisation (smoothness) | 10 | 8 | 17.5 | |

Weight for local standard deviation term | 10 | 12 | 17.5 | |

g_adj_low | Lower bound for rescaling of local standard deviation image | 0.08 | 0.05 | 0.08 |

g_adj_high | Upper bound for rescaling of local standard deviation image | 0.12 | 0.20 | 0.12 |

Weight for area regularisation | 1 | 1 | 10 | |

timeStep | Time step in gradient descent equation | 1 | ||

maxIterations | Maximum number of iterations for tracking contour evolution | 5000 | 2500 | 5000 |

phiUpdate | Frequency of reinitialisation of level-set function with signed distance function | 50 | 10 | 50 |

Parameter in regularisation of gradient magnitude | 0.0001 | |||

Parameter in regularisation of Dirac delta function | 2 |

MiA | CV | GAC | ||

1 | JSC | 0.6103 | 0.7056 | 0.4960 |

MHD | 2.1713 | 0.9120 | 2.5202 | |

2 | JSC | 0.5736 | 0.6276 | 0.3412 |

MHD | 3.1590 | 1.3452 | 4.8249 | |

3 | JSC | 0.6093 | 0.4678 | 0.3561 |

MHD | 2.1612 | 3.0814 | 5.5271 | |

4 | JSC | 0.6741 | 0.4493 | 0.2848 |

MHD | 1.7944 | 3.0597 | 6.5961 | |

5 | JSC | 0.4428 | 0.4243 | 0.2608 |

MHD | 4.5354 | 2.8466 | 6.8695 | |

6 | JSC | 0.7133 | 0.5030 | 0.3716 |

MHD | 0.8024 | 3.3760 | 6.0406 | |

7 | JSC | 0.6623 | 0.3160 | 0.4835 |

MHD | 1.6662 | 4.3808 | 3.9340 | |

8 | JSC | 0.5402 | 0.4014 | 0.4541 |

MHD | 3.5367 | 3.2165 | 4.6687 | |

9 | JSC | 0.5417 | 0.1597 | 0.4175 |

MHD | 2.0813 | 6.5744 | 5.4570 | |

10 | JSC | 0.6877 | 0.2445 | 0.4496 |

MHD | 0.7584 | 5.9506 | 4.7614 | |

avg | JSC | 0.6050 | 0.4299 | 0.3915 |

MHD | 2.2666 | 3.4743 | 5.1200 |