# The cell signaling structure function

Layton Aho  Electrical and Computer Engineering, Drexel University, USA  Mark
Winter  Electrical and Computer Engineering, Drexel University, USA  Marc
DeCarlo  Electrical and Computer Engineering, Drexel University, USA  Agne
Frismantiene  Institute of Cell Biology, Univ. of Bern, Switzerland  Yannick
Blum  Institute of Cell Biology, Univ. of Bern, Switzerland  Paolo Armando
Gagliardi  Institute of Cell Biology, Univ. of Bern, Switzerland  Olivier
Pertz  Institute of Cell Biology, Univ. of Bern, Switzerland  Andrew R. Cohen

###### Abstract

Live cell microscopy captures 5-D $(x,y,z,channel,time)$ movies that display
patterns of cellular motion and signaling dynamics. We present here an
approach to finding spatiotemporal patterns of cell signaling dynamics in 5-D
live cell microscopy movies unique in requiring no _a priori_ knowledge of
expected pattern dynamics, and no training data. The proposed cell signaling
structure function (SSF) is a Kolmogorov structure function that optimally
measures cell signaling state as nuclear intensity w.r.t. surrounding
cytoplasm, a significant improvement compared to the current state-of-the-art
cytonuclear ratio. SSF kymographs store at each spatiotemporal cell centroid
the SSF value, or a functional output such as velocity. Patterns of similarity
are identified via the metric normalized compression distance (NCD). The NCD
is a reproducing kernel for a Hilbert space that represents the input SSF
kymographs as points in a low dimensional embedding that optimally captures
the pattern similarity identified by the NCD throughout the space. The only
parameter is the expected cell radii ($\mu m$). A new formulation of the
cluster structure function optimally estimates how meaningful an embedding
from the RKHS representation. Results are presented quantifying the impact of
ERK and AKT signaling between different oncogenic mutations, and by the
relation between ERK signaling and cellular velocity patterns for movies of
2-D monolayers of human breast epithelial (MCF10A) cells, 3-D MCF10A spheroids
under optogenetic manipulation of ERK, and human induced pluripotent stem
cells .

Index Terms Patterns of cell signaling, metric learning, kinase translocation
reporters, ERK, AKT, embedding live cell microscopy movies, cytonuclear ratio

##  I Introduction

Cell signaling in biology refers to the mechanisms with which a cell interacts
with its environment. A key open challenge in systems biology is to elucidate
and manipulate these cellular / environmental interaction mechanisms. For
example, ERK and AKT are both kinases that feature prominently in the control
of basic cellular functions such as motion, cell cycle, _etc._ [1, 2, 3, 4].
In microscopy, a new class of image-based biosensors called kinase
translocation reporters (KTRs) [5, 6] offer improved imaging of multiple
kinases simultaneously in living cells. For both KTRs and existing nuclear
markers like Histone 2B (H2B), we define the signaling state of a single cell
at a fixed time and imaging channel based on the intensity configuration of
its nuclear voxels w.r.t. the surrounding cytoplasmic voxels. We also define
the notion of functional signals, measured as _e.g._ the instantaneous cell
velocity between adjacent image frames in order to quantify the relationships
between cell signaling and functional cellular outputs.

We propose a tool for discovering patterns of cell signaling in microscopy
movies, with the key steps shown in Figure 1. The input is a collection of $N$
5-D microscopy movies, $(x,y,z,channel,time)$. The movies are segmented to
identify individual cells, and optionally tracked and lineaged. The new cell
signaling structure function (SSF) is an optimal and robust measure of cell
signaling state, in contrast to the widely used but mathematically poorly
behaved cytonuclear ratio[5, 3]. The SSF works for KTR markers, and more
generally for any fluorescence microscopy marker. SSF kymographs are 3-D
images that record at each $(x,y,time)$ cell centroid location the cell
signaling state. The kymographs are generated exactly for 2-D movies, for 3-D
movies we project to two spatial dimensions plus time. We would prefer to
maintain the full 3-D spatial information plus time in the kymograph but there
does not yet exist an effective lossless compression algorithm for 4-D images.
Details of the SSF are given in Section I-A.

The normalized compression distance (NCD) is a metric distance based on
Kolmogorov complexity [7] that uses lossless file compression algorithms to
identify patterns of similarity in digital objects. The Free Lossless Image
Format (FLIF) is a lossless 3-D image compression algorithm based on entropy
coding [8], used here with the NCD to find patterns of spatiotemporal
signaling similarity between kymographs with no _a priori_ knowledge or
training data required[9]. We compress each of the $N$ input movies separately
and together. Movies having common patterns compress better together _vs._
separately. Because the NCD is a metric distance, or equivalently, a positive
semi-definite function obeying the (1) identity, (2) symmetry and (3) triangle
inequality, the resulting $NxN$ pairwise distance matrix defines a
_reproducing kernel Hilbert space_ (RKHS) [10]. We call this RKHS an embedding
space, and by preserving the principal $K$ dimensions of the RKHS ($K<N$) we
optimally preserve the characteristics of the input image similarity space
defined by the NCD. Each input movie is represented by a single point in this
real valued $K-dimensional$ space, $\mathbb{R}^{K}$. The RKHS is a metric
space, meaning that the Euclidean distance between points throughout the RKHS
optimally estimates the NCD between equivalent points in the input image
space. Learning in this RKHS is called semi-supervised spectral learning [11,
12]. The pairwise NCD matrix defining the RKHS is generated with no need for
labeled or training data. Data with known labels embedded in the RKHS can then
be used to train supervised learners with reduced complexity in the
$\mathbb{R}^{K}$ space compared to the input image space.

The Cluster Structure Function (CSF) measures exactly the goodness of fit, or
how meaningful, some clustering or partitioning of a dataset[9]. The CSF is
computed from the compression statistics of the input kymographs individually
and combined as clusters. The size and quantity of 3-D kymographs used here
are so large that compressing clusters together is not tractable. Leveraging
the metric characteristics of the RKHS embedding, we propose a new formulation
of the CSF in the RKHS embedding space as an optimal approximation. Details
are in Section V-E. One unique aspect of the present work is that we are not
focused on classifying unknown data, but rather on identifying common cell
signaling patterns among the kymographs and between different imaging channels
in the same kymograph. We use the CSF throughout to measure how much structure
_e.g._ ERK signaling accounts for in the pattern space partitioned according
to the ground truth labels and we use the NCD directly between imaging
channels in the same movie to quantify relationships between _e.g._ ERK
signaling and velocity. Each SSF kymograph defines an individual pattern of
spatiotemporal signaling and cellular motion, the CSF and NCD together measure
common patterns among kymographs from different experimental conditions.

Figure 1: The normalized compression distance is a metric embedding kernel
that finds spatiotemporal patterns of cell signaling in SSF kymographs. Start
with $N$ 2-D/3-D multi-channel time-lapse live cell microscopy movies.
Segmentation identifies the centroid and radius of each cell. The cell
signaling state is quantified by the pixel intensity of the nuclear pixels
w.r.t the surrounding cytoplasm. The SSF kymograph stores this cell signaling
state at each cell centroid. 3-D image compression is then used as a
normalized pairwise distance metric between kymographs. The resulting distance
matrix defines a reproducing kernel Hilbert space (RKHS) an optimal
representation of the input similarity space defined by the compression-based
similarity measure. Each of the 3-D kymographs is represented as a single
point in this RKHS. Finally, the cluster structure function quantifies how
much structure is associated with a given input channel and optionally track
subset. The only parameter is the range of expected cell radii in the images
(in $\mu m$), and no training data or _a priori_ knowledge is required.

Figure 2 shows a sample 2-D image frame from the ERK-KTR channel (A). This
image is taken from a movie of cells with a mutation (PIK3CA_H1047R) that
exhibits a distinctive pulsing throughout the tissue monolayer. ERK-KTR
activation increases as the nucleus darkens w.r.t. surrounding cytoplasm. The
3-D kymograph (B) writes the SSF activation at each $(x,y,time)$ cell
centroid. A 2-D projection of the 3-D kymograph (C) provides an effective
visualization of the ERK pulsing dynamics as diagonal yellow lines. Finally,
the RKHS embedding of the 147 input movies from six different genetic
conditions, each associated with distinctive patterns of ERK signaling
dynamics, is shown in (D). The principal 3 dimensions $(k1,k2,k3)$ are shown.
The RKHS advantage over other embedding methods is that the Euclidean distance
between points is an optimal low-dimensional estimate of the differences in
signaling patterns between the corresponding spatiotemporal signaling dynamics
observed in the input movies. This makes the RKHS especially well-suited for
subsequent supervised (semi-supervised spectral learning [11]) or unsupervised
(spectral clustering [12]) machine learning algorithms.

