2021

[1]Jean Ollion [2]Maxime Deforet

1]SABIlab, Die, 26150, France 2]Laboratoire Jean Perrin, Sorbonne Université,
Centre National de la Recherche Scientifique, Paris, 75005, France
3]Laboratoire PhysicoChimie Curie UMR168, Institut Curie, Paris Sciences et
Lettres, Centre National de la Recherche Scientifique, Sorbonne Université,
Paris, 75005, France

# DistNet2D: Leveraging long-range temporal information for efficient
segmentation and tracking

<EMAIL_ADDRESS>Martin Maliet  Caroline Giuglaris  Élise Vacher
<EMAIL_ADDRESS>[  [  [

###### Abstract

Extracting long tracks and lineages from videomicroscopy requires an extremely
low error rate, which is challenging on complex datasets of dense or deforming
cells. Leveraging temporal context is key to overcoming this challenge. We
propose DistNet2D, a new deep neural network (DNN) architecture for 2D cell
segmentation and tracking that leverages both mid- and long-term temporal
information. DistNet2D considers seven frames at the input and uses a post-
processing procedure that exploits information from the entire video to
correct segmentation errors. DistNet2D outperforms two recent methods on two
experimental datasets, one containing densely packed bacterial cells and the
other containing eukaryotic cells. It is integrated into an ImageJ-based
graphical user interface for 2D data visualization, curation, and training.
Finally, we demonstrate the performance of DistNet2D on correlating the size
and shape of cells with their transport properties over large statistics, for
both bacterial and eukaryotic cells.

###### keywords:

Segmentation, Tracking, Deep-Learning, Images, Temporal Information

##  1 Introduction

Automated image analysis is revolutionizing the study of cells, enabling
scientists to measure their position, velocity, shape, and fluorescent signal
intensity. When cells are motile, cell segmentation (localizing cell
boundaries) and tracking (associating cells in consecutive images) are
combined to follow cell properties over time boquet2021bioimage . This data
has already revealed fundamental mechanisms in cell motion, especially in
vitro friedl2009collective ; paul2017cancer ; sengupta2021principles ;
wadhwa2022bacterial . Further improving performance of cell segmentation and
tracking will enable us to understand long-term correlation and extend
analysis to more and more complex systems. However, extracting long lineages
demands an extremely low error rate, which can be difficult to achieve with
complex datasets.

Challenges commonly found in biological images stem from: (i) Cell morphology:
The shape and size of cells can differ from one another or be very similar
within a cell population, and this can change over time, making cell
morphology unreliable as a means to segment and track cells accurately. (ii)
Cell boundaries: The boundaries between cells can sometimes be difficult to
distinguish, as they often touch each other. (iii) Cell movement: The movement
of cells can be drastically different from one to the next, or in time, for
instance with abrupt changes of direction of motion. (iv) Cellular events:
Cellular events, such as mitosis or apoptosis, alter lineages and prove
difficult to track using conventional general-object tracking methods. For
such complex datasets, improving the analysis of individual images is
insufficient to overcome the challenges of cell segmentation and tracking. To
increase performance, it is necessary to leverage temporal context by
considering multiple frames simultaneously.

###  1.1 Cell Segmentation

The main challenge in cell segmentation is not only foreground-background
classification, but also the separation of adjacent cells. The emergence of
deep neural networks (DNN) has greatly improved the efficiency and robustness
of segmentation methods. Early DNN-based methods were mostly based on pixel-
wise classification ronneberger2015u , which resulted in sub-efficient
separation of adjacent cells. Another popular family of methods involve two
sequential DNNs: a region proposal network that generates axis-aligned
bounding box candidates followed by a classification network that filters and
classifies them ren2015faster . Although this type of method is very efficient
on natural images, it is complex to train and has been shown to be unable to
cope with certain cellular shapes schmidt2018cell . Other methods do not
directly segment cell instances but predict (in this article, the output
computed by a DNN is referred to as prediction) pixel-wise features as proxy
for segmentation that are subsequently fed to a clustering algorithm:
payer2019segmenting  proposed multi-dimensional embedding with a loss function
that pushed dissimilarity between neighbors, naylor2018segmentation  proposed
to predict the Euclidean distance map (EDM), fed to a watershed algorithm.
Compared to a binary probability map, EDM has the advantage of emphasizing the
boundary between touching cells while being independent of morphology. A
popular method that is similar (but not equivalent) predicts radial distance
between center and boundaries at predefined angles, which limits the
application to convex objects schmidt2018cell . neven2019instance  proposed an
efficient method that jointly predicts, for each cell pixel, an offset vector
pointing to the cell center and a clustering bandwidth. Similarly,
stringer2021cellpose  predicts a normalized offset vector pointing to the cell
center. In case of filamentous bacteria, this method tends to produce over-
segmentation. This problem was reduced in cutler2022omnipose  by predicting
the EDM along with an offset vector pointing to the cell medial axis
(skeleton), defined by the local maxima of the EDM.

###  1.2 Cell Tracking

The most straightforward approach to cell segmentation and tracking runs in
two independent successive steps: object detection followed by temporal
association of detected objects. The recent method DeLTA 2.0 o2022delta  uses,
for the tracking step, a classification neural network to predict the next
cell for each cell. However, because predictions are made independently for
each cell, this method is likely to produce inconsistent results. A two-step
approach can allow long-term temporal information to feed the tracking
algorithm, as segmentation enables data compaction. Notably, ben2022graph
uses graph neural network to model the entire time-lapse sequence, resulting
in very effective tracking. The main drawback of the two-step approach is that
it is directly limited by the accuracy of the segmentation step. In difficult
cases such as high density of similar cells, even a trained expert requires
temporal context to perform accurate segmentation.

###  1.3 Combined Segmentation and Tracking

Temporal information can be leveraged by combining segmentation and tracking
into a single operation. Several recent methods simultaneously train a
bounding box detector with a tracker that associates bounding boxes candidates
between successive frames voigtlaender2020siam ; chen2021celltrack . In the
context of cell tracking, one limitation of this kind of method is that they
have restricted access to the spatial context around the bounding boxes (such
as the position of neighboring objects), which is crucial when cells have
similar aspects. An emerging trend is the prediction of the displacement
vector of each cell between two successive frames as a proxy for tracking
Hayashida2020MPM ; ollion2020distnet ; loeffler2022 ; hayashida2022consistent
, along with a proxy for segmentation or detection. The actual association of
cells is performed in a post-processing step. One advantage of this strategy
is that it enables simultaneous segmentation and tracking of all cells present
in a time window using a single DNN, which likely yields more consistent
results for both tasks. It is noteworthy that Hayashida2020MPM ;
hayashida2022consistent  do not segment cells but only detect their centers.
ollion2020distnet  introduced an attention layer vaswani2017attention  in the
neural network, and showed that it captures long-range spatial information, in
the one-dimensional case of bacterial cells growing in a microfluidic device.
loeffler2022  use a DNN architecture that performs segmentation independently
for each frame and thus cannot leverage temporal context for segmentation.
However, several works have shown that performing joint segmentation and
tracking improves segmentation by leveraging sequential information
payer2019segmenting ; ollion2020distnet . Due to memory limitations, these
methods can only use a small temporal neighborhood, e.g. Hayashida2020MPM ;
ollion2020distnet ; loeffler2022  use pairs of successive frames (t,t+1). More
recently, hayashida2022consistent  have shown that tracking and detection
performance can be improved by using a larger temporal neighborhood of six
frames as well as a carefully designed loss function that penalizes
inconsistencies between detection and tracking.

In this work, we describe DistNet2D, a novel 2D cell segmentation and tracking
method, with a carefully designed DNN architecture that leverages mid- and
long-term temporal context for both tasks. Mid-term temporal context is
incorporated at the input of the DNN: our method typically considers a
15-frame time window, but this size is adaptable to the features of the
dataset and can be much wider if needed. Long-term temporal context is
incorporated through a post-processing procedure that uses information from
the whole video to correct segmentation errors. We compare DistNet2D to two
recent methods (DeLTA 2.0 o2022delta , EmbedTrack loeffler2022 ) on two
experimental datasets that we publish along with this work: a dataset
containing phase contrast images of dense communities of motile bacterial
cells, and a dataset of fluorescence images of adherent migrating eukaryotic
cells. We also adapted the graphical user interface of BACMMAN software
ollion2019high  for 2D data visualization, curation, and training. Finally, we
demonstrate how DiSTNet2D’s performance enables us to correlate the size and
shape of cells with their motion properties over large statistics, for both
bacterial and eukaryotic cells.

##  2 Results

Following the work of ollion2020distnet ; hayashida2022consistent ;
loeffler2022 , we developed a method that performs segmentation and tracking
simultaneously with a single DNN. This strategy has several advantages over
methods that perform the tasks independently. First, it leverages temporal
information for segmentation, improving the accuracy of the results. Second,
it is easier to train and use a single DNN than two separate networks. Our
method is based on a novel DNN architecture that incorporates a sequence of
operations designed to blend the information gathered from the different input
frames, enabling the use of this information for both segmentation and
tracking (see details in Online Methods). Specifically, several frames are fed
to the DNN, which predicts proxies for both segmentation and tracking (Figure
1A). Using a larger time window is expected to increase temporal consistency,
but the number of considered frames is limited by GPU memory. We chose to
consider seven frames: three frames before and three frames after the current
frame. To enable the DNN to use information at a longer time range without
exceeding its memory capacity, we distributed the seven frames across a larger
range, by spacing them apart (see Figure S1 for a diagram). The gap between
considered frames depends on the dataset, in this work we used values of one
and three frames (depending on the dataset). This strategy is referred to as
temporal subsampling.

Figure 1: Method overview. A: Output of the model for a given frame pair. For
each frame, EDM is the map of the Euclidean distance to the edge of each cell,
GDCM is the map of the geodesic distance to the center. For each pair of
frames, in each direction, dX and dY are the cell center displacements from
previous frame for each axis, $P(\textit{Link multiplicity}=0)$ and
$P(\textit{Link multiplicity}>1)$ are the probabilities that the link
multiplicity is zero (no linked cell in the other frame) and strictly greater
than one (several linked cells in the other frame), respectively. In this
dataset, $P(\textit{Link multiplicity}>1)$ is null because it contains no
mitosis or merging cells. Note that only one frame pair ($t$, $t+1$) is
represented, but the model inputs a larger temporal neighborhood and predicts
these maps for more frame pairs. The output includes both forward
($t\rightarrow t+1$) and backward ($t+1\rightarrow t$) predictions. B:
Segmentation procedure: A watershed transform is applied to the EDM using
regional maxima as seeds, which likely produces over-segmentation. Gaussian
function is applied to the predicted GDCM and a watershed algorithm on the
Laplacian transform is used to detect centers, which are used to reduce over-
segmentation (see main text for details). C: Tracking procedure: The centers
of each cell at $t$ are shifted by the predicted displacement between $t$ and
$t+1$ (dX and dY). Each cell at $t$ is associated with the cell at $t+1$ in
which the shifted center falls. Images are from dataset PhC-C2DH-PA14. Method
overview for dataset Fluo-C2DH-HBEC is available in Figure S2.

###  2.1 Segmentation

The network predicts two complementary proxies for segmentation, the Euclidean
distance map (EDM), which aims at identifying the cell shape, and the geodesic
distance to the center map (GDCM), which aims at identifying the cell center.
These proxies are predicted within each cell and then combined to produce
objects’ contours. More formally: let $B$ be the background, $F$ the
foreground and $C_{j}$ the $j^{th}$ cell ($F=\bigcup_{j}C_{j}$), $c_{j}$ its
center, $d$ the Euclidean distance function and $d_{g}$ the geodesic distance
function (see section 3.6); for each pixel $i$:

$EDM_{i}=\left\\{\begin{array}[]{ll}min(d(i,B),d(i,F\setminus C_{j}))\text{ if
}i\in C_{j}\\\ -1\text{ if }i\in B\\\ \end{array}\right.$
$GDCM_{i}=\left\\{\begin{array}[]{ll}d_{g}(i,c_{j})\text{ if }i\in C_{j}\\\
0\text{ if }i\in B\\\ \end{array}\right.$

The medoid center of the cells is used in order to ensure it is contained in
the cell, even for non-convex shapes. For simplicity, it will be referred to
as center in this work.

EDM-based segmentation is efficient even on non-convex cell morphologies, such
as bent bacterial cells. It is performed by applying a watershed algorithm on
EDM. Watershed is naturally limited to positive values (as the background is
set to -1) and seeded from regional maxima of the EDM, which can easily
produce over-segmentation, especially in long cells or cells with complex
shapes that may contain several local maxima.

To reduce over-segmentation, we combined contours with predicted centers.
Centers are segmented by performing a watershed algorithm on the Laplacian
transform of the image that results from the Gaussian function applied to GDCM
(see section 3.6 for details). Two segmented regions in contact are merged if
either one of them or both do not contain a segmented center, or if the ratio
of intensity amplitude of the centers is below a user-defined threshold (see
Figure 1B).

Moreover, we also observed that predicting a unique center per cell improves
EDM prediction especially in distinguishing neighboring cells instances.

###  2.2 Tracking

Tracking is performed using the prediction of the displacement of the cell
centers along the X and Y axis that occurs between two frames. The center of
each cell at $t$ is shifted by its predicted displacement between $t$ and
$t+1$; if the shifted center falls into a segmented cell at $t+1$, then the
two cell instances are associated (see Figure 1C). This is similar to the
procedure used in loeffler2022 .

To assist the tracking procedure and manage more complex cases, we also
predict the link multiplicity category in both the forward ($t\rightarrow
t+1$) and backward ($t+1\rightarrow t$) directions, accounting for the
expected number of links for each cell (Figure S3). The possible values for
forward link multiplicity are: no next cell (the cell will leave the field of
view, or will die), one next cell (regular case, or the cell will fuse with
another cell), multiple next cells (the cell will divide). The possible values
for backward link multiplicity are: no previous cell (the cell has entered the
field of view, or just appeared), one previous cell (regular case, or the cell
has just divided), multiple previous cells (the cell has just fused).
Formally, we predict each time three probability maps: $P(\textit{Link
Multiplicity}=0)$, $P(\textit{Link Multiplicity}=1)$, $P(\textit{Link
Multiplicity}>1)$, summing to 1. We assign the link multiplicity category to
each cell as the multiplicity with the highest median probability within the
cell.

Cells that are predicted to have a single next cell are linked by forward
tracking using the predicted forward displacement. Cells with multiple next
cells (either because of over-segmentation at $t+1$ and not at $t$, under-
segmentation at $t$ and not at $t+1$, or a predicted division event) remain
unlinked after forward tracking. Backward tracking is then applied on unlinked
cells that are predicted to have a single previous cell, using predicted
backward displacement (Figure S4).

Forward and backward tracking allow to assign both merge links —in which
several cells at $t$ are associated to a single cell at $t+1$— and split links
—in which several cells at $t+1$ are associated to a single cell at $t$. When
merge links or split links are not confirmed by the link multiplicity
category, it means they arise from over/under segmentation errors and they
will be corrected in the post-processing stage (see section 2.3). Identifying
merge and split links also enables a finer definition of metrics and a more
accurate diagnosis of origin of errors (incorrect segmentation vs incorrect
linking), see section 2.6 for details.

###  2.3 Post-processing: Segmentation correction

A set of rules was designed to correct over/under segmentation errors using
the tracking information. We especially observed such errors in the PhC-C2DH-
PA14 dataset, which consists of high-speed acquisition videos (typically 100
frames-per-second for a few seconds) of motile rod-shaped bacteria that divide
typically every hour. The invagination of the cell membrane at the center of
the mother cell is the last stage of the cell cycle before the separation of
the two daughter cells. In phase contrast imaging, this invagination appears
as a bright region within the cell body, which is similar to the bright area
that connects two separate cells when they are in contact (Figure S5). It is
virtually impossible to determine from a single frame whether an object
represents a late-dividing long single cell or two adjacent short cells. The
DNN architecture already helps avoid such errors by considering a mid-term
temporal context (typically fifteen frames).

To consider long-term temporal information (which can cover an extended period
of time, up to the duration of the entire video), we analyze the trajectories
obtained during the tracking step. We focus on the merge and split links.
Merge and split links consistent with predicted link multiplicity are
considered mitosis or fusion, and are left untouched. In contrast, merge and
split links that contradict with predicted link multiplicity are suspected
errors. We treat them using the following general principle: if an object has
always been seen as one cell, it should remain as one; however, if it was
detected as two distant cells at any point in the past or future, this
indicates that it should be considered two cells (Figure 2A). This approach is
based on the assumption that errors are rare and can be corrected by looking
at errorless past and future. The high efficiency of our DNN-based combined
segmentation and tracking algorithm supports this assumption.

In practice, for each merge link, we check if all cells before the link are in
contact (two objects are considered to be in contact if the distance between
their contours is lower than a threshold ; for rod-shaped bacteria, an
alignment criterion is also used). If they are, then we merge all of the
cells. Otherwise, we split the objects following the link by applying a
watershed transform on the EDM (Figure 2A-i). Cell fragments are linked to the
previous objects using the same procedure as in section 2.2. If the watershed
algorithm generates more fragments than there are previous objects, fragments
linked to the same previous object are merged. Similarly, for each split link,
if all cells after the link are in contact, then we merge all of the cells.
Otherwise, we split the objects before the link (Figure 2A-ii). Common
examples are depicted in Figure 2B-D.

Figure 2: Post-processing uses temporal information on large timescale to
correct segmentation errors. A: Diagrams presenting the procedure to correct
wrong links. B: Illustrative example of two distinct cells that are
transiently detected as one object (under-segmented in frame $t=65$). Because
objects involved in the merge link are sometimes seen separated, we can assume
the object in frame $t=65$ should be split. C: Illustrative example of one
cell that is transiently detected as two objects (over-segmented in frame
$t=79$). Because the cell is detected as a single object throughout the entire
video (200 frames), except for one frame in which it is seen as two objects,
we can assume that the two objects should be merged into one. D: Illustrative
example of a more complex error that implies more lineages. In frame $t=145$,
one cell is under-segmented and one cell is over-segmented. Here again,
temporal information at the scale of the entire video enables us to correct
the segmentation errors.

###  2.4 Model architecture

Figure 3: Model architecture. The encoder is fed by successive frames (green,
blue and red rectangles) and produces encoded features (green, blue and red
cubes). Features are processed in pairs (corresponding to successive frames)
by the pair blender module, which produces feature pairs. Encoded features and
feature pairs are blended together by the blender module (see Figure S6 for
details). The segmentation extractor generates three segmentation features
corresponding to each frame for both EDM and GDCM, that are decoded by two
distinct decoders to produce images of the same size as the input image.
Likewise, the tracking extractor generates two tracking features corresponding
to each frame pair for both displacement and link multiplicity, that are
decoded by two distinct decoders. For simplicity only three frames have been
represented but we considered seven frames in this work, and only one
segmentation decoder and one tracking decoder are represented instead of two.

The DNN has an encoder-decoder architecture, with a single encoder and one
decoder per output type (EDM, GDCM, displacement and link multiplicity). The
encoder and decoders are shared between frames for better training efficiency.
For segmentation outputs (EDM and GDCM), one prediction per frame is made,
whereas for tracking outputs (displacement and link multiplicity) one
prediction per frame pair is made. For a DNN time window of size
$(N-3)\delta+3$ centered on frame $t$, the $N$ considered frames are the three
central frames ($t-1$, $t$, $t+1$), and $m=(N-3)/2$ frames on each side of the
central frames ($t-1-m\delta$, …, $t-1-2\delta$, $t-1-\delta$ and
$t+1+\delta$, $t+1+2\delta$, …, $t+1+m\delta$) (Figure S1). $2N-4$ frame pairs
are defined as follow (Figure S7):

  * •

$N-1$ frame pairs between consecutive considered frames, for short-range
displacements,

  * •

$N-3$ frame pairs between the central frame (frame $t$) and each other frame
except frames $t-1$ and $t+1$, for mid-range displacements.

Figure 3 displays the global architecture. The detailed architecture of each
box is described in Online Methods 3.1. Encoder and decoders are mainly
composed of residual blocks of two successive convolutions. Between the
encoder and the decoders, we introduced a blending module —a sequence of
operations that blends the encoded features of all frames together— and then
extract one feature per frame (resp. frame pair) that are fed to segmentation
(resp. tracking) decoders. Extraction sequences simply contain a distinct
convolution per frame (resp. frame pair). Resulting tensors are combined with
encoded features (resp. feature pairs). We define the combine operation as a
1x1 convolution applied to the concatenation of two tensors.

The blending-extraction sequence makes the information from the whole time
window available for the prediction of each output at each frame. This
contrasts with loeffler2022  in which segmentation is performed independently
on each other frame.

This architecture has the advantage of allowing proxy predictions to be made
only for the central frame and the frame pair (t, t+1) during the prediction
phase, which improves speed and reduces memory consumption.

###  2.5 Software

The software associated with our method is BACMMAN ollion2019high , an ImageJ
schneider2012nih  plugin that was initially developed for analysis of
bacterial cells growing in microchannels, with displacement along the
microchannel axis. Such data was naturally displayed on kymographs, in which
the horizontal axis represents time. In order to display 2D data, we added the
hyperstack visualization mode (see Figure S8), in which lineage information is
displayed as colored contours. All the features of the graphical user
interface of BACMMAN such as interactive navigation through images, manual
curation, two-way interplay with R/Python for statistical analysis, are thus
available for 2D data. BACMMAN was also augmented with new features:
generation of training sets as well as DistNet2D training and prediction can
now be performed directly from the software. BACMMAN also provides a command-
line interface, enabling its use on a computational cluster.

###  2.6 Evaluation Metrics

Objective metrics for segmentation and tracking were previously introduced in
matula2015cell . In that work, cell tracking results were represented using an
acyclic oriented graph, in which nodes corresponded to the detected cells and
edges represented links (i.e temporal relations) between them. Metrics were
based on the number of operations required to transform the result graph into
the reference graph. Those operations were {split/delete/add} a node and
{delete/add/change the semantics of} an edge (e.g.: a change between a split
link and a normal link). Correspondence between a reference (R) and a result
(S) segmented object were established using the following criterion: $\lvert
R\cap S\lvert\>>0.5\cdot\lvert R\lvert$, which implied that each reference
cell could correspond to one result cell at most. We consider this a
limitation because it does not allow to take over-segmentation into account;
over-segmented cells were thus systematically considered as false positives.
Instead, we used the following criterion: $\lvert R\cap S\lvert\>>0.5\cdot
min(\lvert R\lvert,\lvert S\lvert)\;\textit{OR}\;\lvert R\cap S\lvert\>>C$
with $C$ a user-defined constant that is typically $50\%$ of the average
reference cell size. The second term accounts for cases of under-segmentation
and partial overlap with ground truth, where the relative overlap is too low
but the absolute overlap is significant.

This approach identifies four types of segmentation errors: false positives
(result cells with no reference counterpart), false negatives (reference cells
with no result counterpart), over-segmentation (when $N$ result cells match a
given reference cell, $N-1$ over-segmentation are counted) and under-
segmentation (when $N$ reference cells match a given result cell, $N-1$ under-
segmentation are counted).

We also observed that the procedure proposed in matula2015cell  does not fully
distinguish between tracking and segmentation errors: for instance, over-
segmentation of one object into two parts is counted as a false positive
segmentation error as well as a false negative link. However, if the over-
segmented cell parts were all linked to the correct cell(s), no tracking error
should be counted (Figure S9). We thus developed a procedure inspired by
matula2015cell  to identify tracking errors that are independent of
segmentation errors. In other words, our procedure evaluates the tracking
efficiency per se, given the segmentation errors.

To do so, for each frame pair ($t$, $t+1$), we transform the nodes of the
result graph so that they match with the nodes of the reference graph by
applying four successive operations (Figure S9): splitting under-segmented
cells at $t+1$, splitting under-segmented cells at $t$, merging over-segmented
cells at $t$ and merging over-segmented cells at $t+1$. At each split/merge
operation, links are propagated to the resulting nodes. In the case of
splitting, if this implies linking $M$ nodes at $t$ to $N$ nodes at $t+1$,
where $M>1$ and with $N>1$, links are determined by a simple linear assignment
algorithm that minimizes the distance between cell centers. In the case of
merging, all links are simply added to the resulting node. After these
transformations, all nodes of the transformed result graph correspond to a
single node in the reference graph, except for false positives, which have no
counterparts in the result graph. This enables counting false positive and
false negative links: the former are links found in the transformed result
graph and not in the reference graph, except links automatically added by
splitting an under-segmented object. The latter are links from the reference
graph that are not found in the transformed result graph and that do not
involve false negative objects. Our choice of link propagation in case of
merging during graph transformation would miss some false negative links, that
are thus added to the count: in case a result cell was merged at frame $t+1$
but has no link with cells at $t$ and the corresponding reference cell is
linked, a false negative link is counted. The same applies for a result cell
merged at frame $t$ that has no link with cells at $t+1$.

This procedure is applied for each pair of frames, but this does not tell us
how these errors are distributed among the different cell lineages. This
information is crucial for evaluating an algorithm, as recovering more error-
free cell lineages can be more useful even if it makes more frame-pair-wise
errors mavska2023cell . We therefore counted the number of error-free cell
lineages, allowing for a user-defined tolerance to the frame at which mitosis
is detected.

Lastly, we noticed that many errors arise from cells that are partially out of
bounds. We believe that these errors can be easily removed automatically and
should not be counted as errors. Therefore, our procedure simply ignores
errors that are related to cells that touch edges or lineages that contain at
least one cell that touches an edge.

###  2.7 Evaluations

Table 1: Comparison of DistNet2D, DeLTA 2.0 and EmbedTrack on datasets
PhC-C2DH-PA14 and Fluo-C2DH-HBEC. Segmentation errors is the sum of false
positive, false negative, under- and over-segmentation, divided by the number
of cells in the ground truth. Tracking errors is the sum of false positive and
false negative links divided by the number of links in the ground truth.
Incomplete lineages is the number of lineages with at least one segmentation
or tracking error divided by the number of lineages in the ground truth. For
datasets PhC-C2DH-PA14 and Fluo-C2DH-HBEC respectively, the total number of
cells is 123,057 and 4,273, the total number of links is 145,306 and 4,890,
and the total number of lineages is 561 and 60. As explained in the main text,
only cells that do not touch edges (and lineages with no cell touching edges)
are taken into account.

Dataset | Method | Segmentation | Tracking | Incomplete | Prediction
---|---|---|---|---|---
|  | errors (%) | errors (%) | lineages (%) | time (s/frame)
2 |  |  |  |  |
PhC-C2DH-PA14 | DeLTA 2.0 | $1.1$ | $0.22$ | $21$ | $22$
EmbedTrack | $1.2$ | $0.14$ | $7.1$ | $2.7$
DistNet2D | $0.24$ | $0.00$ | $0.53$ | $0.51$
2 |  |  |  |  |
Fluo-C2DH-HBEC | DeLTA 2.0 | $1.9$ | $0.92$ | $37$ | $0.84$
EmbedTrack | $0.89$ | $0.51$ | $25$ | $1.4$
DistNet2D | $0.49$ | $0.10$ | $5.0$ | $0.27$

DeLTA 2.0 independently performs segmentation and tracking using two
independent U-Net models, with ad hoc procedures specifically designed for
rod-shaped bacteria, such as a skeletonization of the cell body shape to set a
maximal weight at the center of the cell during training. EmbedTrack
simultaneously performs segmentation and tracking using a single DNN, but it
uses a total time window of only 2 frames (t,t+1) and only forward
predictions.

DistNet2D outperforms DeLTA 2.0 and EmbedTrack on all our metrics:
segmentation errors, tracking errors, and incomplete lineages (Table 1).
Notably, DistNet2D achieved perfect tracking accuracy on the bacterial
dataset. Additionally, DistNet2D ran faster than the other two methods on both
datasets.

Table 2: Ablation experiments on dataset PhC-C2DH-PA14 (the total number of cells is 123,057, links is 145,306, and lineages is 561). The DNN time window of 2 corresponds to a simplified version of the DNN that considers only frames $(t,t+1)$. The last case, which corresponds to DistNet2D, has a DNN time window of 15 frames ($N=7$ considered frames and $\delta=3$ using the subsampling definition presented in Figure S1), included post-processing, and achieved the best performance. Ablations  | Results
---|---
Post-processing |  DNN Time window |  Segmentation Errors (%) |  Tracking Errors (%) |  Incomplete Lineages (%)
No |  2 |  $0.65$ |  $0.00$ |  $10$
No |  15 |  $0.47$ |  $0.00$ |  $5.3$
Yes |  2 |  $0.50$ |  $0.00$ |  $0.89$
Yes |  15 |  $0.24$ |  $0.00$ |  $0.53$

We also evaluated the contribution of several components of our method by
performing ablation experiments (Table 2). Switching the DNN time window from
two frames (a single frame pair ($t,t+1$), as in EmbedTrack and DeLTA 2.0) to
fifteen frames, without post-processing, increased segmentation performance
and notably decreased the number of incomplete lineages. This shows that mid-
term context is leveraged for segmentation and brings temporal consistency.
While the two-frame version without post-processing exhibits fewer
segmentation errors compared to EmbedTrack (0.65% versus 1.2%), it also
results in more incomplete lineages (10% versus 7.1%). This difference arises
from the contrasting under-segmentation to over-segmentation ratios observed
in the two methods, with under-segmentations having a more pronounced impact
on incomplete lineages.

By leveraging temporal information across the entire video, post-processing
effectively corrected suspected errors (defined by merge links and split links
that are not confirmed by predicted link multiplicity), leading to a
substantial reduction in both object-wise segmentation errors and, more
notably, lineage-wise errors, enhancing the overall temporal coherence of
predictions. We acknowledge that post-processing might introduce new errors
that will be propagated over entire tracks, but we consider this is an
acceptable drawback considering its great performance in terms of reducing
incomplete lineages. We also concede that post-processing works well because
the DNN already makes few errors. Post-processing and DNN work
synergistically, as shown by the poor performance when both components are
turned off.

We then tested how DistNet2D’s performances in leveraging long-term
information are affected by frame sampling. First, we varied how the seven
considered frames were spread apart, by playing with the parameter $\delta$,
which changed the range of the DNN time window (Figure 4A). We found that a
wider DNN time window improved the segmentation performance (the tracking
efficiency was already very high, even for small $\delta$). This suggests that
DistNet2D benefits from having access to a longer temporal context. The
accuracy over entire lineages was also improved with greater $\delta$.
However, this effect was eliminated after post-processing, possibly because
the number of remaining incomplete lineages is too small (lower than 5).

To further assess the influence of time sampling, we compared the performance
of DistNet2D, EmbedTrack, and Delta 2.0 on subsampled evaluation datasets
(Figure 4B). While the accuracy of segmentation stayed roughly the same for
all methods across the tested range, the accuracy of tracking was sensitive to
subsampling, both at the level of individual links or entire tracks. However,
DistNet2D remained fairly accurate for tracking, unlike the other two methods
which quickly lost their effectiveness. For example, DistNet2D performed just
as well on the 6-fold subsampled dataset as the other two methods did on the
full dataset. Its robustness to subsampling can be explained by the awareness
of mid-range temporal information and by the random subsampling performed
during data augmentation (see Online Methods 3.3).

Figure 4: Temporal insights of the DistNet2D’s performance, calculated with
dataset PhC-C2DH-PA14. A: Influence of DNN time window in predictions. DNN
time window was varied by changing the gap between considered frames
($\delta$) from 1 to 10, while keeping the number of considered frames
constant ($N=7$). B: Robustness to acquisition subsampling. We compared the
three methods (trained with the complete training dataset) on computationally-
generated time subsampled versions of the evaluation dataset. Note that for
DistNet2D, $\delta$ was set to 1 for all points, which explains the difference
in segmentation errors compared to Table 1, where $\delta=3$.

Overall, the ablation experiments and subsampling experiments confirm that
leveraging temporal information improves the segmentation and tracking
performance.

###  2.8 Showcase of DistNet2D’s performance on biological systems

We demonstrate the potential of DistNet2D by applying it to bacterial and
eukaryotic datasets.

####  2.8.1 Monolayer of bacterial cells

We measured the mean-squared-displacement (MSD) of P. aeruginosa cells at the
surface of an agar gel, at a surface fraction of $\phi=0.719$, where the cell
monolayer appears ”jammed”. As shown in section 2.7, DistNet2D was able to
extract long tracks: 3404 1000-frame tracks out of an average of 3620 cells
that were visible in one field of view. Track duration statistics indicates
very few errors (Table S2). The MSD scales approximately with $t$, confirming
diffusive behavior over 3 decades of time (Figure 5A). At lower surface
fraction ($\phi=0.466$), cells were more motile and a larger fraction of them
left the field of view within the duration of the video (1413 1000-frame
tracks out of an average of 2911 visible cells). Accordingly, the behavior at
short time scales was over-diffusive but remained diffusive at longer time
scales.

We also correlated the length of each cell with its speed. At low density,
short cells moved faster than longer cells (Figure 5B). We hypothesize that
this is a signature of single-cell motility, as viscous drag varies
monotonously with the length of a rod kamdar2022multiflagellarity . This trend
disappeared at high density as cells collectively blocked each other,
regardless of their length.

Figure 5: A: Mean-squared-displacement measured on bacterial cells at low
surface fraction ($\phi=0.466$, blue dots) and high surface fraction
($\phi=0.719$, orange dots). Guide lines have exponent 1 (dashed lines) and 2
(solid line). B: Cell speed as a function of cell length, using the same color
code. Error bars are standard errors of the mean (error bars are hidden behind
dots for the high surface fraction dataset). Each dot is the mean speed of all
cells binned by length, with a bin size of $0.088~{}\mu m$ (one pixel).
Averaging is weighted by the duration of trajectories. The solid lines
represent the number of cells in each bin, weighted by the duration of
trajectories (a value of 1 is attributed to a cell tracked for 1000 frames).
C: Histogram of angles between the velocity vector and the major axis of the
HBEC cell body (obtained by ellipse fitting), for five time intervals (0h-15h,
15h-30h, 30h-45h, 45h-60h, 60h-75h). Each interval includes
N=40,939/53,284/76,305/115,860/138,749 data points (respectively). No
chirality is measured in the data. Inset: average cell density for each time
interval. D: Norm of the velocity vector with respect to the angle between the
velocity vector and the major axis of the HBEC cell body, for the same five
time intervals.

####  2.8.2 Eukaryotic cells

The ability of DistNet2D to extract precise cell contours and long
trajectories enables the correlation of migrating speed and direction with
cellular shape. At low density, HBEC cells typically migrate perpendicular to
their major axis (Figure 5C). At higher density, collisions reorient cells,
leading to a less peaked angle distribution. Interestingly, while the largest
displacements are perpendicular to the cell body major axis at low density,
the trends reverses at high density as some cells move along other cells
(Figure 5D).

##  3 Discussion

DistNet2D introduces several novel components to address the challenges of
segmentation and tracking in bioimages. These components include a novel
segmentation method that combines EDM and GDCM, improving the separation of
adjacent cells, while being able to segment a wide range of cell morphologies.
DistNet2D also employs a novel approach to predicting backward and forward
tracking proxies to handle cell division and fusion events and improve
robustness of tracking to under/over-segmentation errors. Long-range temporal
context is leveraged in a novel post-processing stage that corrects incorrect
merge and split links, relying on the predicted link multiplicity. This stage
strongly reduces the lineage-wise error rate by analyzing entire lineages and
correcting them when necessary. This strategy is efficient even if errors are
also generated at the lineage scale as long as the cell-wise error rate before
post-processing is very low.

We developed a series of carefully chosen innovations for the implementation
of the DistNet2D algorithm. Tracking and segmentation proxies are predicted by
a DNN with a novel architecture designed for leveraging mid-range temporal
information for segmentation, while being optimized for size and training
efficiency. This range is further increased thanks to a gapped input strategy
at no additional GPU memory cost. The loss function was chosen in order to
effectively guide the training process by generating gradients of similar
magnitude between segmentation and tracking proxies, regardless of cell size
and displacement amplitude. Carefully designed data augmentation allows
generalization to diverse imaging conditions without requiring retraining. In
particular, we introduced on-the-fly random frame sub-sampling, which improves
robustness to changes in acquisition rate, but also increases tracking
performances by diversifying displacements during training. Each of these
components plays a crucial role in enhancing the performance of DistNet2D,
making it a powerful tool for analyzing biological processes involving cell
movement.

Following ollion2020distnet , we tried to introduce an attention layer at the
Pair Blender module, but this did not improve performance. This is likely
because the two datasets used in the study only contained short-range
displacements, which could be captured by dilated convolutions. However, an
attention layer may be useful for processing datasets with longer-range
displacements or displacements that depend on location, such as in
microfluidic devices robert2018mutation ; dal2020short .

We demonstrated the performance of DistNet2D on two different datasets: a
bacterial dataset, where cells are densely packed, have similar shape and only
differ in size, and an eukaryotic dataset where cells are sparse, change
shape, and divide. Some methods are specifically designed for a precise type
of datasets, with ad hoc procedures. This is the case for DeLTA 2.0, which was
designed for bacterial datasets. Despite its specificity, DeLTA 2.0 under-
performs compared to DistNet2D. A major strength of DistNet2D is its ability
to leverage both mid- and long-range temporal context, unlike other considered
methods: DeLTA 2.0 performs segmentation and tracking separately and does not
try to leverage temporal information for segmentation. EmbedTrack performs
segmentation and tracking jointly, but its DNN does not blend temporal
information and thus does not make it available for all frames, thus
segmentation is performed on each frame independently. In both DeLTA 2.0 and
EmbedTrack, segmentation does not have access to temporal context.

We believe DistNet2D is general and can be used to segment and track any type
of moving object in a 2D setting: living or inert, with changing or with fixed
shape, with division, with fusion, at any surface density. Any type of imaging
can be used: fluorescence, phase contrast, brightfield, etc. The graphical
user interface, BACMMAN, is an ImageJ plugin that enables training, curation,
manual correction, re-training, distribution of the trained DNN weights to
other labs, and data export as tabular data or as label images. Moreover,
BACMMAN is able to handle multiple classes of object simultaneously, for
instance cell membrane/cell nucleus, cells/foci, head/tail. This could be of
particular interest in microbiology, cell biology, soft matter, active matter,
but also ethology. The expansion of DistNet2D to 3D datasets is a promising
avenue for future exploration. However, the increased memory requirements
associated with processing 3D images necessitate further research to ensure
efficient training.

Like any other supervised DNN-based method, DistNet2D training requires a
training dataset that can be cumbersome to generate. To facilitate de novo
dataset creation, BACMMAN includes a DNN-based segmentation method with very
low annotation requirements ludvikova2023near . Automated tracking can be done
in BACMMAN, that includes several tracking methods, such those in Trackmate
tinevez2017trackmate , an ImageJ plugin that BACMMAN is connected to. The
entire pipeline for creating the training dataset (annotations, DNN-based
segmentation, manual correction of segmentation, tracking, manual correction
of tracking), as well as training of DistNet2D, can be performed within
BACMMAN.

DistNet2D, a novel method for segmentation and tracking of cells in bioimages,
effectively addresses the challenges of segmentation and tracking by
leveraging temporal context and employing carefully designed deep learning
architectures. The integrated graphical user interface, BACMMAN, offers a
comprehensive pipeline that streamlines the process of generating training
datasets, training DNN, and applying DNN-based segmentation and tracking,
making it more accessible and practical for a broad range researchers.
Together, DistNet2D and BACMMAN form a versatile framework for analyzing
biological processes involving cell movement.

Acknowledgments This project has received financial support from the ANR
X-BACAMAT project (ANR-21-CE30-0025). We are grateful to Pascal Silberzan and
Charles Ollion for helpful discussions.

## Declarations

JO and MD were involved in conceptualization, planning and supervised the
work. JO developed and implemented DistNet2D. MM and MD generated the
bacterial datasets. EV and CG generated the eukaryotic datasets. MM analyzed
the bacterial data and EV analyzed the eukaryotic data. MM, MD, and JO
participated in the manual correction of training datasets. JO and MM
performed training, evaluations, and ablation experiments. All authors wrote
the manuscript. JO is the founder and director of the company SABILab. Other
authors declare no competing interest. Code will be available upon
publication. This work is licensed under a Creative Commons Attribution 4.0
International License (CC-BY-4.0).

## References

  *   * (1) Boquet-Pujadas, A., Olivo-Marin, J.-C., Guillén, N.: Bioimage analysis and cell motility. Patterns 2(1) (2021)
  * (2) Friedl, P., Gilmour, D.: Collective cell migration in morphogenesis, regeneration and cancer. Nature reviews Molecular cell biology 10(7), 445–457 (2009)
  * (3) Paul, C.D., Mistriotis, P., Konstantopoulos, K.: Cancer cell motility: lessons from migration in confined spaces. Nature reviews cancer 17(2), 131–140 (2017)
  * (4) SenGupta, S., Parent, C.A., Bear, J.E.: The principles of directed cell migration. Nature Reviews Molecular Cell Biology 22(8), 529–547 (2021)
  * (5) Wadhwa, N., Berg, H.C.: Bacterial motility: machinery and mechanisms. Nature reviews microbiology 20(3), 161–173 (2022)
  * (6) Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18, pp. 234–241 (2015). Springer
  * (7) Ren, S., He, K., Girshick, R., Sun, J.: Faster r-cnn: Towards real-time object detection with region proposal networks. Advances in neural information processing systems 28 (2015)
  * (8) Schmidt, U., Weigert, M., Broaddus, C., Myers, G.: Cell detection with star-convex polygons. In: Medical Image Computing and Computer Assisted Intervention–MICCAI 2018: 21st International Conference, Granada, Spain, September 16-20, 2018, Proceedings, Part II 11, pp. 265–273 (2018). Springer
  * (9) Payer, C., Štern, D., Feiner, M., Bischof, H., Urschler, M.: Segmenting and tracking cell instances with cosine embeddings and recurrent hourglass networks. Medical image analysis 57, 106–119 (2019)
  * (10) Naylor, P., Laé, M., Reyal, F., Walter, T.: Segmentation of nuclei in histopathology images by deep regression of the distance map. IEEE transactions on medical imaging 38(2), 448–459 (2018)
  * (11) Neven, D., Brabandere, B.D., Proesmans, M., Gool, L.V.: Instance segmentation by jointly optimizing spatial embeddings and clustering bandwidth. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 8837–8845 (2019)
  * (12) Stringer, C., Wang, T., Michaelos, M., Pachitariu, M.: Cellpose: a generalist algorithm for cellular segmentation. Nature methods 18(1), 100–106 (2021)
  * (13) Cutler, K.J., Stringer, C., Lo, T.W., Rappez, L., Stroustrup, N., Brook Peterson, S., Wiggins, P.A., Mougous, J.D.: Omnipose: a high-precision morphology-independent solution for bacterial cell segmentation. Nature methods 19(11), 1438–1448 (2022)
  * (14) O’Connor, O.M., Alnahhas, R.N., Lugagne, J.-B., Dunlop, M.J.: Delta 2.0: A deep learning pipeline for quantifying single-cell spatial and temporal dynamics. PLOS Computational Biology 18(1), 1009797 (2022)
  * (15) Ben-Haim, T., Raviv, T.R.: Graph neural network for cell tracking in microscopy videos. In: Computer Vision–ECCV 2022: 17th European Conference, Tel Aviv, Israel, October 23–27, 2022, Proceedings, Part XXI, pp. 610–626 (2022). Springer
  * (16) Voigtlaender, P., Luiten, J., Torr, P.H., Leibe, B.: Siam r-cnn: Visual tracking by re-detection. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 6578–6588 (2020)
  * (17) Chen, Y., Song, Y., Zhang, C., Zhang, F., O’Donnell, L., Chrzanowski, W., Cai, W.: Celltrack r-cnn: A novel end-to-end deep neural network for cell segmentation and tracking in microscopy images. In: 2021 IEEE 18th International Symposium on Biomedical Imaging (ISBI), pp. 779–782 (2021). IEEE
  * (18) Hayashida, J., Nishimura, K., Bise, R.: Mpm: Joint representation of motion and position map for cell tracking. In: The IEEE Conference on Computer Vision and Pattern Recognition (CVPR) (2020)
  * (19) Ollion, J., Ollion, C.: Distnet: Deep tracking by displacement regression: application to bacteria growing in the mother machine. In: International Conference on Medical Image Computing and Computer-Assisted Intervention, pp. 215–225 (2020). Springer
  * (20) Löffler, K., Mikut, R.: Embedtrack—simultaneous cell segmentation and tracking through learning offsets and clustering bandwidths. IEEE Access 10, 77147–77157 (2022)
  * (21) Hayashida, J., Nishimura, K., Bise, R.: Consistent cell tracking in multi-frames with spatio-temporal context by object-level warping loss. In: Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision, pp. 1727–1736 (2022)
  * (22) Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A.N., Kaiser, Ł., Polosukhin, I.: Attention is all you need. Advances in neural information processing systems 30 (2017)
  * (23) Ollion, J., Elez, M., Robert, L.: High-throughput detection and tracking of cells and intracellular spots in mother machine experiments. Nature protocols 14(11), 3144–3161 (2019)
  * (24) Schneider, C.A., Rasband, W.S., Eliceiri, K.W.: Nih image to imagej: 25 years of image analysis. Nature methods 9(7), 671–675 (2012)
  * (25) Matula, P., Maška, M., Sorokin, D.V., Matula, P., Ortiz-de-Solórzano, C., Kozubek, M.: Cell tracking accuracy measurement based on comparison of acyclic oriented graphs. PloS one 10(12), 0144959 (2015)
  * (26) Maška, M., Ulman, V., Delgado-Rodriguez, P., Gómez-de-Mariscal, E., Nečasová, T., Guerrero Peña, F.A., Ren, T.I., Meyerowitz, E.M., Scherr, T., Löffler, K., et al.: The cell tracking challenge: 10 years of objective benchmarking. Nature Methods, 1–11 (2023)
  * (27) Kamdar, S., Ghosh, D., Lee, W., Tatulea-Codrean, M., Kim, Y., Ghosh, S., Kim, Y., Cheepuru, T., Lauga, E., Lim, S., et al.: Multiflagellarity leads to the size-independent swimming speed of bacteria. arXiv preprint arXiv:2212.13632 (2022)
  * (28) Robert, L., Ollion, J., Robert, J., Song, X., Matic, I., Elez, M.: Mutation dynamics and fitness effects followed in single cells. Science 359(6381), 1283–1286 (2018)
  * (29) Dal Co, A., van Vliet, S., Kiviet, D.J., Schlegel, S., Ackermann, M.: Short-range interactions govern the dynamics and functions of microbial communities. Nature ecology & evolution 4(3), 366–375 (2020)
  * (30) Ludvikova, L., Simon, E., Deygas, M., Panier, T., Plamont, M.-A., Ollion, J., Tebo, A., Piel, M., Jullien, L., Robert, L., et al.: Near-infrared co-illumination of fluorescent proteins reduces photobleaching and phototoxicity. Nature Biotechnology, 1–5 (2023)
  * (31) Tinevez, J.-Y., Perry, N., Schindelin, J., Hoopes, G.M., Reynolds, G.D., Laplantine, E., Bednarek, S.Y., Shorte, S.L., Eliceiri, K.W.: Trackmate: An open and extensible platform for single-particle tracking. Methods 115, 80–90 (2017)
  * (32) Charbonnier, P., Blanc-Féraud, L., Aubert, G., Barlaud, M.: Deterministic edge-preserving regularization in computed imaging. IEEE Transactions on image processing 6(2), 298–311 (1997)
  * (33) Friedman, J., Tibshirani, R., Hastie, T.: The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer, ??? (2009)
  * (34) Van Tulder, G.: elasticdeform: elastic deformations for n-dimensional images (2021). URL https://doi. org/10.5281/zenodo 4569691 (2022)
  * (35) Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 (2014)
  * (36) Deforet, M.: Long-range alteration of the physical environment mediates cooperation between pseudomonas aeruginosa swarming colonies. Environmental Microbiology 25(8), 1451–1464 (2023)

## Online Methods

###  3.1 Model Architecture

Encoder module is a sequence of downsampling operations, convolutions, and
residual convolution operations (RConv). A downsampling operation is composed
of 2x2 max-pooling operations concatenated with a 2x2-strided convolution.
RConv operations are composed of two successive convolutions and the second
activation is applied on the sum of the input of the block and the output of
the second convolution. For all convolutions, the activation function is ReLU.
Decoder module is composed of transpose convolutions for upsampling, RConv
operations and convolutions. Encoder and decoder operations are detailed in
Table S1.

Several modules are used to blend the encoded features, which organisation is
depicted in Figure 3. The Combine operation concatenates the inputs and
applies a convolution with a kernel size 1. Frame pairs are generated by the
Pair Blender module that is a Combine operation, with a convolution kernel
size of 5 and the same number of filters as the encoded features. Feature
blending module is depicted in Figure S6. It consists of three Combine
operations applied on encoded features and frame pairs, followed by three
RConv operations. The number of filters of all the operation of the Feature
blending module depends on the number of considered frames, $N$, and is
calculated using the following formula: $(F\times N)/2$, $F$ being the number
of filters of the encoded features ($F=128$ in this work).

###  3.2 Loss functions

Link multiplicity was trained using the categorical cross-entropy loss
function with balanced class frequencies. For EDM, GDCM, and displacement, the
Pseudo-Huber loss function —a smooth approximation of the Huber loss— was used
with a delta parameter of 1 charbonnier1997deterministic ;
friedman2009elements . It combines the advantages of the L2 squared loss and
L1 absolute loss by being strongly convex near the target but not sensitive to
large values. This results in gradients of comparable magnitude for the three
predicted features, regardless of the displacement amplitude, or object size
and shape. For all predicted features except EDM, losses were computed only
within cells.

###  3.3 Data augmentation

We applied several random intensity and geometrical transformations detailed
below, as well as random subsampling to simulate larger displacements. The DNN
time window was drawn in $[7,\;67]$ for the PhC-C2DH-PA14 dataset ($N=7$,
$\delta$ in $[1,\;16]$), and $[7,\;27]$ for the Fluo-C2DH-HBEC dataset ($N=7$,
$\delta$ in $[1,\;6]$). Data augmentation was computed on-the-fly during
training.

####  3.3.1 Intensity transformations

We randomly scaled the intensity in order to improve generalization and
robustness to the presence of long tails in intensity distribution, which is
very common in phase-contrast and fluorescence images. Per-image normalization
to range $[0,\;1]$ was done using random values in the range of percentiles
$[0.01,\;20]$ for the lower bound and $[80,\;99.9]$ for the upper bound. In
order to improve robustness to loss-of-focus, we applied a Gaussian blur with
a random standard-deviation in $[1,2]$ pixels for 50% of the images. In order
to simulate variations in signal to noise ratio, we added Gaussian noise with
a random level. We also performed a random elastic deformation of the
histogram as well as in o2022delta .

####  3.3.2 Geometrical transformations

Images were randomly cropped in tiles, with a +/- 10% random zoom applied
simultaneously with cropping. Tile size was set to $192x192$ for PhC-C2DH-PA14
and $384x384$ for Fluo-C2DH-HBEC. Cropping coordinates were also randomly
shifted between frames to simulate stage drift/instability and increase
displacement variability. Shift was drawn in $[0,\;10]$ pixels for PhC-C2DH-
PA14 and $[0,\;20]$ pixels for Fluo-C2DH-HBEC dataset. We also applied random
horizontal and vertical flip, $90^{\circ}$ rotation as well as elastic
deformation van2022elasticdeform  as in ronneberger2015u , except that we set
a null deformation at edges in order to avoid edge artifacts.

###  3.4 Datasets

####  3.4.1 Fluo-C2DH-HBEC

This dataset consists of fluorescence images of actin-GFP labelled HBEC cells
migrating and dividing on a glass slide, with a time interval of 5 minutes,
and a spatial resolution of 0.908 µm/pixel. Training set consists of one video
of 366 frames of 900x700 pixels containing 18,344 cell observations.
Evaluation set consists of two videos, the first one with $790\times 590$
pixels, 99 frames and 2,568 cell observations, and the second one with
$900\times 700$ pixels, 35 frames and 2,461 cell observations.

####  3.4.2 PhC-C2DH-PA14

This bacterial dataset is made of phase contrast images of Pseudomonas
aeruginosa (strain PA14) cells, arranged as a monolayer on an agar gel, with a
time interval of 10 ms, and a spatial resolution of 0.088 µm/pixel, acquired
at a frame rate of 100 frames-per-second. At this timescale, cell division can
be neglected (division time is typically one hour). The dataset is composed of
a training set of two 200-frames videos, of $700\times 950$ pixels and
$875\times 435$ pixels, containing in total about 920,000 bacteria
observations. The evaluation set is a 200-frame video of $574\times 515$
pixels. This dataset contains many cells that are partially included in the
video, which results in many irrelevant errors. To reduce those errors,
predictions were made on the whole video, but evaluated on a $478\times 364$
pixels subset.

###  3.5 Training

Models were trained for 1000 epochs of 200 steps using the Adam optimizer
kingma2014adam  with learning rate $2.10^{-4}$. Learning rate was decreased by
a factor 2 on plateau of 80 epochs. Mini-batch size was set to 24 for
PhC-C2DH-PA14 and 16 for Fluo-C2DH-HBEC.

###  3.6 Algorithmic details on segmentation

GDCM is the proxy for center segmentation. Geodesic distance is a measure of
distance that takes cell shape into account, it corresponds to the shortest
path within the mask of the segmented cell. Ground truth GDCM was computed
using fast marching method from scikit-fmm.

Centers are segmented by performing a watershed algorithm on the Laplacian
transform of scale $\sigma$ of the image that results from the Gaussian
function ($d\rightarrow e^{-\frac{d^{2}}{2\sigma^{2}}}/\sqrt{2\pi}\sigma$)
applied to the predicted GDCM. $\sigma$ is a constant of typically 25% of the
object thickness and is limited in the range $[1,\;3]$ pixels. Segmented
centers with aberrant size (outside $[50\%,\;200\%]$ of the expected size) or
eccentricity $>0.9$ are filtered out.

###  3.7 Experimental details

####  3.7.1 PA14 growth

Bacterial cells were grown overnight in LB medium at 37°C with aeration. Agar
plates of swarming medium (47 mM Na2HPO 4, 22 mM KH 2PO 4, 8.5 mM NaCl, 1 mM
MgSO 4, 0.1 mM CaCl 2, 5 g/L casamino acids (Bacto, BD)) were obtained via
addition of agar to reach a 0.5% mass fraction deforet2023long . Overnight
bacterial suspension was washed twice in PBS buffer and diluted 1000-fold. 2
µL of the washed suspension were added to inoculate an agar plate, which was
then flipped and placed inside a 37°C microbiological incubator overnight.

####  3.7.2 PA14 observation

Swarming plates were placed inside an stage-top 37°C incubation chamber
(Okolab) mounted on an IX-81 Olympus inverted microscope. Phase contrast
videos of 10 seconds at 100 frames per seconds at 40x magnification were
acquired using a Blackfly S camera (FLIR), with a resolution of 0.088
µm/pixel. The field-of-view was cropped from $2448\times 2048$ to $1000\times
1000$ pixels to reach a frame rate of 100 frames-per-second. The data for
training and evaluation sets were obtained with PA14 wild-type. The data for
validation (Figure 5A-B) was obtained with strain PA14 $\Delta$PilA, which is
flagellated but non-piliated (gift from Dominique Limoli’s laboratory in
University of Iowa). This strain yields more stereotypical long-term diffusive
behavior in monolayers than the wild-type.

####  3.7.3 HBEC growth

Human bronchial epithelial cells (HBECs) (gift from John Minna’s laboratory in
Dallas, TX) were cultured in supplemented keratinocyte serum-free medium with
L-glutamine (Keratinocyte-SFM with L-glutamine; Gibco). Cells were maintained
at 37°C under 5% CO2 partial pressure and 95% relative humidity.

####  3.7.4 HBEC observation

HBECs were seeded at low density in wells of a glass bottom 12-well plate in
which adequate substrates (glass slide for training and evaluation dataset,
PDMS for validation sets) were previously deposited (approximately 26,200
cells/cm2 i.e., 100 000 cells per well). They were left to incubate for a few
hours until they were fully attached. Prior to imaging, cell monolayers were
rinsed with PBS and 5 mL of fresh medium was added. Time-lapse multi-field
experiments were performed on automated inverted microscopes equipped with
thermal and CO2 regulations. Experiments lasted typically several days. The
displacements of the sample, the illumination sequences and the camera
acquisitions were computer-controlled by Metamorph software (Universal
Imaging). Experiments were performed with a 10X objective in fluorescence
microscopy. The images of the same field of view were taken every 5 min. We
used a Leica DM IRB with camera R600 (Qimaging, field of view of $1344\times
1100$ pixels, with a resolution of 0.908 µm/pixel) for training and evaluation
sets, and a Olympus IX71 with camera Prime BSI (Teledyne Photometrics, field
of view of $1024\times 1024$ pixels, with a resolution of 1.3 µm/pixel) for
validation (Figure 5C-D).

###  3.8 Computation

All computations (training, evaluation, validation) were performed on
computers equipped with an NVIDIA RTX A6000 48GB graphics card, 128GB RAM,
Intel Xeon W-2255 running under Microsoft Windows 11, or AMD Ryzen 7 5800X
running under Ubuntu Linux 22.04. Prediction times given in Table 1 were
calculated on the latter. Data analysis for Figure 5 was done with MATLAB (The
Mathworks) from the CSV files extracted from BACMMAN, which contain all
relevant cellular data (track ID, frame ID, XY coordinates, major axis length,
major axis orientation).

## Supplementary tables and figures

Supplementary Table S1: Details of operations by layer in encoder and decoder modules. As described in the main text, Down stands for downsampling, RConv for residual convolution, Conv for 2D convolution, Up for upsampling and BN for batch-normalization. 2D Spatial dropout with 0.2 proportion is added to underlined operations. Numbers in parentheses after the operations are kernel size, number of output filters, and in case of RConv, dilation rate applied to the second convolution. Note that for the decoder associated to displacement, there are two distinct upsampling operations for each axis at layer 1.  Model | Layer |  Encoder |  Decoder
---|---|---|---
D2 | 1 |  Down(3x3, 32) |  Up(4x4, 1)
D2 | 2 |  Conv(3x3, 32), Conv(5x5, 32), Down(3x3, 128) |  Up(4x4, 32), ResConv(3x3, 32), ResConv(3x3, 32)
D2 | 3 |  RConv(5x5, 128, 2), RConv(5x5, 128, 3), RConv(5x5, 128, 4), RConv(5x5, 128, 2), RConv(5x5, 128, 3), Conv(5x5, 128), BN |  Conv(3x3, 64), RConv(3x3, 64), Conv(3x3, 64), BN
D3 | 1 |  Down(3x3, 32) |  Up(4x4, 1)
D3 | 2 |  Conv(5x5, 32), Down(3x3, 64) |  Up(4x4, 32), RConv(3x3, 32), RConv(3x3, 32)
D3 | 3 |  RConv(3x3, 64), RConv(3x3, 64), Down(3x3, 192) |  Up(4x4, 64), RConv(3x3, 64), RConv(3x3, 64)
D3 | 4 |  RConv(5x5, 128, 2), RConv(5x5, 128, 3), RConv(5x5, 128, 4), RConv(5x5, 128, 2), RConv(5x5, 128, 3), Conv(5x5, 128), BN |  Conv(3x3, 64), RConv(3x3, 64), Conv(3x3, 64), BN
| Low density | High density
---|---|---
Surface fraction | 0.446 | 0.719
Average number of cells per image | 2911 | 3620
Number of tracks | 7788 | 3758
Number of tracks analyzed for Figure 5A-B ($>10$ frames) | 7620 | 3675
Number of complete tracks (1000 frames) | 1413 | 3404
Number of incomplete tracks ($<1000$ frames) | 6375 | 354
Number of track ends* | 109 | 5
Supplementary Table S2: Statistics of track duration. *The last line is the
number of track ends that are neither at the last frame nor within a distance
of 3 $\mu m$ from the image boundary. These track ends are biologically
irrelevant and are therefore considered segmentation errors. Notably, there
are only five such cell segmentation errors in the high-density video, out of
approximately 1000 frames $\times$ 3620 cells $=3.6\cdot 10^{6}$ cell
instances. Manual review confirmed these errors are temporary under-
segmentations (case B of Figure 2), involving two pairs of cells.
Supplementary Figure S1: Diagram of the subsampling procedure, presenting the
definition of number of considered frames ($N$), gap between considered frames
($\delta$), DNN time window ($(N-3)\delta+3$). Frames immediately before and
after the central frame are always included as they are used to compute the
displacements. Supplementary Figure S2: Prediction results on Fluo-C2DH-HBEC.
A: input and output of the model. EDM is the map of the Euclidean distance to
the edge of each cell, GDCM is the map of the geodesic distance to the center,
dX and dY are the cell center displacement for each axis, above: from $t$ to
$t+1$ and below from $t+1$ to $t$. No link probability is the probability that
a cell has linked cell in the other frame and multiple link probability is the
probability that a cell has several linked cells in the other frame. For a
better readability, contours of segmented cells are represented in the 8
rightmost images. Note that only one frame pair ($t$, $t+1$) is represented,
but the model inputs a larger temporal neighborhood and predicts these maps
for more frame pairs (see main text). B: Tracking result: contours of cells
are represented by a dashed line and contours of assigned previous cells are
represented by a solid line. Arrows represent the predicted displacement from
$t$ to $t+1$. Supplementary Figure S3: Diagram presenting the forward link
multiplicity and backward link multiplicity, defined for cells at time $t$.
Supplementary Figure S4: Diagram presenting forward and backward tracking. The
first case is a single link between two objects. The second case is either an
over-segmentation at time $t$, an under-segmentation at time $t+1$, or a
predicted merging event. The third case is a predicted split event. The fourth
case is either an over-segmentation at time $t+1$, or an under-segmentation at
time $t$. Supplementary Figure S5: Cells exhibit a brighter region at the
middle of the body, where the invagination of the cell membrane precedes the
complete separation of the two daughter cells (depicted by the magenta
arrowheads). The contact between two separate cells look similar (cyan
arrowhead). This confusion can only be resolved when the two separate cells
actually move apart (yellow arrowhead). Supplementary Figure S6: Blender
module architecture. Encoded features (green, blue and red cubes) and feature
pairs (turquoise and magenta cubes) are respectively combined. Combined
feature and feature pairs are then combined together. Three residual
convolutions are then applied to combined features. For clarity only three
frames have been represented but we considered seven frames in this work.
Supplementary Figure S7: Diagram for the definition of the $2N-4$ frame pairs
used as input of the DNN. Supplementary Figure S8: Adaptation of BACMMAN’s GUI
for hyper-stack visualization. Left: example of hyper-stack visualization. On
can select objects directly on the image and their tracks are displayed as
coloured contours. Left panel displays a part of the GUI of BACMMAN, that
enables manual curation of segmentation and tracking, as well as visualization
and navigation of sub-populations of interest (selections), displayed here in
orange and green. Those selections can be defined from an external statistical
processing software such as Python, R or Matlab. Supplementary Figure S9: CTC
vs DistNet2D metrics. Three examples of segmentation errors (corresponding to
examples shown in Figure 2) and how these errors are counted in CTC and in
DistNet2D. The first column represents the ground truth. The second column
represents the incorrect prediction. The third column represents the error
counting according to matula2015cell . The last column represents the error
counting used in this work. Note that the three examples do not depict all
possible error types, but instead highlight how some segmentation errors are
also counted as tracking errors in CTC but not in DistNet2D. In the first case
(under-segmentation error), we do not count any tracking error: original edges
are not deleted at split operation, but modified by minimizing distances. In
the second case (over-segmentation error), we do not count any tracking error
if both over-segmented cells are linked to a unique cell at previous and next
frames. In the third case (segmentation errors involving two lineages), CTC
allows a reference object to match with only one single predicted object.
Under-segmentation is thus not detected. On the contrary, DistNet2D detects
here one under-segmentation, one over-segmentation, and no tracking error.

## Supplementary videos

Supplementary Video S1: Training video #1 of dataset PhC-C2DH-PA14. To reduce
the total video size, the image stack was time-subsampled 2-fold, but the
video still plays in real time (50 fps instead of 100 fps). Supplementary
Video S2: Training video #2 of dataset PhC-C2DH-PA14. To reduce the total
video size, the image stack was time-subsampled 2-fold, but the video still
plays in real time (50 fps instead of 100 fps). Supplementary Video S3:
Evaluation video of dataset PhC-C2DH-PA14. To reduce the total video size, the
image stack was time-subsampled 2-fold, but the video still plays in real time
(50 fps instead of 100 fps). Supplementary Video S4: Training video of dataset
Fluo-C2DH-HBEC. Supplementary Video S5: Evaluation video #1 of dataset
Fluo-C2DH-HBEC. Supplementary Video S6: Evaluation video #2 of dataset
Fluo-C2DH-HBEC. Supplementary Video S7: Video of the low density monolayer of
PA14 $\Delta$PilA cells, used for Figure 5A-B. To reduce the total video size,
only the first second of the 10-second dataset was included and the image
stack was time-subsampled 2-fold. The video still plays in real time, at 50
fps instead of 100 fps. Supplementary Video S8: Video of the high density
monolayer of PA14 $\Delta$PilA cells, used for Figure 5A-B. To reduce the
total video size, only the first second of the 10-second dataset was included
and the image stack was time-subsampled 2-fold. The video still plays in real
time, at 50 fps instead of 100 fps. Supplementary Video S9: Video of the HBEC
cells, used for Figure 5C-D. To reduce the total video size, the image stack
was time-subsampled 4-fold. The entire video lasts 75 hours.


# Introduction of Machine Learning for Astronomy (Hands-on Workshop)

Yu Wang  ICRA, Dip. di Fisica, Università di Roma “La Sapienza”, Piazzale Aldo
Moro 5, I-00185 Roma, Italy  ICRANet, Piazza della Repubblica 10, I-65122
Pescara, Italy  INAF – Osservatorio Astronomico d’Abruzzo,Via M. Maggini snc,
I-64100, Teramo, Italy.  Rahim Moradi  ICRA, Dip. di Fisica, Università di
Roma “La Sapienza”, Piazzale Aldo Moro 5, I-00185 Roma, Italy.  ICRANet,
Piazza della Repubblica 10, I-65122 Pescara, Italy  INAF – Osservatorio
Astronomico d’Abruzzo,Via M. Maggini snc, I-64100, Teramo, Italy  Mohammad H.
Zhoolideh Haghighi  Department of Physics, K.N. Toosi University of
Technology, Tehran P.O. Box 15875-4416, Iran  School of Astronomy, Institute
for Research in Fundamental Sciences (IPM), Tehran, 19395–5746, Iran  Fatemeh
Rastegarnia  ICRANet, Piazza della Repubblica 10, I-65122 Pescara, Italy  Dip.
di Fisica e Scienze della Terra, Università degli Studi di Ferrara, Via
Saragat 1, I–44122 Ferrara, Italy  Department of Physics,Faculty of Physics
and Chemistry, Alzahra University,Tehran, Iran<EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This article is based on the tutorial we gave at the hands-on workshop of the
ICRANet-ISFAHAN Astronomy
Meeting111https://indico.icranet.org/event/2/page/12-data-science-in-
relativistic-astrophysics-hands-on-workshop. We first introduce the basic
theory of machine learning and sort out the whole process of training a neural
network. We then demonstrate this process with an example of inferring
redshifts from SDSS spectra. To emphasize that machine learning for astronomy
is easy to get started, we demonstrate that the most basic CNN network can be
used to obtain high accuracy, we also show that with simple modifications, the
network can be converted for classification problems and also to processing
gravitational wave data.

##  1 Introduction

Machine learning is no more a new concept as currently it is involved in many
aspects of life. Every time you pick up the phone and unlock it using face
recognition and translate your voice into text when chatting, all have
contributions from machine learning and in terms of the professional field,
the first time that machine learning made a major contribution was in 2016,
when Google’s AlphaGo playing Go defeated the world champion LiShishi and Ke
Jie (Silver et al., 2016, 2017, 2018; Schrittwieser et al., 2020), Alphago
played some of the moves that confused human beings and these are commonly
used now in human-human combat. If you often play go, you will find that many
human players play with a machine style. In other words, humans accept the
changes brought by machines. Training alphaGo used convolutional neural
networks (Fukushima & Miyake, 1982; LeCun et al., 1989). Today our examples
use a similar network. Last year Google’s Deepmind team announced that protein
folding prediction could also be successfully solved in this way and this had
posed a major problem in biology for almost 50 years. (Jumper et al., 2021).
Also they immediately used the trained machine to predict the protein
structure of the COVID19 virus to help in the development of drugs (Arora &
Bist, 2020; Jumper et al., 2021).

Figure 1: Neural network composed of various neurons. features and labels are
the input data represented by grey dots, the hidden layers constructed by
neurons and activation functions are represented by blue dots.

Returning to astronomy, its data records spatial and temporal information,
which is very close to image recognition and speech recognition. The pictures
of galaxies seen by telescopes, the spectra of detected radiation, and the
structure of the universe generated by simulations, all these can be
categorized into image analysis. The evolution of celestial bodies, the
lightcurves of explosions, and simulations of changes at each time step, can
be used with techniques similar to voice recognition. Broadly speaking, the
only axes of the four-dimensional universe we live in are time and space, so
the data is no exception with time and space. What machine learning has to do
is to receive this data and then give an answer.

Abstracted to mathematics, machine learning can be seen as a map, mapping data
to answer. For deep learning (see e.g. LeCun et al., 2015; Goodfellow et al.,
2016; Dong et al., 2021, and references therein), this map consists of many
layers of neurons. Deep learning is a subset of machine learning, which has
developed rapidly in the last decade. Machine learning is also a subset of
artificial intelligence (AI), which has been talked about for many years. The
AI that surpasses human beings in science fiction movies is defined as strong
AI (Strauß, 2018; Butz, 2021), the development of strong AI has only reached
the stage of thesis. But according to our estimation, now the computer can
probably simulate millions of neurons. The human brain has hundreds of
millions of neurons. According to the speed of computer development, when in
about 20 years, the computer can reach the complexity of the human brain, then
strong AI will be really practical.

This article is based on our hands-on tutorial in the workshop. In section 2
we introduce the principle and implementation of neural networks, in section 3
we illustrate it with an example of inferring redshift by SDSS spectra, and in
section 4 we demonstrate that the application can be extended to
classification and other sources by simple changes to the neural network. This
article does not cover the code explanation, running and troubleshooting in
the workshop.

##  2 Basics

Figure 2: Demonstration of the gradient descent process, the loss decreases
following the curve line from the red region to the blue valley via the change
of $\omega$ at each training step, eventually, the $\omega$ approaches
$\omega^{*}$. This figure is reproduced from (Amini et al., 2018).

Figure 1 shows an example of a neural network. The two leftmost grey dots
represent the input, which is passed into the middle neural network, composed
of three layers, and each layer contains four neurons. If we know the correct
answer, we can define a loss function to compare the difference between the
machine’s answer and the real answer. For example, if we input the LIGO data,
we can calculate whether it contains the gravitational wave signal, and the
parameters of the binary system.

Our input can generally be represented as a vector or matrix, corresponding to
each of our neurons (each blue dot in figure 1) doing a matrix operation, as
shown in this equation,

$\vec{y}=A(\omega\cdot\vec{x}+b)$ (1)

where $\vec{x}$ is the input, multiplied by weight $\omega$, which is a matrix
and added by an offset $b$, then the nonlinear function $A$ acts on the
result. When we train the machine, we are actually fitting this $\omega$.
There are thousands of parameters in each $\omega$, and the whole network has
hundreds of neurons, so the total number of parameters reaches the order of a
million. To fit such a large number of parameters, a large amount of data
input is required. The more data input, the more accurate the fit, the better
the prediction ability of the machine. Usually we say the quality of data
defines the upper limit of the accuracy, and the algorithm of the network
always tries to improve the accuracy to reach the upper limit.

Considering that we have a network and many data, how can we train the
machine? Here we take supervised learning as an example, which is the most
common training method at this stage. Supervised learning requires that the
data provided has the correct answer. It is like teaching a child to
distinguish a dog from a cat, you show him a cat picture, and he answers dog,
then you correct him “you are wrong, it should be a cat”. After many times of
repeated training, the child masters the method of distinguishing between cats
and dogs. Mathematically expressed, it is to find a set of parameters $\omega$
that makes the neural network output the same answer as the real answer, that
is, the loss function is close to 0. Figure 2 represents the loss
corresponding to a certain $\omega$. The initial $\omega$ corresponds to a
very high loss in the red region. For each input of data, or each training, we
seek the gradient of the loss, our $\omega$ moves along the gradient by one
small step in the inverse direction, with the step size controlled by $\eta$.
After many training sessions, $\omega$ step by step eventually goes to a very
small loss, the blue valley on the figure.

Figure 3: Normally, a complete training process includes initiation, forward
propagation, loss computation, backpropagation and the final convergence. The
left side shows some examples of each step.

Figure 4: The spectra of quasars. There are several emission/absorption lines
present in the spectra, including Ly$\alpha$ (121.6 nm), $\rm CIV$ (154.9 nm),
$\rm CIII$ (190.9 nm), $\rm MgII$ (279.6 nm), $\rm H\beta$ (486.2 nm) and $\rm
H\alpha$ (656.3 nm) as well as a $\rm CIV$ line with a broad absorption
feature. The spectral number, visually inspected redshift and the signal to
ratio are labelled on each figure, respectively. This figure is reproduced
from Rastegarnia et al. (2022).

We go over the procedure, as shown in figure 3. First we design a neural
network and initialize the parameters. Then we do the forward propagation to
make a prediction. By comparing with the ground truth, we obtain the loss.
Then the error is back propagated to the entire network to refine the
parameters. We repeat this training process by feeding a lot of data until
convergence is reached.

##  3 Example

Figure 5: Structure of one dimensional CNN. The quasar spectrum is input as a
one-dimensional array, which goes through the convolutional layer of kernel
size = 500, 200, 10 respectively to search for the global and local pattern.
The fully connected layers output the redshift.

Our example follows Rastegarnia et al. (2022), that is to infer the
spectroscopic redshift of quasars by deep learning. The relevant code and data
can be found in https://github.com/YWangScience/Isfahan-workshop-2021/ and
https://www.kaggle.com/datasets/ywangscience/sdss-iii-iv respectively.

The quasar spectra are retrieved from Sloan Digital Sky Survey (SDSS) Data
Release 16 and quasars only (DR16Q) (Hutchinson et al., 2016; Lyke et al.,
2020), which is one of the best samples to train the neural networks, with
more than $700,000$ quasars, and includes $326,535$ quasars that are visually
inspected. The spectra are standardized by fitting and then extrapolating to
4618 data points (pixels) uniformly distributed in $\log\lambda$, where
$\lambda$ is the wavelength in the range of $360~{}\rm nm-1032.5~{}\rm nm$.
Then the spectra are normalized using the Zero-Mean Normalization method
(Jayalakshmi & Santhakumaran, 2011) as features. The corresponding redshifts
are exported as labels. We take $90\%$ of samples for training and $10\%$ for
testing. Figure 4 shows two spectra with the emission/absorption lines marked.

In the field of machine learning, the most commonly used network architecture
is the convolutional neural network (CNN), and many complex networks are built
on top of the CNN. For our simple 1D data, a very deep network is rarely
needed, and our experiments also demonstrate that building a CNN network with
about $10$ layers is sufficient to obtain accurate redshift. We also tested
that increasing the complexity of the network by a factor of 10 or more, or
using a state-of-the-art network, yields an accuracy improvement of only $\sim
1\%$, so in this tutorial article we use the most basic CNN.

When designing CNN networks, we need to start from the actual physical
problem. For the problem of deriving the redshift from the spectrum, we learn
that the combination of these patterns can estimate the redshift:

  1. 1.

The global shift of spectrum;

  2. 2.

The shift of the emission and absorption lines at different redshifts;

  3. 3.

Some specific signals may appear at given redshifts.

Hence, in the convolutional part, we specially construct a large size filter
of $500$ pixels covering $>10\%$ of the data to capture the global shift of
the spectrum, and in series a middle size $200$ pixels and a small size filter
of $10$ pixels to capture the small and minor shift and those specific
signals. The network architecture is shown in figure 5.

The code of the whole process is written according to to figure 3,
corresponding to:

  1. 0.

He uniform (He et al., 2015) is adopted to initialize the network weights;

  2. 1.

forward propagation follows our net defined in figure 5;

  3. 2.

we adopt the mean squared error as the loss

$MSE(x)=E((x-x^{*})^{2})$ (2)

where $x$ is the prediction and $x^{*}$ presents the ground truth.

  4. 3.

Adam is selected to perform the backpropagation;

  5. 4.

we consider convergence when training cannot reduce the loss.

After dozens of epochs of training, this simple network obtains an accuracy of
$>97\%$ for the redshift of the testing samples.

##  4 Extend the network by simple changes

Figure 6: Example of the gravitational wave data, the blue curve shows the
original data, the orange curve is the window for Fourier transformation, and
the green curve is the filtered signal that only contains $20-500$ Hz data.

Here by changing some lines of the net, we turn to train the net to classify
the SDSS objects, and for simplicity of the tutorial, only involve two classes
of quasar and galaxy.

This classification example uses the same input as the redshift example, a
series of one-dimensional data. The outputs are different, The redshift
example finally gives a redshift value, while here we output two values,
corresponding to the probability of each class. In order to standardize the
output, a LogSoftmax function

$p(x_{i})=\log\left(\frac{\exp(x_{i})}{\sum_{j}\exp(x_{j})}\right),$ (3)

is applied at the last layer. Hence all values are between 0 and 1, and the
sum of all equals one. For instance, the output ($0.2,0.8$) indicates that the
source belongs to the second class. The loss function is also different, and
the resulting $q(x_{i})$ (predicted classification) will be adopted and
together with the labels p(xi) (true ground classification) to compute the
cross entropy as the loss

$H(p,q)=\sum_{i}p(x_{i})\log q(x_{i})$ (4)

With the above two simple changes to training the net following the same
procedure of our first redshift example, this simple net is capable of making
the classification of more than $80\%$ accuracy for quasar and galaxy.

With some simple modifications, this net can also be used for gravitational
wave detection. Similar to the SDSS spectrum, the gravitational wave signal
forms a one-dimensional data structure, that the time versus the strain, or
the frequency versus the amplitude if changing from the time domain to the
frequency domain, see figure 6. The output will be a value indicating the
existence of gravitational wave. So we hardly need to modify the network
structure, just to modify the length of some parameters according to the
length of the input data. Our test using the data provided from Kaggle
competition (https://www.kaggle.com/c/g2net-gravitational-wave-detection/)
shows such a simple net obtains an accuracy of more than $>80\%$, only about
$\sim 5\%$ less accurate than the top one network. We also noticed in this
competition that the top three networks did not use very complex networks but
put much thought into pre-processing the data.

##  5 Conclusion

Analyzing and processing data is the backbone of astronomy, and machine
learning has arguably advanced the most. As a relatively new concept for
astronomy, machine learning may be difficult for astronomers to start with. It
occurred to us to illustrate through this article that, thanks to the easy-to-
use programming frameworks developed by the industry, machine learning in
astronomy does not require very difficult programming theories and
capabilities. In fact, the main development of machine learning in the
industry is in the processing of images and voices, which are actually
recorded in a form no different from astronomical data, both in temporal and
spatial data formats. These similarities allow us to apply many of the neural
networks matured in the industry to astronomical data processing. In
conclusion, machine learning as a tool is easy to use, but astronomers and
astrophysicists need brilliant ideas to use this tool flexibly.

## References

  * Amini et al. (2018) Amini, A., Soleimany, A., Karaman, S., & Rus, D. 2018, arXiv preprint arXiv:1805.04829
  * Arora & Bist (2020) Arora, K., & Bist, A. S. 2020, Aptisi Transactions on Technopreneurship (ATT), 2, 120
  * Butz (2021) Butz, M. V. 2021, KI-Künstliche Intelligenz, 35, 91
  * Dong et al. (2021) Dong, S., Wang, P., & Abbas, K. 2021, Computer Science Review, 40, 100379
  * Fukushima & Miyake (1982) Fukushima, K., & Miyake, S. 1982, in Competition and cooperation in neural nets (Springer), 267–285
  * Goodfellow et al. (2016) Goodfellow, I., Bengio, Y., & Courville, A. 2016, Deep learning (MIT press)
  * He et al. (2015) He, K., Zhang, X., Ren, S., & Sun, J. 2015, in Proceedings of the IEEE international conference on computer vision, 1026–1034
  * Hutchinson et al. (2016) Hutchinson, T. A., Bolton, A. S., Dawson, K. S., et al. 2016, AJ, 152, 205, doi: 10.3847/0004-6256/152/6/205
  * Jayalakshmi & Santhakumaran (2011) Jayalakshmi, T., & Santhakumaran, A. 2011, International Journal of Computer Theory and Engineering, 3, 1793
  * Jumper et al. (2021) Jumper, J., Evans, R., Pritzel, A., et al. 2021, Nature, 596, 583, doi: 10.1038/s41586-021-03819-2
  * Jumper et al. (2021) Jumper, J., Evans, R., Pritzel, A., et al. 2021, Proteins: Structure, Function, and Bioinformatics, 89, 1711
  * LeCun et al. (2015) LeCun, Y., Bengio, Y., & Hinton, G. 2015, nature, 521, 436
  * LeCun et al. (1989) LeCun, Y., Boser, B., Denker, J. S., et al. 1989, Neural computation, 1, 541
  * Lyke et al. (2020) Lyke, B. W., Higley, A. N., McLane, J. N., et al. 2020, ApJS, 250, 8, doi: 10.3847/1538-4365/aba623
  * Rastegarnia et al. (2022) Rastegarnia, F., Mirtorabi, M. T., Moradi, R., Vafaei Sadr, A., & Wang, Y. 2022, MNRAS, 511, 4490, doi: 10.1093/mnras/stac076
  * Schrittwieser et al. (2020) Schrittwieser, J., Antonoglou, I., Hubert, T., et al. 2020, Nature, 588, 604, doi: 10.1038/s41586-020-03051-4
  * Silver et al. (2016) Silver, D., Huang, A., Maddison, C. J., et al. 2016, Nature, 529, 484, doi: 10.1038/nature16961
  * Silver et al. (2017) Silver, D., Schrittwieser, J., Simonyan, K., et al. 2017, Nature, 550, 354, doi: 10.1038/nature24270
  * Silver et al. (2018) Silver, D., Hubert, T., Schrittwieser, J., et al. 2018, Science, 362, 1140, doi: 10.1126/science.aar6404
  * Strauß (2018) Strauß, S. 2018, Big Data and Cognitive Computing, 2, 16