Figure 2: Embedding SSF kymographs in a metric, reproducing kernel Hilbert
space enables improved visualization and quantification of cell signaling
dynamics from time-lapse microscopy movies. Spectral learning uses the NCD
matrix from pairwise comparison among a collection of SSF kymographs to define
a reproducing kernel Hilbert space (RKHS). Each kymograph is represented as a
point in the low dimensional RKHS that captures optimally the patterns of
visual similarity from the input images. Timelapse microscopy captured 147
movies from 6 experiments showing ERK signaling in 2-D monolayers of human
breast epithelial (MCF10A) cells (A). The cell signaling structure function
measures the voxel intensities across the nucleus w.r.t. the cytoplasm for
each detected cell, and stores that value at the (x,y,t) centroid in a 3-D
kymograph (B). Pairs of 3-D kymograph, as shown in (B) with a single image
frame from (A) rendered on top of the cell centroid signaling-state values,
are input to the FLIF 3-D compression to compute the pairwise NCD matrix. A
maximum intensity projection of the 3-D kymograph (B) to a single spatial
dimension creates the 2-D kymograph (C), allowing easy visualization of the
diagonal yellow stripes that show waves of ERK signaling activation
propagating throughout the tissue. The red line in (C) indicates time of the
image frame in (A,B). In the RKHS (D) each marker represents a single input
kymograph. The movies / kymographs are grouped clearly by the different
mutations that are associated with distinct changes in signaling patterns (D).
There are 6 different classes corresponding to different mutations, so the
RKHS requires a 6-D embedding for classification and for the cluster structure
function (CSF) reported in (D) as [mean, standard deviation] of the per-
cluster optimality deficiency for ERK signaling in the 6-D RKHS. See also
supplementary movie 1.

The introduction continues below with a description of the cell signaling
structure function, and a very brief review of the related literature. Section
II details the results of our analysis for 2-D monolayers of human breast
epithelial (MCF10A) cells from six different oncogenic mutations associated
with distinctive changes in cell signaling patterns, for human induced
pluripotent colonies under self-renewing and differentiating conditions, for
optogenetic excitation of MCF10A cells cultured in 3-D synthetic breast
spheroids, and for a synthetic phantom dataset. Section III gives some
concluding observations on the approach and suggests avenues for future
research. Section V details the methods for live cell imaging, for
segmentation and tracking, and for computing the SSF and the SSF kymographs.

###  I-A The cell signaling structure function

A Kolmogorov structure function measures exactly the _“the goodness-of-fit of
an individual model with respect to individual data”_[13]. The first structure
function, proposed by Kolmogorov in a 1974 talk, modeled any digital object
$x$ using a finite set $S$ containing $x$, and was defined as
$H(x)=\min\limits_{S}(log(|S|),x\in S$. The idea is that if there are no “
_simple special properties_ ”[13] that would allow $x$ to be specified within
$S$ more efficiently than encoding its ordinal within the set then $S$ is an
optimally meaningful representation for $x$. We propose here a new cell
signaling structure function (SSF) that quantifies exactly the signaling state
of each cell at each time and imaging channel. The cell signaling state is
defined by the intensity of the nuclear voxels w.r.t. surrounding cytoplasmic
voxels. The nuclear intensity is defined by $n$, the surrounding cytoplasmic
intensity is defined by $c$, with $c,n\in[0.0,1.0]$. The cell signaling
structure function (SSF) is a 2-D vector valued function. The first dimension
is non-zero when $c>n$, the second dimension is non-zero when $n>c$. We take
as the model $M$ the case $c=n$, an absence of signal, and the SSF is defined
as

$\displaystyle H_{SSF}$ $\displaystyle=<c-n,0>$ $\displaystyle c>n$ (I.1)
$\displaystyle=<0,n-c>$ $\displaystyle c<n$ $\displaystyle=<0,0>$
$\displaystyle c=n.$

The SSF starts at $0$ when the model fits perfectly with no signaling activity
($c=n$) and increases to its maximum value as signaling activity increases to
$c=1,n=0$. This definition differs from both the original Kolmogorov structure
function and the recently proposed cluster structure function[9] in that is
2-dimensional. Another difference between the SSF compared to the Kolmogorov
structure function and the CSF is that the SSF is not intended to optimize
model parameters but rather used directly as a measure of model fit indicating
cell signaling activation levels. Another difference is that the SSF is
normalized to $[0.0,1.0]$ in each dimension. For convenience, we encode
$H_{SSF}$ into a signed scalar value $\in[-1.0,1.0]$,

$H_{SSF}=c-n.$ (I.2)

We use the formulation in I.2 throughout the remainder. Our implementation of
the SSF uses a Laplacian of Gaussian (LoG) blob enhancing filter [14] whose
response scales naturally to $[-1.0,1.0]$, as described in section V-C. Figure
3 shows the SSF in comparison to the cytonuclear ratio currently used to
quantify KTR activation levels. Using a phantom image (A), we simulate KTR
activation (B), with the SSF activation level (black line) precisely matching
the ground truth activation percentage (yellow line). The green line (B)
represents the unnormalized LoG response, and the red line shows the
cytonuclear response, illustrating why the use of the cytonuclear ratio
requires careful human supervision for KTR signal quantification. In (C), the
mapping between cytoplasmic and nuclear intensity is shown for the cytonuclear
ratio, and in (D) for the SSF. Note how the linearity and normalization of the
SSF enhances the color space mapping considerably compared to the cytonuclear
ratio. For the 2-D monolayer movies of human breast epithelial (MCF10A) cells
we compared the CSF results obtained from the SSF on the ERK channel to the
results obtained using the cytonuclear ratio, and found the SSF resulted in a
significantly lower CSF value indicating that the SSF extracted significantly
more structure or meaningful information compared to the cytonuclear ratio
($p<3e-10$), see Section II-A.

######  Theorem 1.

$H_{SSF}$ is a positive semi-definite function, making it an optimal mapping
from the 2-D $(cytoplasmic,nuclear)$ intensity space to a single 1-D
activation value.

###### Proof.

Consider the case $c\geq n$ in Eqn. I.1, given $c\in[0.0,1.0],n\in[0.0,c]$,
$c-n$ is equivalent to the Euclidean distance between the 1-D point $(c-n)$
and the point $(0)$ and the result follows from the Euclidean distance being
positive semi-definite. Combine the case $c<n$ using the sum of two positive
semi-definite functions is positive semi-definite, making $H_{SSF}$ a
reproducing kernel Hilbert space (RKHS) embedding of the $(c,n)$. The
optimality of the mapping from $(c,n)\mapsto H_{SSF}$ follows from the RKHS
properties [10]. ∎

######  Remark 1.

Being a positive semi-definite function is equivalent to being a metric
distance function. The cell signaling structure function can then be
interpreted as a metric distance between the model of no signaling activation
($c=n$) and our given $(c,n)$ configuration. Being positive semi-definite
seems a necessary condition for any structure function, but this has not been
expressed explicitly in previous structure function definitions [9, 13, 15].
$\Diamond$

Figure 3: The cell signaling structure function (SSF) is an optimal measure of
nuclear vs. cytoplasmic voxel intensity that varies linearly with the nuclear
intensity w.r.t. the surrounding cytoplasm. Phantom image (A) shows dark
nucleus against bright background. This represents a cell that is fully
activated w.r.t. a KTR signal. The region between the red and magenta circles
is used to estimate the nuclear intensity, the region between the red and blue
circles is used to estimate the cytoplasmic intensity. Adding 40% shot noise
to the phantom image (B) such that 20% of the nuclear pixels are set to white
and 20% of the cytoplasmic pixels are set to black simulates KTR signal
decrease, representing a cell at 60% activation. Varying the activation levels
(C) shows the SSF accurately captures KTR signal activation across variations
in nuclear and cytoplasmic intensity. The activation measure in (C) is shown
with the ground truth varying linearly from 0.0 to 1.0 as simulated signal is
added to the nucleus and subtracted from the cytoplasm. Error bars in (C)
represent standard error over 100 trials. Unlike the cytonuclear ratio
(CNratio), the SSF varies linearly from [0,1] as kinase translocation reporter
moves from nucleus to cytoplasm. The cytonuclear response (D) varies non-
linearly with nuclear intensity. The cytonuclear ratio goes to infinity as
nuclear intensity goes to zero. The SSF response (E) varies linearly with
nuclear and cytoplasmic intensity and is bounded on $[-1.0,1.0]$.

###  I-B Related Literature

The two key contributions presented here are the the cell signaling structure
function (SSF), and the use of the metric normalized compression distance
(NCD) [7] with the FLIF lossless 3-D image compression to embed SSF kymographs
into a reproducing kernel Hilbert space that preserves the metric
characteristics of the distance measure. The SSF improves on current
approaches for measuring cell signaling activity from kinase translocation
reporters (KTRs). KTRs are fluorescent reporters that utilize
nucleocytoplasmic shuttling to measure kinase activity in a single cell[5, 6].
KTRs have been developed for many signaling pathways, such as ERK (ERK-KTR),
AKT (Fox01-FP), and Cd2k (DBH)[6], a cell cycle indicator. The current
approach for measuring the estimated signal activity from the KTR is to
calculate the ratio of cytoplasmic versus nuclear fluorescence intensity (C/N
ratio) from the reporter [5, 6, 1, 16, 3, 2]. The cytonuclear ratio is non-
linear and unbounded, a poor choice for a measure of cell signaling activation
necessitating careful supervision to select the valid response range . In
contrast, the SSF is an optimal 1-D representation of the 2-D
$(cytoplasmic,nuclear)$ intensity space as shown in theorem 1.

The most widely used alternatives to our approach include neural network
latent spaces and / or non-linear embedding methods like t-Stochastic Neighbor
Embedding (t-SNE) and UMAP. We use the metric normalized compression distance
to generate a pairwise distance matrix between each of $N$ kymographs. The
normalized eigenvectors and eigenvalues of this $NxN$ distance matrix define
an embedded space $\mathbb{R}^{N}$ and we say the NCD is a _reproducing kernel
for a Hilbert space_ (RKHS)[10]. With the RKHS we can use the Euclidean
distance in the embedded space as a valid distance w.r.t. the input space.
This is a key advantage of our method compared to non-linear neural network
latent spaces or non-linear embeddings like t-SNE or UMAP where a distance
metric does not exist. Spatial reasoning and quantitative visualization of the
data in the embedded space, as in our clustering plots, is effective here
because of the metricity of the RKHS. We also use the characteristics of the
RKHS to compute the cluster structure function in the embedded space, only
possible because of the valid inner product (metric distance) provided. The
RKHS combined with the NCD provide a potent tool for pattern discovery in 5-D
images.

We use the NCD to estimate a normalized measure of the _“most concise
description of any and all differences”_[7], here between all pairs of SSF
kymographs, with no required _a priori_ domain knowledge or training data.
This approach is unique in its unsupervised, metric RKHS representation. The
SSF measures cell signaling for each cell in each image frame as the intensity
configuration of nuclear voxels w.r.t. the surrounding cytoplasmic voxels. We
form SSF kymographs, writing the SSF value on a given channel at the
$(x,y,time)$ centroid location for each cell centroid. Other approaches
include CODEX, finding similar patterns of cell signaling from individual cell
trajectories of cytonuclear ratios input to a convolutional neural network
(CNN) and the resulting CNN features are embedded using t-SNE for
visualization and classification[16]. Another approach uses the ARCOS
algorithm, segmenting and tracking collective signaling events by thresholding
cytonuclear intensities against a minimum size of collective events[2]. These
approaches are intended to identify specific types of signaling events and are
useful in combination with the more general signaling pattern discovery
techniques proposed here.

##  II Results

###  II-A ERK and AKT signaling in 2-D+time monolayer of human breast
epithelium

Live movies of human breast epithelial cell monolayers were captured in six
different imaging experiments. Each movie contains ERK-KTR and AKT-KTR
channels, as well as an H2B (nuclear marker). Each experiment captured four or
five movies from each of wildtype cells plus five different mutations for a
total of 147 movies. The movies contained an H2B nuclear marker, as well as
ERK and AKT KTR reporters on three separate channels. The movies were captured
at 5 minutes per frame for 24 hours. Each frame was $1024x1024$ pixels. Figure
2 shows an example image frame from this data (A), the 3-D kymograph (B), the
2-D projected kymograph (C) and the embedding of the 147 movies into a 3-D
RKHS (D). The mutations chosen for this application exhibit distinct patterns
of ERK signaling, as seen in 2C as diagonal yellow stripes of ERK activation
across the monolayer. These distinct signaling patterns are captured by the
FLIF compression as can be seen by the segregation of the representative
points for each movie in the RKHS, 2D. We compute the cluster structure
function (CSF) in the RKHS embedding, preserving the principal six dimensions
of the RKHS embedding as there are six ground truth classes. The CSF is
computed against the ground truth classes, in the RKHS as the distance from
each point to the centroid of the mutation associated with that movie, as in
Section V-E. For ERK, the CSF is $[1.1,0.23]$ (the first value is the mean CSF
across the 147 kymographs, the second value is the standard deviation). For
AKT, the CSF is $[1.3,0.18]$, significantly higher ($p<1e-10$) compared to the
ERK CSF, indicating that AKT is significantly less informative compared to ERK
w.r.t. the observed signaling patterns.

Cellular velocity kymographs were also generated for each movie, writing the
average velocity for each cell from time $t-1$ to $t$ and $t$ to $t+1$ at each
$(x,y,t)$ cell centroid. The velocity is normalized to the maximum cellular
velocity gate value used in the tracking algorithm (V-B). For each of the 147
movies, we compute the NCD between the ERK kymograph and the velocity
kymograph. 4A shows the resulting NCDs the ERK and velocity kymographs as a
multiple comparison test grouped by experimental condition. Lower values of
the NCD indicate that the ERK signal is more predictive of the cellular
velocity pattern. For the PIK3CA_H1047R and PIK3CA_E545K mutations, there is a
significant increase in the relationship between ERK Signaling and cellular
velocity ($p<0.003$).

###  II-B ERK signaling in colonies of human stem cells

Live movies showing colony development of human induced pluripotent stem cells
(iPSCs) were captured. Self-renewing cells are iPSCs that divide to produce
other iPSCs, differentiated cells are formed as the cells progress towards a
neural fate. The movies were labeled with ERK-KTR and H2B (nuclear reporter).
Ten 2-D movies of self-renewing colonies and ten of differentiated colonies
were captured. SSF kymographs were generated for both ERK and the pairwise NCD
matrix was generated. Supplementary figure 1 shows 2-D projections of the ERK
kymographs for the 10 differentiated movies (A) and the self-renewing movies
(B). The CSF for these 20 movies was $[1.1,0.2]$, and the embedding did not
indicate clear separation between ERK signaling patterns for self-renewing
_vs._ differentiated movies. Additionally, there was no statistically
significant change in the relation between ERK and cellular velocity between
the self-renewing and differentiated colonies, as shown in Figure 4B.

Figure 4: Quantifying the relationship between ERK signaling dynamics and
cellular velocity patterns. Multiple comparison plots showing normalized
compression distance between SSF kymographs for MCF10A human breast epithelial
cell movies grouped by mutation (A) and induced human pluripotent stem cells
(hPSCs) grouped by self-renewing vs. differentiated (B). ERK SSF and
normalized velocity are computed for every cell in every movie and written
into an $(x,y,time)$ 3-D kymograph. The normalized compression distance (NCD)
between ERK and velocity kymographs measures how informative the ERK dynamics
are in representing the velocity. Lower NCD values (as in e.g. PIK3CA_H1047R)
indicate that there is a stronger relationship in that condition between ERK
signaling patterns and cellular velocity. Whiskers indicate 95% confidence
intervals from standard error of the mean. 148 movies from 5 experiment dates
were analyzed for the MCF10A cells (A), 10 movies each from self-renewing and
differentiated were analyzed for hPSCs (B). In the MCF10A cells, the
PIK3CA_H1047R and PIK3CA_E545K show significant increases in the relationship
between ERK and velocity compared to wild type. In the hPSC cells, there was
no significant change in the relationship between ERK and velocity between
self-renewing and differentiated cells.

###  II-C Optogenetic excitation of 3-D+time human breast epithelial spheroids

Optogenetic excitation of ERK signaling of 3-D breast spheroids during live
imaging is used to quantify the relationship between ERK signaling patterns
and the resulting cellular velocity. Ten movies were captured, showing each
3-D spheroid prior to and during opto-genetic excitation. These movies were
used previously [3], finding when cells are pulsed with optogenetics (from 4
to 8 hours), cells do less collective rotation. Here we explore more closely
the relation between ERK signaling patterns and the resulting cellular
velocity. The ERK and velocity SSF kymographs are shown in supplementary
figure 2. Using the same method as in section II-A finds that the spatial ERK
signaling pattern is highly linearly correlated with the resulting cell
velocity pattern ($\rho=0.99,p<1e-12$).

Ten movies of MCF10A spheroids were imaged at 5 minutes per frame in 3-D
(512x512x135 voxels) for 12 hours (144 frames). ERK and H2B were imaged. Eight
movies had four hours of non-excitation, followed by four hours of excitation
every 30 minutes and then four hours of non-excitation. Movie 2021719_6h had
six hours of non-excitation, followed by six hours of excitation every 30
minutes and then two hours of non-excitation. Movie 20210721_2h had two hour
intervals for pre-, during- and post-excitation. Each movie had separate
kymographs for pre- and excitation. Movie 20210720_4h imaged a spheroid that
shifted partially out of frame between pre- and excitation. The resulting
distances in the embedded space were large due to the change in visual
appearance. Interestingly, the magnitude of the Euclidean distances in the
RKHS for velocity and ERK still follow the same linear relationship seen with
the other nine movies.

Figure 5 shows average velocity pre- and excitation. Figure 3 shows the 20
pre- and excitation movies in the NCD-RKHS. Each movie is represented by one
point for pre- and one for excitation. In Fig 5 A and B, note the lines
connecting these points. The length of these lines represents the visual
difference in signaling patterns between each pre- and excitation movie pair
as measured by the NCD. We find that the length of each line in the ERK RKHS
is highly linearly correlated with the length of the same line in the velocity
RKHS ($\rho=0.99,p=1e-12$). The RKHS here is key in enabling the Euclidean
distance in the embedded space to measure quantitative relationships on the
input image space.

(a)

(b)

(c)

Figure 5:  ERK signaling patterns are significantly correlated with cell
velocity. Each of ten movies $X_{1}..X_{10}$ are split into 20 pre-excitation
and excitation movies
$\\{X_{1}^{pre},X_{1}^{excite}...X_{10}^{pre},X_{10}^{excite}\\}$. We generate
the 3-D SSF kymographs for ERK and velocity for each movie and then compute
the pairwise NCD matrix and its resulting embedding in a 2-D RKHS $(k1,k2)$.
Because there are two classes (pre- and during-excitation), the CSF value for
each embedding is computed in the 2-D space using the ground truth labels to
show slightly more structure from the ERK (a) _vs._ velocity (b) embeddings.
Compute the Euclidean distance in the 2-D embedding space between each pre-
and excitation movie pair for ERK, $d_{ERK}$ and velocity, $d_{velocity}$,
indicated by the connecting lines. The linear correlation between the
difference magnitude for each pre- vs. excitation movie in ERK vs. velocity
(c) is $\rho=0.99$ ($p=1e-12$). The change in ERK signaling pattern induced by
optogenetic manipulation is strongly predictive of the resulting change in
cellular velocity.

###  II-D Synthetic spatio-temporal signaling patterns

Synthetic kymographs simulating 2-D+time constant velocity cellular motion for
three simulated classes were created to characterize spatio-temporal pattern
extraction. Each synthetic kymograph comprises ten cell trajectories moving
with random constant velocity per trajectory normally distributed randomly
generated from a per class mean $\in[1,3,5]$ and standard deviation
$\sigma=0.5$ . Figure 6 shows maximum intensity projections of sample phantom
kymographs for each velocity class (A, B, and C). At each cell centroid the
velocity value is recorded on one channel and a random value $\in[1,255]$ is
recorded on a second channel. 100 synthetic kymographs were generated per
class.

Pairs of the synthetic 3-D kymographs are input to the FLIF 3-D compression to
compute the pairwise NCD matrix, first using the velocity channel. This is
repeated with the random value channel to create two RKHS embeddings, one for
constant velocity values and one for random values. The cluster structure
function (CSF) is computed in the RKHS embeddings, preserving the principal
three dimensions of the RKHS embedding as there are three ground truth
classes. The CSF is computed against these ground truth classes. For the
constant velocity kymographs (D), the CSF is $[0.91,0.052]$ (the first value
is the mean CSF across the three clusters, the second value is the standard
deviation). For the random value kymographs (E), the CSF is $[1.3,0.24]$,
indicating that embedding preserves some class structure but captures less
meaningful information compared with the constant velocity kymographs.

Figure 6:  Synthetic kymographs characterize spatio-temporal pattern
extraction. Kymographs were generated from simulated 2-D+time cellular motion
for constant velocity tracks. The kymograph stores the velocity on channel one
and a random value on channel two at each cell centroid. Three different
velocity classes were simulated (A, B, and C). This embedding for the velocity
signaling channel (D) shows clear separation among the classes. The random
channel embedding (E) preserves some class structure but is a significantly
less meaningful embedding as measured by the cluster structure function (CSF).
The CSF is reported as the mean and standard deviation of the per cluster
optimality deficiency against the ground truth labels.

##  III Discussion

Current approaches to finding spatiotemporal patterns of cell signaling in
timelapse microscopy movies start with an empirical observation of the
pattern, using the computational tools to measure or classify the observed
pattern. The cytonuclear ratio used to measure KTR activation requires a
careful post-processing of the measurements. One class of approaches uses
neural networks and/or non-linear embedding, but the non-linearity of these
methods makes spatial measurements from the resulting embeddings unreliable.
Another class of approaches uses features, segmenting the movie to extract the
observed patterns and generate measurements. While both types of approaches
remain valuable, there is a need for unbiased techniques that identify
patterns with no _a priori_ knowledge in order to quantify patterns that are
not easily visible to the human eye.

The approach presented here is a general purpose tool for identifying and
measuring spatiotemporal signaling patterns in live cell microscopy. The SSF
is a new Kolmogorov structure function that measures signaling activity for
any configuration of nuclear voxel intensities w.r.t. the surrounding
cytoplasm. The SSF is an optimal mapping from the 2-D $(cytoplasmic,nuclear)$
intensity space to a 1-D activation measure on [-1.0,1.0]. The currently
widely used cytonuclear ratio [5, 3] is a non-linear and unbounded function
that is not an effective measure of cell signaling state. Posing the SSF as a
Kolmogorov structure function[13, 9] defines it as an optimally meaningful
measure of signaling activation. The SSF kymographs represent the
spatiotemporal signaling patterns for a given movie (or subset) in a lossless
manner for 2-D+time movies, and via a spatial projection for 3-D+time movies.
We leverage here the effective FLIF lossless 3-D image compression algorithm
as a metric distance (NCD)[17, 18], and use that to embed each movie in a
reproducing kernel Hilbert space (RKHS)[10]. This RKHS preserves optimally in
the low dimensional embedding space the characteristics of the image
similarity measured by the FLIF compression among the input images. This
includes the Euclidean distance throughout the space, a property used here to
estimate the cluster structure function (CSF) for data that would otherwise be
too large for effective compression. The embedding into the RKHS is completely
unsupervised, once the data is in the RKHS it is well suited for unsupervised
clustering or for training supervised learners on subsets of the input data.

For the 2-D human breast epithelial cell (MCF10A) movies, the RKHS embedding
captured the pattern differences among six different oncogenic mutations
associated with distinct signaling patterns, with the CSF measuring the
structure was significantly higher due to ERK compared to AKT. The NCD kernel
measured directly the relation between ERK signaling and cellular velocity
patterns, and found significant increases in the relationship for the
PIK3CA_H1047R and PIK3CA_E545K mutations. For colonies of induced human
pluripotent stem cells, it was found that there was no significant difference
in ERK dynamics between self-renewing and differentiating colonies, and no
change in the relation between ERK and velocity. Using optogenetic excitation
while filming live in the microscope, with 3-D MCF10A organoids we see that
the magnitude of the change in ERK signaling pattern during excitation is
almost perfectly predictive of the change in cellular velocity $\rho=0.99$
($p=1e-12$). Finally, a synthetic dataset using velocity and random signaling
shows the sensitivity in measuring the relation between motion patterns and
signaling dynamics.

Moving forward, an extension of FLIF to 4-D would allow processing of 3-D+time
kymographs with no need for dimensionality reduction. Another area for
possible improvement is in the quantization step. The current approach is to
quantify the SSF kymograph from $[-1.0,1.0]$ to an 8-bit unsigned integer
representation on $[1,255]$. In other work using the NCD with timelapse
microscopy movies has benefited from varying the number of quantization
symbols[19, 17]. An optimization search on the number of quantization symbols
for the movies like those processed here would be time consuming but possible
in future work. Another computationally demanding but useful feature would be
to generate random reference kymographs by writing random values at the cell
centroid locations and iterating to quantify the contribution of motion
patterns _vs._ signaling dynamics in the embedding structure as measured by
the CSF.

##  IV Software and data availability

All of the software tools used are available free and open source, see
https://git-bioimage.coe.drexel.edu/opensource/ssfCluster. The image data
together with segmentation and tracking results can be viewed interactively at
https://leverjs.net/ssfCluster. The LEVERSC 4-D WEBGL viewer [20] renders 3-D
kymographs and images, and the web API also supports downloading metadata and
results directly.

##  V Methods

###  V-A Live cell imaging

####  V-A1 Human breast epithelium (MCF10A) monolayers

Wild-type human mammary epithelial cells MCF10A cells were a gift of Joan S.
Brugge, Harvard Medical School, Boston, MA. AKT1-E17K, PIK3CA E545K, PIK3CA
H1047R knockin, and MCF10A-PTEN deletion (-//-) knockout derivatives of
parental MCF10A cell line[21] was a gift of Ben Ho Park, Johns Hopkins
University, USA. ErbB2 overexpressing MCF10A cell line was generated by
lentiviral transduction of pHAGE-ERBB2 construct (a gift from Gordon Mills &
Kenneth Scott, Addgene plasmid #116734[22]). Transduction was performed in the
presence of 8 $\mu$g//ml polybrene (TR1003, Sigma) in MCF10A WT cells already
expressing H2B-miRFP703 and ERK-KTR-mTurquoise2 biosensors. Cells were
selected with 5 $\mu$g//ml puromycin (P7255, Sigma).

MCF10A cells were cultured in growth medium composed by DMEM:F12 supplemented
with 5% horse serum, 20 ng/ml recombinant human EGF (Peprotech), 10 mg/ml
insulin (Sigma-Aldrich/Merck), 0.5 mg/ml hydrocortisone (Sigma-Aldrich/Merck),
200 U/ml penicillin and 200 mg/ml streptomycin. All the experiments were
carried out in starvation medium consisting of DMEM:F12 supplemented with 0.3%
BSA (Sigma-Aldrich/Merck), 0.5 mg/ml hydrocortisone (Sigma-Aldrich/Merck), 200
U/ml penicillin and 200 mg/ml streptomycin. Cells were growth factor and serum
starved by removing growth medium, washing the monolayers 2 times with PBS and
adding starvation media.

The stable nuclear marker H2B-miRFP703 was a gift from Vladislav Verkhusha
(Addgene plasmid #80001)[23], and subcloned in the PiggyBac plasmid
pPBbSr2-MCS. ERK-KTR-mTurquoiose2 and ERK-KTR-mRuby2 sequences were
synthesized (GENWIZ) by fusing the ERK Kinase Translocation Reporter (ERK-
KTR)[5] CDS with mTurquoiose2[24] and mRuby2[25] CDSs, respectively.
FoxO3a-mNeonGreen sequence was synthesized (GEN- WIZ) by fusing the 1-1188
portion of the homo sapiens forkhead box O3 a (FoxO3a) CDS with mNeonGreen
CDS, a green fluorescent protein derived by Branchiostoma lanceolatum[26].
ERK-KTR-mTurquoiose2 and FoxO3a-mNeonGreen were cloned in the PiggyBac
plasmids pMP-PB, pSB-HPB (gift of David Hacker, Lausanne[27]) or pPB3.0.Blast,
an improved PiggyBac plasmid generated in Olivier Pertz’s lab. For stable DNA
integration PiggyBac plasmids were transfected together with the helper
plasmid expressing the transposase[28]. To generate cell lines stably
expressing nuclear marker and biosensors, transfection was carried out with
FuGene (Promega).

Stable clones expressing biosensors were selected using Puromycin (P7255,
Sigma), Blasticidin S HCI (5502, Tocris), and Hygromycin B (sc-29067, Lab
Force) and imaging-based screening.

MCF10A cells and knock-in/out derivatives were plated on Fibronectin (PanReac
AppliChem) coated (0.25ug/cm2) 96 well 1.5 glass bottom plates (Cellvis) at 30
000 cells/ well density and allowed to adhere and form monolayer in growth
media. Cells were starved for 48h before starting the experiments. In drug
perturbation experiments, starved cells were imaged for 5h, then indicated
drugs or vehicle (DMSO) control was added and imaging was resumed for 15h.

Imaging was done on an epifluorescence Eclipse Ti inverted fluorescence
microscope (Nikon) controlled by NIS-Elements (Nikon) with a Plan Apo air 20X
(NA 0.8) objective. Laser-based autofocus was used throughout the experiments.
Image acquisition was performed with an Andor Zyla 4.2 plus camera at a 16-bit
depth. Illumination was done with a SPECTRA X light engine (Lumencor) with the
following excitation and emission filters (Chroma): far red (miRFP703): 640nm,
ET705/72m; red (mRuby2): 555nm, ET652/60m; green/yellow (mNeonGreen): 508nm,
ET605/52; cyan (mTurquoise2): 440nm, HQ480/40.

####  V-A2 Human induced pluripotent stem cells

Maintenance WA09 (H9) hESC line was purchased from WiCell (wicell.org) and
maintained in Essential 8 flex medium (A2858501, Thermo Fisher Scientific) on
hESC-qualified growth factor reduced Geltrex-coated (A1413302, Thermo Fisher
Scientific) 6 well plates. Cells were split into 6 well plates at 1:10 ratio
when cells become confluent using 0.5 mM EDTA. Medium was changed according to
the E8 flex protocol. Time-lapse imaging and differentiation of hESCs
Established multi-colour hPSC expressing ERKKTR-mClover, ORACLE OCT4tdtomato
and H2BmiRFP were plated onto Geltrex-coated 24-well plates (P24-1.5H-N,
CellVis) a day before imaging supplemented with Essential 8 flex medium and
maintained in the incubator. Colonies were imaged using a Nikon Ti2 with a
Yokogawa CSU-W1 spinning disk system. ERKKTR-mClover, ORACLE-OCT4tdtomato and
H2BmiRFP were captured every 5min using 488, 561 and 642nm laser respectively
using Nis Elements Nd acquisiation modality using 2x2 large image with a 10%
overlap and optimal path blending. To initiate differentiation to
neuroectodermal lineage, essential 8 flex medium was changed to PSC neural
induction medium (A1647801, Thermo Fisher Scientific) 4h prior to imaging and
changed every other day.

####  V-A3 Optogenetic manipulation of 3-D MCF10A spheroids

Optogenetic stimulation of mammary acini

Mammary acini were grown from wild-type human female mammary epithelial MCF10A
cells. The cells were stably modified by using the PiggyBac transposon system
to express H2B-miRFP703, ERK-KTR-mRuby2 and Lyn-cytoFGFR1-PHR-mCit (OptoFGFR),
as previously described[3]. For acini formation, MCF10A single-cell
suspensions were mixed with 4 volumes of growth factor-reduced Matrigel
(Corning) at 4∘ C and spread evenly on the surface of glass bottom cell
culture plates at a concentration of $1.4x10^{4}\frac{cells}{cm^{2}}$. The
acini were cultured in DMEM/F12 supplemented with 2% horse serum, 20 ng/ml
recombinant human EGF, 0.5 mg/ml hydrocortisone, 10 mg/ml insulin, 200 U/ml
penicillin and 200 mg/ml streptomycin. Horse serum, insulin and EGF were
removed after 3 days of culture. For live imaging, 25 mM Hepes was added to
the medium prior to mounting on the microscope. Images of acini were acquired
on an epifluorescence Eclipse Ti2 inverted fluorescence microscope (Nikon)
equipped with a CSU-W1 spinning disk confocal system (Yokogawa) and a Plan Apo
VC 60X water immersion objective (NA = 1.2). For time-lapse imaging, laser-
based autofocus was used. Images were acquired with a Prime 95B or a Prime BSI
sCMOS camera (both Teledyne Photometrics) at 16-bit depth. Temperature, CO2
and humidity were controlled throughout live imaging with a temperature
control system and gas mixer (both Life Imaging Services). The following
lasers were used for excitation: 638 nm for far red/miRFP and 561 nm for
red/mRuby2. For the optogenetic stimulation with OptoFGFR, acini were
illuminated with wide field blue light (470 nm LED) for 100 ms at 50% LED
intensity at defined time points during spinning disc time-lapse imaging. Both
microscopes were controlled by NIS elements (Nikon).

For the optogenetic stimulation with OptoFGFR, acini were illuminated with
wide field blue light (470 nm LED) for 100 ms at 50% LED intensity at defined
time points during spinning disc time-lapse imaging. Live z-stack images of
individual acini were acquired every 5 minutes with a 0.6 µm z resolution.

###  V-B Segmentation and tracking

The segmentation and tracking used here are based on our previously developed
LEVER (lineage editing and validation) tools for segmentation and tracking[29,
30, 31, 32, 33, 34, 35]. The segmentation includes a non-local means
denoising[36], followed by a thresholding and a separation of touching cells.
The thresholding uses the SSF values from the H2B channel with an empirical
threshold of 0.01 that is used across all of the movies analyzed to date. The
thresholding and cell detection uses the H2B channel to identify cell
centroids. For each additional image channel (e.g. ERK-KTR, AKT-KTR), we add
the maximum magnitude of the positive and negative LoG responses from the
thresholded regions identified on the H2B channel, (as in section V-C) and
then use a watershed transform on the combined SSF channel images to separate
among touching cells. This allows the KTR signal, if available, to assist in
the most challenging of the segmentation tasks, separating adherent cells. The
multi-temporal association tracking (MAT) [30, 35] can be used to filter cells
for inclusion in the SSF kymograph, but is used here primarily to
automatically identify and correct segmentation errors [19, 29, 32]. The key
outcome from the segmentation and tracking are the $(x,y,z,time)$ centroid
locations for each cell. The approach is unsupervised, requiring no training
data and taking as its only parameter a range of radii to be used with the LoG
GPU filter to compute the SSF (section V-C), set at $[4:0.5:6]\mu m$ for all
of the human stem and cancer cells movies analyzed here.

###  V-C Computing the cell signaling structure function

The cell signaling structure function (SSF) measures the cell signaling state
as the intensity of the nuclear voxels w.r.t. the surrounding cytoplasm,
useful for any cell imaging protocol, including the powerful new collection of
biosensors known as kinase translocation reporters (KTRs) [5]. Current state-
of-the-art approaches to computationally analyzing these KTR signals rely on
the ratio of cytoplasmic to nuclear intensity, the cyto-nuclear ratio. The
cyto-nuclear ratio is non-linear, making it a poor choice for reporting cell
signaling state. This limitation has been noted in previous work, requiring
careful selection of acceptance regions for the KTR signals. An alternative
approach to computing KTR activation follows from associating image channel
intensities with reporter concentrations. Figure 3 compares the SSF to the
cyto-nuclear ratio. We consider images with intensity values normalized to
[0,1]. Given the average cytoplasmic pixel intensity $c$ and the average
nuclear pixel intensity $n$, define the cell signaling structure function
$H_{SSF}$ as in eqn. I.2. The function $H_{SSF}$ is on [-1.0,1.0], and varies
linearly with the nuclear and cytoplasmic intensities. In addition to having a
linear response to variations in cytoplasmic and nuclear reporter
concentrations, $H_{SSF}$ is equal to zero for cells that are unlabeled ($n$ =
0 and $c$ = 0). The $H_{SSF}$ can be robustly approximated via convolution
with a Laplacian of Gaussian (LoG) filter. The LoG filter simultaneously
estimates $n$ and $c$ while computing Eqn. I.2. The value of $H_{SSF}$ is
found at each cell centroid and unlike the nuclear/cytoplasmic ratio does not
require an accurate segmentation boundary. Finding the cell centroid is an
easier task compared to finding the nuclear boundary [34, 33].

The Laplacian of Gaussian (LoG) is a widely used blob enhancing filter [36].
The LoG combines the smoothing of a Gaussian filter with the edge enhancement
of the Laplacian operator. We define the SSF using the LoG response on each
imaging channel $\lambda\in\Lambda$ evaluated at each cell centroid and
radius. Here we use the LoG response combined across all imaging channels for
the segmentation. Any negative LoG response on the H2B channel is considered
foreground. Details of the segmentation are given in the methods section. The
cell signaling structure function $H_{SSF}$ can be written as

$H_{SSF}(\lambda\in\Lambda)=\frac{LoG(x_{c},y_{c},z_{c},\lambda,t,r_{LoG})}{{LoG}_{ref}(r_{SEG})}.$
(V.1)

$LoG_{ref}$ is the maximum value obtained filtering a zero-intensity spheroid
of radius $r_{SEG}$ against a full intensity background. The locations
$x_{c},y_{c},z_{c},t$ are cell centroids identified by the segmentation
algorithm. In practice, we omit the normalization by $LoG_{ref}$ as this is a
constant term due to the multi-resolution LoG implementation as described
below, and allow the normalization to [-1.0,1.0] to happen during the
quantization step as in Section V-D.

The Laplacian of Gaussian filter is a blob enhancing filter that combines a
Gaussian smoothing with a Laplacian edge detection. We recently developed a
GPU-based implementation of the Laplacian of Gaussian filter that works in 3-D
using NVIDIA’s CUDA parallel programming toolkit [14]. This filter is scale
invariant, meaning that its output remains similar across different radii.
This implementation uses axis-aligned spherical approximations to compute the
blob response efficiently at every voxel in the image. It is possible to
compute the LoG response for arbitrary elliptically oriented blobs, but the
extra computational requirements have not been needed for the applications
considered to date. The LoG filter takes a single parameter of radius. The
radius of the filter relates to the standard deviation of the underlying
kernel as $r=\sigma\ast\sqrt{d}$ where $d$ is 2 for 2-D images and 3 for 3-D
images. Define a Gaussian kernel $G(x,y,z)$,

$G(x,y,z)=\frac{1}{\sqrt{(2\
\pi)^{d}\sigma_{x}^{2}\sigma_{y}^{2}\sigma_{z}^{2}}}\exp^{-\frac{1}{2}((\frac{x}{\sigma_{x}})^{2}+(\frac{y}{\sigma_{y}})^{2}+(\frac{z}{\sigma_{z}})^{2})}.$
(V.2)

Then write the scale invariant Laplacian of Gaussian as

$LoG(x,y,z)={((x,y,z)}^{T}\Lambda^{-2}(x,y,z)-d)G(x,y,z)$ (V.3)

where

$\Lambda=\left[\begin{matrix}\sigma_{x}&0&0\\\ 0&\sigma_{y}&0\\\
0&0&\sigma_{z}\\\ \end{matrix}\ \ \right].$ (V.4)

V.3 is efficient to compute because we omit the covariance terms so
$LoG(x,y,z)$ can be computed as a combination of 1-D LoG and Gaussian kernels
across the d-dimensional image as

$LoG(x,y,z)=LoG(x)G(y)G(z)+LoG(y)G(x)G(z)+LoG(z)G(y)G(x)$ (V.5)

This implementation is faster compared to computing a full 3-D kernel as would
be required for non-diagonal covariance matrices representing non-axis aligned
ellipses. The response $LoG(x,y,z)$ is normalized so that the kernel always
sums to zero (even for kernels that protrude from the image) reducing
filtering artifacts at image boundaries.

###  V-D Quantizing SSF kymographs

We define the _SSF kymograph_ as a 3-D image whose axes are $(x,y,time)$. The
SSF kymograph uses the cell centroid locations obtained from the segmentation
algorithm, and optionally uses the tracking to select a subset of trajectories
(see Section V-B). At each cell centroid location, the SSF value for that cell
at the given spatiotemporal location and channel is written. For 2-D+time
movies, the 3-D SSF kymograph is an exact representation. For 3-D+time movies,
as in the optogenetic excitation of 3-D human breast epithelial cell (MCF10A)
spheroids, a maximum intensity projection along the $z$ axis reduces the
spatial dimensions to 2-D so our kymographs can be compressed using the 3-D
FLIF compressor. For kymographs where the minimum cell radius in pixels is
larger than 2 to avoid occlusions we compute the kymographs by down-sampling
in each spatial dimension by a factor of 2 to reduce compression time
requirements. Given kymographs $M_{1}$ and $M_{2}$, we compute the normalized
compression distance as
$NCD(M_{1},M_{2})=\frac{Z(M_{1}||M_{2})-\min(Z(M_{1}),Z(M_{2}))}{\max(Z(M_{1}),Z(M_{2}))}$
where $Z(M_{1})$ is the size in bytes of the FLIF compression file for $M_{1}$
and $||$ is the concatenation operator[7]. The FLIF compression requires 8-bit
input images, so the kymographs must be quantized. Quantization is also
required for color mapping the SSF values so the kymographs can be visualized,
as in Figure 2B,C and the supplementary figures. For a given experiment
containing $N$ 5-D movies, we generate $N$ 3-D kymographs and extract all non-
zero SSF voxels from the $N$ kymographs into a multiset $V$. The quantization
bins are then defined as 254 linearly spaced boundaries on
$(\mu(V)-\sigma(V),\mu(V)+\sigma(V))$, where $\mu(V)$ is the average SSF value
and $\sigma(V)$ is the standard deviation. The quantization maps the SSF
values to $[1,255]$, reserving the value of zero to indicate an absence of
signaling. When assigning a colormap to the the quantized value, the ‘parula’
colormap [36] that is widely available has a bimodal appearance helpful in
visualizing the 2-D SSF as a combined color signal and has been used for
visualizing kymographs throughout. The FLIF compression algorithm is color
agnostic, relying only on entropy calculations on the 8-bit quantized
kymographs for pattern detection. The FLIF compression also supports RGB input
images, allowing us to compute the NCD between kymographs containing up to
three imaging channels simultaneously, but the results to date only utilized
single channel kymographs as input to the NCD.

###  V-E The cluster structure function in the RKHS

The cluster structure function measures how meaningfully a given partition
into clusters represents an input dataset. The theory is Kolmogorov
complexity, an absolute measure of information content within and between
digital objects. For brevity, we show here only the compression approach used
to compute the CSF and omit the theory background, for details consult the
recent papers and the textbook[37, 9, 18]. Our input is a collection of $N$
5-D $(x,y,z,channel,time)$ microscopy movies $X=(x_{1},...,x_{N})$. Given a
partitioning of $X$ into $K$ clusters as $X={Y_{1},...,Y_{K}}$, with the
$i^{th}$ cluster $Y_{i}=(y_{1},...y_{m})$. For each $y_{i}$, write the
optimality deficiency $\delta(Y,y_{i})=Z(Y)+log(|Y|)-Z(y_{i})$, where $Z(Y)$
is the size in bytes of the compressed $Y$ and $|Y|$ is the cardinality of Y.
The average of the optimality deficiency for each of the $m$ elements in the
cluster is computed, and that average is taken again across all $K$ clusters
to compute a mean and standard deviation of the per cluster optimality
deficiencies. When used for 2-D images or smaller 3-D images, the approach is
to minimize this CSF in order to select the optimal value of $K$, addressing
in an absolute sense the optimal number of clusters in the given dataset.

We compute and use the CSF in this work using a different approach than
originally proposed[9]. For the 2-D MCF10A human breast epithelium monolayer
movies, there are 24 or 25 movies from each genetic condition. To compute
$Z(Y)$ for this data would require compressing the concatenation of all those
movies, intractable due to the image sizes involved. Instead of compressing
the movies in each cluster together, we use the pairwise NCD results to embed
the movies in an RKHS and then compute the CSF in the RKHS as follows. For
each cluster in the RKHS $Y^{\prime}=(y^{\prime}_{1},...y^{\prime}_{m})$, we
compute the centroid of the points in that cluster as a point in
$\mu_{Y^{\prime}}=\frac{\sum{y^{\prime}_{i}}}{m}$. For each optimality
deficiency, we use the Euclidean distance between each of the points
representing individual movies for the optimality deficiency,
$\delta(Y^{\prime},y_{i}^{\prime})=|y_{i}-\mu_{Y^{\prime}}|$. We compute the
CSF for every input kymograph, and then test statistically significant
differences between CSF values, as in ERK vs. AKT for the 2-D MCF10A human
breast epithelial cells, using the non-parametric Wilcoxon signed rank test
for significant differences of paired medians [38] for the $p$-values reported
throughout.

Our application of the CSF in the present work is also used to evaluate how
meaningful a given embedding is, rather than the originally proposed use for
selecting the optimal number of clusters. Instead of clustering in the RKHS
space, we use the ground truth labels for each movie as the partitioning of
the data, and choose as $K$ the true number of classes from the experimental
conditions. Then for the given ground truth partitioning and $K$ we compute
the CSF as a criterion function to evaluate the effectiveness of a particular
input, _e.g._ as in evaluating how meaningful ERK vs. AKT in inducing a given
structuring of the data for the human breast epithelial cells in Section II-A.

##  VI Acknowledgements

The authors are grateful to Ben Ho Park for providing AKT1E17K, PIK3CAE545K,
PIK3CAH1047R knockin, and PTEN(-/-) knockout MCF10A - derived cell lines[21].
The authors are also grateful to Rafael Carazo Salas and Sungmin Kim from the
Univ. Bristol, UK, for the imaging constructs used in the hIPSC dataset. This
project has received funding from the European Union’s Horizon 2020 research
and innovation program under the Marie Skłodowska-Curie grant agreement No
896310 to Agne Frismantiene. We acknowledge support of the Microscopy Imaging
Center of the University of Bern (https://www.mic.unibe.ch/). Portions of this
work were supported by Human Frontiers Science Grant RGP0043/2019
(Pertz/Cohen).

## References

  * [1] Paolo Armando Gagliardi et al.  “Collective ERK/Akt activity waves orchestrate epithelial homeostasis by driving apoptosis-induced survival”  In _Developmental Cell_ 56.12 Elsevier BV, 2021, pp. 1712–1726.e6  DOI: 10.1016/j.devcel.2021.05.007
  * [2] Paolo Armando Gagliardi et al.  “Automatic detection of spatio-temporal signaling patterns in cell collectives”  In _Journal of Cell Biology_ 222.10 Rockefeller University Press, 2023  DOI: 10.1083/jcb.202207048
  * [3] Pascal Ender et al.  “Spatiotemporal control of ERK pulse frequency coordinates fate decisions during mammary acinar morphogenesis”  In _Developmental Cell_ 57.18 Elsevier BV, 2022, pp. 2153–2167.e6  DOI: 10.1016/j.devcel.2022.08.008
  * [4] Jeremy E. Purvis and Galit Lahav  “Encoding and Decoding Cellular Information through Signaling Dynamics”  In _Cell_ 152.5 Elsevier BV, 2013, pp. 945–956  DOI: 10.1016/j.cell.2013.02.005
  * [5] Sergi Regot et al.  “High-Sensitivity Measurements of Multiple Kinase Activities in Live Single Cells”  In _Cell_ 157.7 Elsevier BV, 2014, pp. 1724–1734  DOI: 10.1016/j.cell.2014.04.039
  * [6] Takamasa Kudo et al.  “Live-cell measurements of kinase activity in single cells using translocation reporters”  In _Nature Protocols_ 13.1 Springer ScienceBusiness Media LLC, 2017, pp. 155–169  DOI: 10.1038/nprot.2017.128
  * [7] Rudi Cilibrasi and Paul M. Vitányi  “Clustering by Compression”  In _IEEE Transactions on Information Theory_ 51.4 Institute of ElectricalElectronics Engineers (IEEE), 2005, pp. 1523–1545  DOI: 10.1109/tit.2005.844059
  * [8] Jon Sneyers Wuille and Pieter  “FLIF: Free lossless image format based on MANIAC compression”  In _2016 IEEE International Conference on Image Processing (ICIP)_
  * [9] Andrew R. Cohen and Paul M. Vitányi  “The cluster structure function”  In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ Institute of ElectricalElectronics Engineers (IEEE), 2023, pp. 1–13  DOI: 10.1109/tpami.2023.3264690
  * [10] Jonathan H. Manton and Pierre-Olivier Amblard  “A Primer on Reproducing Kernel Hilbert Spaces”  In _Foundations and Trends in Signal Processing_ 8.1-2 Now Publishers, 2015, pp. 1–126  DOI: 10.1561/2000000050
  * [11] S.. Kamvar, D. Klein and C.. Manning  “Spectral learning”  In _International Joint Conference of Artificial Intelligence_ , 2003, pp. 561–566
  * [12] Andrew Y. Ng, Michael Jordan and Yair Weiss  “On Spectral Clustering: Analysis and an algorithm”  In _Advances in Neural Information Processing Systems 14_ , 2002
  * [13] Nikolai Vereshchagin and Paul M. Vitányi  “Kolmogorov’s Structure Functions and Model Selection”  In _IEEE Transactions on Information Theory_ 50.12 Institute of ElectricalElectronics Engineers (IEEE), 2004, pp. 3265–3290  DOI: 10.1109/tit.2004.838346
  * [14] Wait et al.  “Hydra image processor: 5-D GPU image analysis library with MATLAB and python wrappers”  In _Bioinformatics_ 35.24, 2019, pp. 5393–5395  DOI: 10.1093/bioinformatics/btz523
  * [15] Paul Vitányi  “Meaningful Information”  In _IEEE Transactions on Information Theory_ 52.10 Institute of ElectricalElectronics Engineers (IEEE), 2006, pp. 4617–4626  DOI: 10.1109/tit.2006.881729
  * [16] Marc-Antoine Jacques et al.  “CODEX, a neural network approach to explore signaling dynamics landscapes”  In _Molecular Systems Biology_ 17.4 EMBO, 2021  DOI: 10.15252/msb.202010026
  * [17] Andrew R. Cohen et al.  “Automatic Summarization of Changes in Biological Image Sequences Using Algorithmic Information Theory”  In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 31.8 Institute of ElectricalElectronics Engineers (IEEE), 2009, pp. 1386–1403  DOI: 10.1109/tpami.2008.162
  * [18] A.. Cohen and P… Vitanyi  “Normalized Compression Distance of Multisets with Applications”  In _Pattern Analysis and Machine Intelligence, IEEE Transactions on_ 37.8, 2015, pp. 1602–1614  DOI: 10.1109/TPAMI.2014.2375175
  * [19] A.. Cohen et al.  “Computational prediction of neural progenitor cell fates”  In _Nat Methods_ 7.3, 2010, pp. 213–218
  * [20] Winter, Mark, Cohen and Andrew R.  “LEVERSC: Cross-Platform Scriptable Multichannel 3-D Visualization for Fluorescence Microscopy Images”  In _Frontiers in Bioinformatics_ 2, 2022  DOI: 10.3389/fbinf.2022.740078
  * [21] John P. Gustin et al.  “Knockin of mutant PIK3CA activates multiple oncogenic pathways”  In _Proceedings of the National Academy of Sciences_ 106.8 Proceedings of the National Academy of Sciences, 2009, pp. 2835–2840  DOI: 10.1073/pnas.0813351106
  * [22] Patrick Kwok-Shing Ng et al.  “Systematic Functional Annotation of Somatic Mutations in Cancer”  In _Cancer Cell_ 33.3, 2018, pp. 450–462.e10  DOI: https://doi.org/10.1016/j.ccell.2018.01.021
  * [23] Daria M. Shcherbakova et al.  “Bright monomeric near-infrared fluorescent proteins as tags and biosensors for multiscale imaging”  In _Nature Communications_ 7.1 Springer ScienceBusiness Media LLC, 2016  DOI: 10.1038/ncomms12405
  * [24] Joachim Goedhart et al.  “Structure-guided evolution of cyan fluorescent proteins towards a quantum yield of 93%”  In _Nature Communications_ 3.1 Springer ScienceBusiness Media LLC, 2012  DOI: 10.1038/ncomms1738
  * [25] Amy J. Lam et al.  “Improving FRET dynamic range with bright green and red fluorescent proteins”  In _Nature Methods_ 9.10 Springer ScienceBusiness Media LLC, 2012, pp. 1005–1012  DOI: 10.1038/nmeth.2171
  * [26] Nathan C. Shaner et al.  “A bright monomeric green fluorescent protein derived from Branchiostoma lanceolatum”  In _Nature Methods_ 10.5 Springer ScienceBusiness Media LLC, 2013, pp. 407–409  DOI: 10.1038/nmeth.2413
  * [27] Sowmya Balasubramanian, Florian M. Wurm and David L. Hacker  “Multigene expression in stable CHO cell pools generated with the piggyBac transposon system”  In _Biotechnology Progress_ 32.5 Wiley, 2016, pp. 1308–1317  DOI: 10.1002/btpr.2319
  * [28] Kosuke Yusa et al.  “A hyperactive piggyBac transposase for mammalian applications”  In _Proceedings of the National Academy of Sciences_ 108.4 Proceedings of the National Academy of Sciences, 2011, pp. 1531–1536  DOI: 10.1073/pnas.1008322108
  * [29] Mark Winter et al.  “Vertebrate Neural Stem Cell Segmentation, Tracking and Lineaging with Validation and Editing”  In _Nature Protocols_ 6.12, 2011, pp. 1942–52  DOI: 10.1038/nprot.2011.422
  * [30] Winter et al.  “Axonal transport analysis using Multitemporal Association Tracking”  In _Int J Comput Biol Drug Des_ 5.1, 2012, pp. 35–48  DOI: 10.1504/IJCBDD.2012.045950
  * [31] Apostolopoulou et al.  “Non-monotonic Changes in Progenitor Cell Behavior and Gene Expression during Aging of the Adult V-SVZ Neural Stem Cell Niche”  In _Stem Cell Reports_ 9.6, 2017, pp. 1931–1947  DOI: https://doi.org/10.1016/j.stemcr.2017.10.005
  * [32] Winter et al.  “Computational Image Analysis Reveals Intrinsic Multigenerational Differences between Anterior and Posterior Cerebral Cortex Neural Progenitor Cells”  In _Stem Cell Reports_ 5.4, 2015, pp. 609–620  URL: http://www.sciencedirect.com/science/article/pii/S2213671115002362
  * [33] M. Winter et al.  “Separating Touching Cells using Pixel Replicated Elliptical Shape Models”  In _IEEE Trans Medical Imaging_ 38.4, 2018, pp. 883–893
  * [34] Maska et al.  “The Cell Tracking Challenge: 10 years of objective benchmarking”  In _Nat Methods_ , 2023  DOI: 10.1038/s41592-023-01879-y
  * [35] Chenouard et al.  “Objective comparison of particle tracking methods”  In _Nat. Methods_ 11, 2014  DOI: 10.1038/nmeth.2808
  * [36] R Gonzales, R Woods and S Eddins  “Digital Image Processing Using MATLAB”  Knoxville, TN: Gatesmark Publishing, 2009
  * [37] Ming Li and Paul Vitányi  “An Introduction to Kolmogorov Complexity and Its Applications”  Springer New York, 2008  DOI: 10.1007/978-0-387-49820-1
  * [38] Sergios Theodoridis and Konstantinos Koutroumbas  “Pattern Recognition, Fourth Edition”  Elsevier, 2009  DOI: 10.1016/b978-1-59749-272-0.x0001-2

(a)

(b)

Supplementary Figure 1:  2-D projections of 3-D ERK SSF kymographs for 10
differentiated (a) and 10 self-renewing (b) movies. The $X$ axis represents
the spatial dimension, and is obtained by taking a maximum intensity
projection along the $Y$ axis. The time axis in each panel spans 865 frames at
5 minutes per frame, 72 hours. The full dataset can be seen here:
https://leverjs.net/ssfCluster/HSC.

(a)

(b)

Supplementary Figure 2:  2-D projections of 3-D ERK (a) and velocity (b) SSF
kymographs for 3-D mammary acini, spheroids of human female mammary epithelial
MCF10A cells . The original 3-D+time movies are rendered to a 3-D SSF
kymograph using a maximum intensity projection along the $Z$ spatial axis. For
the 2-D rendering shown here, the spatial component is again projected via
maximum intensity along the $Y$ axis. Each movie is labeled by the time
(hours) before optogenetic excitation, the time (hours) that optogenetic
excitation lasts (pulses every 30 minutes), and the age of the organoid
(days). The dashed vertical lines indicate the beginning and end of the
optogenetic excitation. The 2-D kymographs shown here are useful for human
visualization, the full 3-D kymographs are input the FLIF compression
algorithm to compute the pairwise NCD to generate the reproducing kernel
Hilbert space embedding (5). The optogenetic excitation dataset can be seen
here: https://leverjs.net/ssfCluster/optoGenetic.  Supplementary Figure 3:
2-D projections of 3-D SSF kymographs for movies showing 2-D+time movies of
live monolayers of human breast epithelial cells (MCF10A). These 24 kymographs
are from one of six imaging experiments, showing five oncogenic mutations plus
wild type (one each per row), also shown as an RKHS embedding of Fig. 2D. The
first PIK3CA_1047R column is the same kymograph shown in Fig. 2B,C. The full
dataset contains 147 movies, and can be viewed here:
https://leverjs.net/ssfCluster/MCF10A_2D.  Supplementary Movie 1:  Animated
version of Figure 2. A timelapse movie showing ERK-KTR signaling in a
monolayer of human breast epithelial cells (MCF10A) from the PIK3CA_H1047R
mutation with cellular activation indicated by dark nuclei against bright
cytoplasm clearly propagating across the image (left panel). The 3-D kymograph
(center panel) is the input to the FLIF 3-D compression used with the
normalized compression distance to define the RKHS embedding, shown here with
the current image frame overlaid in gray. The 2-D kymograph (right) panel is a
projection of the 3-D kymograph to facilitate human visualization, with the
current timepoint indicated by the red line and the signaling patterns clearly
visible as diagonal yellow stripes of activation across the monolayer (Link).
Supplementary Movie 2:  Rotating view of 3-D kymograph from Figure 2B. The 3-D
SSF kymograph stores the SSF value at each cell centroid location
$(x,y,time)$. The normalized compression distance finds patterns of similarity
between pairs of SSF kymographs. Shown here is the SSF kymograph for the
PIK3CA_H1047R oncogenic mutation that contains distinct signaling patterns,
(Link).