Mechanically-driven Stem Cell Separation in Tissues caused by Proliferating
Daughter Cells

Johannes C. Krämer1, Edouard Hannezo2, Gerhard Gompper1, and Jens
Elgeti1$\star$

1 Theoretical Physics of Living Matter, Institute of Biological Information
Processing and Institute for Advanced Simulations, Forschungszentrum Jülich,
52425 Jülich, Germany,

2Institute of Science and Technology Austria, 3400 Klosterneuburg, Austria.

⋆<EMAIL_ADDRESS>
October 6, 2023

## Abstract

The homeostasis of epithelial tissue relies on a balance between the self-
renewal of stem cell populations, cellular differentiation, and loss. Although
this balance needs to be tightly regulated to avoid pathologies, such as tumor
growth, the regulatory mechanisms, both cell-intrinsic and collective, which
ensure tissue steady-state are still poorly understood. Here, we develop a
computational model that incorporates basic assumptions of stem cell renewal
into distinct populations and mechanical interactions between cells. We find
that the model generates unexpected dynamic features: stem cells repel each
other in the bulk tissue and are thus found rather isolated, as in a number of
in vivo contexts. By mapping the system onto a gas of passive Brownian
particles with effective repulsive interactions, that arise from the generated
flows of differentiated cells, we show that we can quantitatively describe
such stem cell distribution in tissues. The interaction potential between a
pair of stem cells decays exponentially with a characteristic length that
spans several cell sizes, corresponding to the volume of cells generated per
stem cell division. Our findings may help understanding the dynamics of normal
and cancerous epithelial tissues.

##  1 Introduction

Tissue renewal through cell division to balance the constant loss of cells is
a hallmark of life in multicellular organisms. It is widely accepted that for
most tissue types, stem cells (SC) play a key role in this complex process [1,
2]. Stem cells have the potential to proliferate and self renew over long
timescales, continuously generating various types of differentiated cells
required for physiological tissue function. Such tissue maintenance via small
populations of stem cells requires fine-tuned fate choices to ensure not only
a constant and well-defined ratio of cellular types, but also stable tissue
size and overall cellular numbers.

Equally important, however, is a homogeneous positioning of stem cells across
the tissue, such that lost cells can quickly be replenished. If all stem cells
were to segregate from progeny, tissue function and maintenance could be
compromised. Indeed, in a lattice model without mechanical interactions,
proliferating and non-proliferating cells have been shown to require
reversible differentiation in a manner dependent on relative local cell
numbers to avoid non-homeostatic behavior, such as phase separation of cell
types [3]. Other types of fate control that have been recently explored
theoretically and experimentally are an effective fate determinant field of a
diffusible chemical. It is homogeneously supplied but consumed by stem cells,
thus, providing a biochemical negative-feedback mechanism for stem cell
density, which can also reproduce homogeneous distribution of stem cells
within the tissue [4].

In recent years, mechanical forces have been shown to be crucial regulators
for tissue growth and homeostasis. A key concept is the homeostatic pressure,
which is the pressure at which cellular division and apoptosis are balanced,
under the generic assumption that these processes depend on mechanical forces.
Based on this concept, competition between different tissues [5], tissue
fluidization [6], negative bulk homeostatic pressure [7], interface dynamics
during competition [8, 9, 10], and the evolution of tissue [11] could be
explained. However, these studies always considered situations without the
possibility of conversion between cellular types and change in self-renewal
potential, which is by definition key to understand stem cell dynamics in
tissues. This raises the question of whether purely mechanical interactions,
which must occur in any confluent tissue with cell renewal and loss, could be
a factor in the regulation of stem cell dynamics and positioning.

Here, we study the interaction and distribution of stem cells in self-renewing
tissues in the presence of mechanical interactions only. As a minimal
hierarchical scheme of tissues with stem and differentiated cells (Fig. 1), we
take a classical model [12], where stem cells (SC) asymmetrically divide into
a stem and a transient amplifying (TA) cell. The TA cells can undergo a number
$N_{TA}$ of symmetric divisions before terminally differentiating into non-
proliferative cells (TD). Although this is an oversimplified picture [13], we
show that it already gives rise to complex dynamics and provides a
computational framework which can be readily generalized to more realistic
stochastic models of cell differentiation. For instance effects of the stem
cell environment, the so called niche, on fate decisions [14, 15] or self-
controlled fate decisions such as in "open niche" [16], where all cells have
the same differentiation potential.

A \begin{overpic}[width=433.62pt]{SC_Division_Scheme.png}
\put(-5.0,63.0){\captiontext*{}} \end{overpic}

B \begin{overpic}[scale={0.4}]{Tissue_Snapshot.png}
\put(-5.0,90.0){\captiontext*{}} \end{overpic}

Figure 1:  Stem cell (SC) division and differentiation scheme leads to well-
separated SC arrangement. (a) Schematic illustration of growth, division, and
differentiation implemented in the two-particle growth model. SCs grow by an
active growth force $F^{(G)}$ and divide when the two particles of a cell are
separated by a threshold distance $R_{ct}$. A SC divides always asymmetrically
in one SC daughter and one TA daughter cell. TA cells still grow, but their
maximal number of divisions are limited. They always divide symmetrically, and
in the final step both daughters are TD cells. TA and TD cells are removed
with the same, finite apoptosis rate $k_{a}$. (b) In bulk simulations, SCs are
distributed all over the tissue, and rarely at close distance. Each SC is
surrounded by a halo of TA and TD cells. The main mass is made up by non-
growing TD cells. The same color scheme as in Fig. 1 is applied with finer
gradation for TA cells, which can divide up to five times. The simulation was
performed in a squared two-dimensional box with boxlength $100~{}R_{PP}$.
After equilibration, the average cell density is $\rho\sim 1.58/R_{PP}$, and
the snapshot was taken after 25 generations.

Implementing the renewal scheme into the two particle growth model [17] (sec.
2 and App. B), we find that the resulting tissue is characterized by well
separated SCs surrounded by a small cloud of TA cells in a sea of TD cells in
two (see Fig. 1), and three dimensions (see sec. 6). The cause of the
homogeneous distribution of stem cells in the tissue on larger distances lies
in the of outflux of cells from the stem cells, and short-range mechanical
interactions. Surprisingly, we find that this repulsion and distribution of SC
in the tissue can be effectively described by a simple thermodynamic picture
with SCs behaving like soft repulsive colloids in a thermal bath. The length
scale of the soft repulsive potential is set by the volume of cells generated
per stem cell division.

In section 3, we first discuss the dynamics of single stem cells to determine
the length scale of the interaction, and to study how the cloud of their
progeny around them impacts their movements. Subsequently, in section 4, we
study the dynamics of SC pairs. From their distance distribution function, we
infer the amplitude of the effective repulsive potential that is mediated by
the outward flow of differentiated cells that they each produce. In section 5,
we show that the distribution of stem cells in a bulk tissue can be described
by this effective pair interaction. In three dimensions, we find that the
repulsive interaction of stem cells is also preserved (see sec. 6). Last, we
demonstrate independence of the homeostatic pressure controlled division
mechanism, by replacing it with a stochastic one, in section7. This shows that
our results are primarily driven by mechanical interactions and proliferation
generated outflow of cells.

##  2 Simulation model

Each cell in the simulation model consists of two particles, which interact
via a repulsive force to model growth. When the particle distance reaches a
size threshold, a cell divides and two new particles are placed at random next
to the old ones. Particles of different cells interact via an attractive and a
repulsive force, to model adhesion and volume exclusion, respectively. The
length $R_{PP}$, up to which inter-cellular interactions are present, defines
the length scale of our model and depicts the cell size. TA and TD cell
removal is modelled in a stochastic process with fixed rate $k_{a}$, which
defines the timescale of the model. Note that the division rate $k_{d}$ is not
fixed and will follow from the other model parameters. Random and dissipative
forces are included by dissipative particle thermostat. Cells follow a a
classical model [12], see Fig. 1, and have all identical interaction
parameters. Solely stem cells do not die ($k_{a}=0$), and TD cells do not
actively grow ($G=0$) to prevent divisions. More details on the model can be
found in App. B and previous works [17, 7, 8, 11, 10].

##  3 Isolated stem cell

To understand the separation of stem cells in bulk tissue, we begin our
analysis by studying the dynamics of a single SC and its offspring. We start
with an isolated stem cell and shortly thereafter find it surrounded by a
cloud of progeny that maintains a certain size throughout the simulation time.
(see Fig. 2).

A\begin{overpic}[scale={0.4}]{CellColonyInteractionRadiusSnapshot.png}\put(-3.0,97.0){\captiontext*{}}\end{overpic}

B\begin{overpic}[scale={0.4}]{2SC-interaction-
snapshot.png}\put(-2.0,97.0){\captiontext*{}}\end{overpic}

C \begin{overpic}[scale={0.4}]{CellColonyDistribution_char-length-
twinaxis.eps}\put(1.0,61.0){\captiontext*{}}\end{overpic}

D\begin{overpic}[scale={.4}]{CellColonyDistribution_density_isolated.eps}\put(0.0,74.0){\captiontext*{}}\end{overpic}

E\begin{overpic}[scale={0.4}]{TwoSCpaircorrelation_rdf.eps}\put(-1.0,65.0){\captiontext*{}}\end{overpic}

F\begin{overpic}[scale={0.4}]{TwoSCpaircorrelation_betaV0.eps}\put(-2.0,75.0){\captiontext*{}}\end{overpic}

Figure 2:  Interaction parameter determination. (a) Snapshot of isolated SC
and (b) two interacting SCs. (c) Cell number generated in each SC division
cycle, i. e. sum of TA cell and latest TD cells, for the upper limit
$N_{p}=2^{N_{TA}+1}$ (open circles, dotted line) and measured from simulation
by averaging over time $\langle N_{p}\rangle$ (filled circles, solid line) as
function of $N_{TA}$. The latter is used to determine the characteristic
progeny distance $L_{p}$, which is used on the right y-axis label. (d) Cell
density of TA (solid line) and all daughter cells (dashed line) around SC,
where $L_{p}$ derived from average measured progeny number $N_{p}$ is marked
by dotted line, for $N_{TA}=6$. (e) pair correlation $g(\Delta r)$ as function
of SC distance $\Delta r$ for a pair of SCs. $N_{TA}$ is coded in color. Fit
shows interaction via repulsive potential and obtained interaction strength
$\beta V_{0}$ is shown in (f) as a function of $N_{TA}$. Multiple independent
long simulation runs of these systems for each $N_{TA}$ are used to determine
the interaction parameters.

Assuming equal division rates $k_{d}$ for SC and TA cells and neglecting TA
cell loss on average $2^{N_{TA}}$ growing cells are present at a time.
Further, $2^{N_{TA}}$ TD cells are produced in each cycle. Thus, each stem
cell division can produce at most $2^{N_{TA}+1}$ progeny. However, due to the
stochastic loss rate $k_{a}$ of TA and TD cells in our model (which can occur
either due to apoptosis in monolayers, or out-of-plane delamination in
multilayered tissues such as the epidermis) we arrive at a lower effective
number $\langle N_{p}\rangle<2^{N_{TA}+1}$(see Fig. 2). On average, the stem
cell is thus surrounded by a circular arrangement of its offspring, with a
characteristic length $L_{p}=\sqrt{\langle N_{p}\rangle/\rho\pi}$, where
$\rho$ denotes the cell density (see Fig. 2). Note that $k_{d}$ is not fixed
in our simulation model, but controlled by pressure. An increase of TA cell
generations can thus decreases the proliferation rate and effective cell
number further (see App.C Fig. 7).

However, the average hides an important aspect of the dynamics. While in some
instances the SC is indeed in the center of the cell mass, it is also often
found at the boundary due to the stochastic nature of the orientation of
asymmetric divisions. When the SC is located at the boundary, the divisions of
the TA cells result in an effective propulsion of the SC and persistent
directed motion (App. D Fig. 8).

Correspondingly, the mean squared displacement (MSD) of the stem cell shows a
regime of active motion up to about the lifetime of one generation, and a
regime of activity enhanced long-term diffusion thereafter (see Fig. 4 and
App. D Fig. 9). In analogy to active Brownian particles (ABPs) or Run’n’Tumble
particles, the MSD of stem cells can be described in terms of a translation
diffusion coefficient $D_{t}$, a characteristic rotational time $\tau_{R}$,
which confirms decorrelation of the SC motion after about one cell generation,
and a propulsion velocity $v_{0}$ (see App. D.1). The fitted values for
$\tau_{R}$ are about one half to one apoptosis times and correspond to the
decorrelation of directed motion. For the propulsion velocity we find a
maximum for $N_{TA}=6$, indicating an optimal population size for sped-up
motion.

We quantify the active outbursts of SC motion by the displacement distribution
of the stem cell, i. e. the probability that the SC undergoes a given
displacement within a fraction of a generation (_here_ $0.2/k_{a}$ see App. D
Fig. 11). Heavy tails in the distribution display the persistent motion during
outburst. When a stem cell is displaced from the center of mass of the cell
population the heaviness of the tail increases strongly (see App. D Fig. 11).

##  4 Stem cell interaction

In order to derive an effective pair interaction, we initialize the simulation
with two stem cells in a large periodic box. Each stem cell generates its own
progeny populations, which form an aggregate via which stem cells effectively
repel each other (see Fig. 2). Loss of cells in the region between stem cells
can lead to two separate cell populations, both behaving as described in the
previous section. In simulations with periodic boundary conditions, they
eventually collide again and SC rejection can be observed yet again. We
quantify the stem cell repulsion by measuring the pair correlation function

$g(r)=\left\langle\sum_{m,n\neq
m}^{N_{SC}}\Theta_{H}\left(r-r_{mn}(t_{i})\right)\Theta_{H}\left(r_{mn}(t_{i})-(r+\delta
r)\right)/\Omega(r)\right\rangle$

in an annulus ranging from $r$ to $r+\delta r$, where $\Theta_{H}(x)$ is the
Heaviside function and
$r_{mn}(t_{i})=\allowbreak\left|\boldsymbol{r}_{m}(t_{i})-\boldsymbol{r}_{n}(t_{i})\right|$
the stem cell pair distance, measured applying minimum image convention. The
average combines data from different times and independent runs to improve
statistics and the sum runs over pairs of SC. In order to normalize $g(r)$,
and to take into account the radial bin size, we divide by a geometrical
scaling constant $\Omega(r)=\allowbreak
L_{XY}^{2}/\allowbreak[4(N_{SC}-1)\allowbreak((r+\delta r)^{2}-r^{2})]$. The
probability to find two SCs at short distances is strongly reduced, and
reaches a plateau at larger distances, where the cell clouds are out of
contact again (see Fig. 2).

Next we utilize an analogy to thermal equilibrium to derive an effective
repulsive pair potential $V(r)=V_{0}\exp\left\\{-r/L_{rep}\right\\}$ of very
soft colloids. The repulsion length $L_{rep}$ is set to the characteristic
length $L_{p}$ of progeny around the SC. In equilibrium, the pair correlation
function follows from Boltzmann statistics

$G(r)=G_{0}\exp\left\\{-\beta V_{0}\exp\left\\{-r/L_{p}\right\\}\right\\},$
(1)

with $\beta=1/k_{B}T_{eff}$ and effective temperature $T_{eff}$. The prefactor
$G_{0}$ is obtained from the normalization condition $\int G(x)dx=1$. We
determine the effective potential strength $\beta V_{0}$ by fitting $G(r)$ to
$g(r)$ and find good agreement of this simple model with our data. Results are
shown in Figure 2+2. The decrease of $\beta V_{0}$ with $N_{TA}$ corresponds
to a softening of the potential, which could be due to an increasingly
asymmetric distribution of TA cells.

##  5 Confluent tissue maintained by stem cells

We perform tissue simulations for different $N_{TA}$, and vary the number of
stem cells $N_{SC}$ in the simulation such that the average total cell density
is approximately the same for different simulations. We compare these tissue
simulation with Brownian dynamics (BD) simulations of thermal colloid
particles interacting via the proposed effective repulsion. The interaction
parameters, obtained in section 4, are used without further adaptions.
Snapshots for the mechanical model and the thermal colloid model are shown in
Figure 3 and 3, respectively.

Clearly, SCs are well separated, and homogeneously distributed over the
system. To quantify this separation, we measure again the pair correlation
function for both systems (see Fig. 3). The pair correlation function reveals
a strong depletion of stem cells from the vicinity of another on length scales
set by the amount of progenitors they generate. Interestingly, we find that
the colloid system (dotted lines) reproduces well the results from the tissue
simulations (solid lines) even in confluent tissues. Note that this requires a
rather soft colloidal interaction on length scales much larger than the (stem)
cell size.

A\begin{overpic}[scale={0.4}]{Tissue_Snapshot_SCcircles.png}\put(0.0,90.0){\captiontext*{}}\end{overpic}

B\begin{overpic}[scale={0.4}]{BD_Snapshot_SCcircles.png}\put(-2.0,90.0){\captiontext*{}}\end{overpic}

C\begin{overpic}[scale={0.4}]{TissueRadialDistribution.eps}\put(-1.0,65.0){\captiontext*{}}\end{overpic}

A\begin{overpic}[scale={0.4}]{Tissue-StemCell-
ClusterAnalysis_CDF.eps}\put(0.0,72.0){\captiontext*{}}\end{overpic}

Figure 3:  SC separation in bulk tissue is explained by equilibrium model with
simple repulsive interaction. (a) Snapshot of the tissue simulation for
$N_{TA}=7$ and (b) respective thermal colloid simulation. The characteristic
progeny length scale is displayed by circles. (c) The SC and particle pair
correlation functions, measured from simulations of the two-particle growth
model (solid lines) and thermal colloids, interacting with the assumed
repulsive interaction (dotted lines), are shown for color coded $N_{TA}$. (d)
Cluster size distribution function for SCs in tissues (dots, dashed line) and
thermal colloids (cross, dotted line) for color coded $N_{TA}$. Distance
threshold for clustering is set to $d=2L_{p}$, much larger than the cell size,
in order to see any clustering at all. The number of SCs was chosen such that
the cell density in all simulations is $\sim 1.58$ cells/$R_{PP}$, and the
number of thermal colloids equals the number of SCs. Interaction parameters of
colloids were extracted from single and two SC simulations. Note the
astonishing agreement between both systems.

We further quantify the stem cell separation with a cluster analysis. We
consider that SCs which are found at a distance smaller than $d$ belong to the
same cluster, and calculate the cluster size distribution function
$\mathcal{N}(n)=p(n)\cdot n$, where $p(n)$ is the probability of a cluster of
$n$ SCs, and average cluster size $\langle n\rangle$. If we take the distance
threshold $d$ equal to the cells interaction range $R_{pp}$, as commonly done
for cluster analysis, we identify close to no clusters at all
($\mathcal{N}(n=2)=0.024$ and $\mathcal{N}(n=3)=0.00014$, for $N_{TA}=5$). By
choosing a larger cutoff $d=2L_{p}$ on the basis of their progeny, we find
some small clusters, but still not any large clusters (see Fig. 3 and App. E
Fig. 13), a feature indicative of SC separation, which is also reproduced in
the colloidal systems.

A\begin{overpic}[scale={0.4}]{SCdynamics-MSD-Single-Confluent-
BD.eps}\put(0.0,72.0){\captiontext*{}}\end{overpic}

B\begin{overpic}[scale={0.4}]{SCdynamics-structure-
factor-Q.eps}\put(0.0,74.0){\captiontext*{}}\end{overpic}

F\begin{overpic}[scale={0.4}]{Structure_factor_relaxation_doublecellsize.eps}\put(-4.0,75.0){\captiontext*{}}\end{overpic}

D\begin{overpic}[scale={0.4}]{Structure_factor_relaxation_interaction.eps}\put(0.0,75.0){\captiontext*{}}\end{overpic}

Figure 4:  Dynamics of SC and BD simulations. (a) Mean squared displacement
(MSD) for isolated SC (green) and SCs in tissue (blue), both with $N_{TA}=4$,
and respective thermal colloid system (yellow). In the tissue model, we find
superdiffusive SCs. Their displacement gets reduced with increased cell
density in confluent tissues. Time scale of thermal colloid system was aligned
with the short term diffusion of SC in tissues. (b) Dynamic structure factor
(DSF) as function of wavevector for $t=0$ and $t=0.02k_{a}$. Relaxation of the
DSF as function of time for (c) small structures ($q=3.08/R_{PP}$ or
$\lambda_{q}=2.03R_{PP}$, i. e. twice the cell diameter) and (d) large
structures ($q=0.69/R_{PP}$ or $\lambda_{q}=9.11R_{PP}$, i. e. comparable to
"weakly" interacting stem cells). Small structures relax at comparable
timescale, whereas large structures relax faster in the tissue model due to
the superdiffusive behavior of stem cells.

To obtain more insight in the structural dynamics of both systems, we
additionally calculate the MSD of SC, and the intermediate scattering function

$S(q,t)=\left\langle\sum_{m=1}^{N_{SC}}\sum_{n=1}^{N_{SC}}e^{i\boldsymbol{q}\boldsymbol{r}_{m}(t)}e^{-i\boldsymbol{q}\boldsymbol{r}_{n}(0)}\right\rangle,$

where $\langle\cdot\rangle$ denotes a time average and we calculate the over
$S(q,t)$ with wave vector $q=|\boldsymbol{q}|$ (see Fig. 4).

The MSD of SCs in tissues shows the same superdiffusive behaviour as found for
an isolated SC in free space (see Fig. 4). However, the displacement is
reduced in confluent tissue due to the increased cell density. We find, that
this affects mainly the translational diffusion coefficient and propulsion
velocity, while the reorientation remains in the same order of magnitude as
for isolated SC (see App.D Fig. 10). The thermal colloid system on the other
hand shows normal diffusion, as should be. Because the Boltzmann description
does not provide any timescales, we use the short-term diffusive limit of
colloids and SCs in tissues to fix the timescale. The short-term diffusion
coefficient can also be calculated from the relaxation of the intermediate
scattering function as $D_{eff}=-\log\\{S(q,t)/S(q,0)\\}/(q^{2}t)$ in the
short time limit. After aligning the time scale of the thermal colloid system,
we find that this effective diffusion is approximately the same for both
systems (see App. E Fig. 14). For large $q$, $D_{eff}$ is constant, while for
small $q$, finite-size effects of the simulation area come into play.

For zero time lag ($S(q,\Delta t=0)$), we find a plateau, corresponding to the
autocorrelation term ($n=m$). With a finite time lag, the scattering function
starts to decay with $q$ (see Fig. 4 for $\Delta t=0.02k_{a}$). For tissues
and colloids, we show the relaxation of small and large structures in Fig. 4
and Fig. 4, respectively. Because we match the timescale of both systems for
the short-term diffusion, we find that small structures relax at similar
times. Still, SCs in tissue show faster relaxation compared to the BD system.
For larger structures (small $q$), the SCs in tissue relax significantly
faster compared to the thermal colloids due to their superdiffusive motion.
The enhanced relaxation with increased structure size holds for all $q$ as can
be seen from the half life relaxation time (see App. E Fig. 14).

A \begin{overpic}[scale={0.4}]{box_30_SC_3_RadialDistributionFunction.eps}
\put(0.0,64.0){\captiontext*{}} \end{overpic}

B
\begin{overpic}[scale={0.45}]{ConfluentTissue_SC_3_TA_10_snapshot_50.00tapo.png}
\put(10.0,60.0){\captiontext*{}} \end{overpic}

Figure 5:  Stem Cell repulsion qualitatively remains in three dimensions. (a)
Pair correlation function for stem cells in three dimensions shows reduced
probability to find SCs at short distances, and plateau at longer distances.
(b) Snapshot for a three dimensional tissue with $N_{TA}=10$ after 50
generations. All simulations are performed with 27 SCs in a cubic box with box
length $30~{}R_{PP}$.

##  6 Stem cell repulsion in three dimensions

Many biological tissues are not restricted to two dimensions, but extend to
the third dimension, which is especially the case in cancer, where cells grow
out of plane and form spheroidal tumors. We show that the effect of repulsion
by progeny is not limited to two dimensions and can also be found in tissue
simulations in three dimensions (see Fig. 5). Qualitatively, the pair
correlation function shows the same repulsion of SC at a distance
corresponding to the volume of their progenitors. Note that here the geometric
normalization of the pair correlation function is the volume of a sphere
instead of the area of a circle around the stem cell.

##  7 Independence of homeostatic pressure controlled divisions

Stem cells interacting effectively via their (growing) progeny is a robust
effect in tissues with mechanical cell-cell interactions. We show this by
replacing the division after reaching a size threshold, which introduces a
feedback on the homeostatic pressure, by stochastic divisions with fixed rate
$k_{d}$. This rate is set to the measured division rate in the size threshold
division model for an exemplary tissue with $N_{TA}=6$. Details on the
implementation can be found in Appendix F.

The distribution of TA cells around the SC remains almost unchanged (see Fig.
6), thus, also the SC pair correlation function of both models remains the
same (see Fig.6). The SC dynamics, measured by the MSD, remains unchanged as
well (see App. F Fig. 16). Thus, the mechanical cell-cell interactions in the
model and outflux of cells from the SC are generating an effective SC
repulsion, but not the feedback of growth and division on pressure in the
tissue.

B\begin{overpic}[scale={0.4}]{StochsticDivisions_snapshotstoch.png}\put(-5.0,88.0){\captiontext*{}}\end{overpic}

D\begin{overpic}[scale={0.4}]{StochsticDivisions_RDF.eps}\put(0.0,70.0){\captiontext*{}}\end{overpic}

Figure 6: Relevance of homeostatic pressure growth control. (a) Snapshot for
deterministic TA cell differentiation ($N_{TA}=6$) in the model with
stochastic division mechanism. (b) Pair correlation function for simulation
model with division after reaching a size threshold (red) and stochastic
divisions (blue).

##  8 Discussion and conclusion

In this work, we have studied the dynamics of stem cells in self-renewing
tissues driven by mechanical interactions. Due to the constant outward flow of
progeny generated by stem cells, we find that SCs effectively repel each other
and spontaneously organize into dispersed structures with a corona of their
progeny. This provides an updated and more biophysical view of the classical
concept of proliferating units in tissues, in which size and shape of self-
renewing domains are flexible, but still organized by stem cells which are
robustly present within each unit [18].

Interestingly, we find that the interaction between stem cells can be largely
described in a simplified fashion with help of an effective equilibrium
approach. The assumption purely repulsive interaction and determination of the
interaction length from the spatial requirements for the progeny produced by
each stem cell division allows us to characterize SCs like very soft colloids,
which interact via this effective repulsive interaction – much larger than the
cell size. Further, we have shown that this simplified model reproduces the
structural results of the two-particle growth model in confluent tissues. The
pair correlation function of the thermal colloids is remarkably similar to
that of actively growing cells. The dynamics however are different. We show
evidence of an active and persistent self-propulsion force acting on SCs due
to the surrounding growing cells. We systematically characterize these, in
particular in terms of the statistics and active nature of their persistent
random walks. As a result of this proliferation-driven movements, the
confluent tissue shows faster relaxation characteristics on long time and
length scales. The only detectable difference between the actively growing
tissue and the thermal colloid system was found in the dynamic structure
factor. It shows a speedup of large structure relaxations in the tissue model
due to its active nature, and thus confirms that the active system can be
mapped surprisingly well to an thermal equilibrium system.

The derivation of an effective interaction using nearest neighbor distribution
functions and mapping it on particle systems may lack sensitivity to certain
characteristics of the interaction, especially with increasing density [19,
20]. However, they are still able to capture the main influences of the
interaction. Further, replacing the passive Brownian particle system by an
active one, is very likely to fail as the continuously decrease of activity
enhanced motion cannot be derived from first principles of infinitely diluted
systems [21].

Small populations of SCs dispersed within the tissue were also shown to be
particularly susceptible to tumor formation in mouse epidermis [22] as well as
to have large contribution upon wound healing [23], which could also be
investigated via our simulation framework. Furthermore, spermatogenic stem
cells are also found to be well separated in the planar basal layer in mouse
testis, and show a striking similarity to the stem cell arrangement found in
our simulations [24]. Other models of the formation of epithelial protrusions
and location of stem cells on the tips of these have been studied previously
using a particle-based model [25]. By adding substrate deformability and cell
type dependent adhesion forces, the model displayed a stable stratified
epithelium with stem cells located at protrusions. It will be interesting to
combine these models.

The homogeneous distribution of stem cells arising from our purely mechanical
feedback model may complement previous studies [3, 4], in which this was
achieved with help of a reversible differentiation scheme or biochemical
feedback mechanisms based on diffusible signals. However, in reality, a
combination of several of these mechanisms will likely be relevant, also
depending on the tissue type. In particular, stochastic cell fate models may
be of great interest in the future and could facilitate the description of the
plasticity of cells during tissue repair [26, 27]. Also, there exists evidence
that more than one SC type maintains some tissues [28, 29] or that cell fates
are reversible processes [30, 3], which could be readily integrated in our
simulation framework. In the presence of local niches, it has also been shown
that the ratio of proliferation to random cell movements determines the
effective number of stem cells via a mechanism of stochastic competition for
space, although this has not been modelled via explicit mechanical simulations
[31, 32]. Furthermore, it would be interesting to see if a purely mechanical
model is able to capture the alternation from a quiescent long-lived SC to
produce rapidly dividing short-lived progeny during wound repair [33, 34, 35].

Finally, our results might contribute to a better understanding of the growth
and dynamics of cancer. In the cancer stem cell (CSC) hypothesis, it is
assumed that malignant tissue follows the same underlying stem cell dynamics.
Therapies, which target fast growing tissue, fail, since the stem cells evade.
A better understanding of the dynamics of CSCs could improve their detection
and the design of therapies. Our results would suggest, that CSCs are isolated
and well spread throughout the tumor. To further tackle these questions, it
will be necessary to combine our model, with models of cellular competition
[5, 7, 8, 9, 11], and taking the cell cycle explicitly into account might be
essential [36, 37]

## Acknowledgements

#### Funding information

JE and JK gratefully acknowledge financial support from the Initiative and
Networking Fund (IVF) via the grant number ERC-RA-004. Simulations were
performed with computing resources granted by RWTH Aachen University under
project ’rwth0475’.

#### Data accessibility

Source code, simulation data, and analysis scripts will be made publicly
available on zenodo.org (10.5281/zenodo.8410957) after acceptance of the
manuscript.

##  Appendix A Supplementary animations of tissue simulations

Supplementary animations of tissue simulations can be found next to the source
code, simulation data, and analysis scripts on _zenodo.org_
(10.5281/zenodo.8410957). Uploaded animations are:

  * •

_S1 Video_ Confluent tissue. Animation of confluent tissue shown in Fig. 1.

  * •

_S2 Video_ Isolated stem cell. Animation of isolated SC performing random walk
with long persistent segments (see Fig. 2 and sec. D).

  * •

_S3 Video_ Stem cell interaction. Animation of two SCs interacting with each
other via their progeny (see Fig. 2).

  * •

_S4 Video_ Interaction in Tissue and BD simulations. Animation of confluent
tissue and thermal colloids (see Fig. 3 and Fig. 3).

##  Appendix B Simulation model

The described SC model is implemented in the two-particle growth (2PG) model
of Refs. [17, 7] and has been adapted in Refs. [9, 10, 11].

Each cell is described by two particles, which repel each other via an active
growth force

$F_{ij}^{(G)}=\frac{G}{\left(r_{ij}-r_{0}\right)^{2}}\hat{\boldsymbol{r}}_{ij}$
(2)

with unit vector $\hat{\boldsymbol{r}}_{ij}$, distance $r_{ij}$ between the
two particles, a constant $r_{0}$, and growth strength $G$. For non-growing
cells, like TD cells, the growth strength is set to zero, $G_{TD}=0$, and for
the sake of simplicity SC and TA cells have the same growth force,
$G_{SC}=G_{TA}$.

To prevent overlap of cells, particles of different cells interact via a soft
repulsive volume-exclusion force

$\boldsymbol{F}_{ik}^{(V)}=f_{0}\left(\frac{R_{PP}^{5}}{r_{ik}^{5}}-1\right)\hat{\boldsymbol{r}}_{ij}\text{,
for }r_{ik}<R_{PP}$ (3)

with exclusion strength $f_{0}$ and interaction length $R_{PP}$, which sets
the length scale of our simulations. Further, cells in contact interact via an
attractive adhesion force of the form

$\boldsymbol{F}_{ik}^{(V)}=-f_{1}\hat{\boldsymbol{r}}_{ik}\text{, for
}r_{ik}<R_{PP}$ (4)

with adhesion strength $f_{1}$.

We employ a dissipative particle dynamics-type thermostat [38], with an
effective temperature $T$, to account for energy dissipation

$F^{(D)}_{ij}=\gamma\cdot\left(1-\frac{dr}{rt}\right)^{2}\cdot\left(\boldsymbol{v}_{ij}\cdot\boldsymbol{r}_{ij}\right)$
(5)

and random fluctuations

$F^{(R)}_{ij}=\sigma\xi\cdot\hat{\boldsymbol{r}}_{ij}$ (6)

where $\gamma$ and $\sigma$ are related to fulfill the fluctuation-dissipation
theorem [39]. Also, background dissipation is taken into account as

$\boldsymbol{F}^{B}_{i}=-\gamma_{b}\boldsymbol{v}_{i}.$ (7)

A self-consistent velocity-Verlet algorithm [40] is implemented to integrate
the equations of motion, and all simulations were performed with periodic
boundary conditions.

Cell division is performed, when the two particles of one cell are separated
by a critical threshold $R_{ct}$. A new particle is randomly placed near each
original particle, and the two new particle pairs form the two daughter cells.
Differentiation is implemented at the time of division as described above, and
the number of divisions of TA cells is tracked by an internal variable.

TA and TD cells are stochastically removed from the simulation with a finite
apoptosis rate $k_{a}$. For simplicity, the same rate is chosen for different
cell types. Time is measured in terms of the inverse apoptosis rate, and
referred to as "generation".

The standard parameter set for our simulations are given in Tab. 1.

Table 1: Simulation parameters of the standard tissue. Parameter | Symbol | Value
---|---|---
Time step | $\Delta t$ | 0.001
Pair potential interaction range | $R_{PP}$ | 1
Cell expansion constant | $r_{0}$ | 1
Division distance threshold | $R_{ct}$ | 0.8
New cell particle initial distance | $R_{d}$ | $10^{-5}$
Growth-force strength | $G$ | 50
Mass | $m$ | 1
Intra-cell dissipation coefficient | $\gamma_{c}$ | 100
Inter-cell dissipation coefficient | $\gamma_{t}$ | 50
Background dissipation coefficient | $\gamma_{c}$ | 0.1
Apoptosis rate | $k_{a}$ | 0.01
Noise intensity | $k_{B}T$ | 0.1
Repulsive cell-cell potential coefficient | $f_{0}$ | 2.39566
Attractive cell-cell potential coefficient | $f_{1}$ | 3.0

##  Appendix C Division rate of isolated SC

Assuming equal division rates $k_{d}$ for all types of proliferating cells, we
can derive an average cell number from number balance for each cell type

$\displaystyle n_{SC}$ $\displaystyle=\langle n_{SC}\rangle=const.$ (8)
$\displaystyle\langle n_{i}\rangle$
$\displaystyle=2^{i-1}\left(\frac{k_{d}}{k_{d}+k_{a}}\right)^{i}n_{SC}$
$\displaystyle\langle n_{TD}\rangle$
$\displaystyle=\left(\frac{2k_{d}}{k_{d}+k_{a}}\right)^{N_{TA}}\frac{k_{d}}{k_{a}}n_{SC},$

where $i$ denotes the TA cell cycle. $k_{d}$ is obtained from fitting these
equations to the cell numbers obtained from simulations of isolated SCs. With
increasing $N_{TA}$, the division rate decreases (see Fig. 7).

Figure 7: Division rate $k_{d}$ as function of $N_{TA}$ for proliferating
cells in isolated SC simulations.

##  Appendix D Stem cell motion

Here, we give some additional insight in the SC motion of isolated SCs and SCs
in tissues, which are also discussed in the main text.

An exemplary single stem cell snapshot for $N_{TA}=4$ is shown, where the SC
trajectory over ten generations illustrates the motion of the SC, and
accompanied The full trajectory of the SC and center of mass of the respective
single stem cell simulation along both directions of the simulation plane (see
Fig. 8).

A \begin{overpic}[scale={0.4}]{Snapshot_singleSC_trajectory.png}
\put(0.0,95.0){\captiontext*{}} \end{overpic}

B \begin{overpic}[scale={0.4}]{singleSC-trajectory-full.eps}
\put(0.0,65.0){\captiontext*{}} \end{overpic}

Figure 8:  (a) Snapshot of cell population generated by a single SC with
trajectory of SC over the ten cell generations (red line). (b) Trajectory of
isolated SC (solid line) and center of mass (CoM) of cell population (dotted
line) along $\hat{e}_{x}$ (red) and $\hat{e}_{y}$ (blue). Inset shows time
close-up to highlight how the center of mass follows the stem cell motion.
Both, SC and CoM trajectory perform a random walk with surprisingly long
persistent segments, a reminiscent of Run’n’Tumble-like motion. Simulation
with $N_{TA}=4$ is shown.

###  D.1 Mean squared displacement of isolated stem cells and stem cells in
tissues

For active Brownian particles (ABPs) in two-dimensions the mean squared
displacement $MSD$ is given as [41, 42]

$MSD=4~{}D_{t}\Delta t+\frac{v_{0}^{2}\tau_{R}^{2}}{2}\left[\frac{2\Delta
t}{\tau_{R}}+e^{-2\Delta t/\tau_{R}}-1\right],$ (9)

with translational diffusion coefficient $D_{t}$, characteristic time scale
for rotational diffusion $\tau_{R}$, and propulsion velocity $v_{0}$. In the
short-time limit $dt\ll\tau_{R}$ the effective diffusion is given by
$MSD_{dt\ll\tau_{R}}\approx 4D_{t}dt$. For long-time diffusion the effective
diffusion is enhanced and the MSD found to be
$MSD\approx(4D_{t}+v_{0}^{2}\tau_{R})~{}dt$, again linear in time. Between
these two regimes, ballistic motion is observed, and non-linear terms enters
the MSD, $MSD\approx 4D_{t}dt+v_{0}^{2}dt^{2}$.

In the tissue simulations, the stem cells are not ABPs, but move actively due
to propulsion by progeny. However, we can quantify the dynamics of the stem
cell motion in terms of $D_{t}$, $\tau_{R}$, and $v_{0}$ to gain a deeper
understanding. Also different models of self-propelling particles, like ABPs
and Run’n’Tumble particles have been shown to be equivalent on larger scales
[43, 44], so that we can describe them using the same theory, without making
strong assumptions on the underlying propulsion mechanisms.

A\begin{overpic}[scale={0.4}]{StemCellMobility_SingleSC_MSD_NTA_fitted.eps}\put(0.0,60.0){\captiontext*{}}\end{overpic}

B\begin{overpic}[scale={0.4}]{StemCellMobility_SingleSC_MSD_NTA_Dtrans.eps}\put(-2.0,72.0){\captiontext*{}}\end{overpic}

A\begin{overpic}[scale={0.4}]{StemCellMobility_SingleSC_MSD_NTA_tauR.eps}\put(-5.0,72.0){\captiontext*{}}\end{overpic}

B\begin{overpic}[scale={0.4}]{StemCellMobility_SingleSC_MSD_NTA_v0.eps}\put(0.0,70.0){\captiontext*{}}\end{overpic}

Figure 9:  (a) Mean squared displacement for isolated SCs with different
$N_{TA}$ (color coded) and fit of ABP model (red dotted lines). (b)
Translational diffusion coefficient $D_{t}$, (c) characteristic rotational
diffusion time $\tau_{R}$, and (d) propulsion velocity $v_{0}$ obtained from
fit function in (a).

From the MSDs for simulations of SCs with varying number of TA cell cycles
$N_{TA}$ (see Fig. 9), one can see that the displacement does vary little with
a change of $N_{TA}$. In particular the short-time diffusion decreases with
$N_{TA}$, which could be due to blocking of cell motion at higher density.
This resembles the decrease of the translational diffusion coefficient $D_{t}$
with $N_{TA}$. The long-time diffusion first increases with $N_{T}A$, but then
starts to decrease again. Possibly, in increasing number of TA cells around
the stem cell initially increases the propulsion force on the stem cell, but
saturates at higher number of progeny, which form a more dense halo around the
stem cell for more TA cell cycles. The non-monotonic fit parameter for the
self propulsion velocity, with maximum at $N_{TA}=6$, and increasing
characteristic rotational diffusion time $\tau_{R}$ show, that the decrease in
propulsion velocity does cause the effective decrease of the long-time
diffusion with $N_{TA}>6$. The characteristic rotational diffusion time
$\tau_{R}$ is found to be in the range from half to more than one apoptosis
time, which corresponds to a loss of orientation after times larger than one
cell generation.

A\begin{overpic}[scale={0.4}]{StemCellMobility_tissueSC_MSD_NTA_fitted}\put(0.0,57.0){\captiontext*{}}\end{overpic}

B\begin{overpic}[scale={0.4}]{StemCellMobility_tissueSC_MSD_NTA_Dtrans.eps}\put(0.0,70.0){\captiontext*{}}\end{overpic}

A\begin{overpic}[scale={0.4}]{StemCellMobility_tissueSC_MSD_NTA_tauR.eps}\put(-5.0,70.0){\captiontext*{}}\end{overpic}

B
\begin{overpic}[scale={0.4}]{StemCellMobility_tissueSC_MSD_NTA_v0.eps}\put(-5.0,75.0){\captiontext*{}}\end{overpic}

Figure 10:  (a) Mean squared displacement for SCs in tissue with different
cell density (color coded) controlled by the stem cell number, with
$N_{TA}=6$. (b) Translational diffusion coefficient $D_{t}$, (c)
characteristic rotational diffusion time $\tau_{R}$, and (d) propulsion
velocity $v_{0}$ obtained from fit function in (a).

In tissues, we still can apply the MSD theory to fit to the SC displacement
(see Fig. 10). With increasing cell density, controlled by an increase in SC
numbers, we find that the translational diffusion coefficient $D_{t}$ and
propulsion velocity $v_{0}$ decrease, and the SC motion freezes out. The
rotational diffusion time is not as strongly effected by the cell density, and
approximately half the value compared to the case of an isolated SC, since it
is mainly affected by the apoptosis rate.

###  D.2 Displacement distribution

To quantify the persistent motion further, we calculate the displacement
distribution function and find heavy tails. The heaviness of the tails
increases, when the SC is located further away from the center of mass in the
boundary region of the cell population (see Fig. 11).

A
\begin{overpic}[scale={0.4}]{DisplacementDistribution_NTA_F1_3_DT_20_dr_sc_log.eps}
\put(0.0,62.0){\captiontext*{}} \end{overpic}

B
\begin{overpic}[scale={0.4}]{DisplacementDistribution_NTA_F1_3_DT_20_dr_sc_toCoM_log.eps}
\put(0.0,60.0){\captiontext*{}} \end{overpic}

Figure 11:  (a) Displacement distribution for an isolated stem cell with
different $N_{TA}$ (green color gradient) over normalized displacement. Red
dotted line displays normal displacement (Gaussian). Heavy tails emerge from
long consistent segments of random walk. (b) Displacement distribution as in
(a) for isolated SC with $N_{TA}=5$, where color gradient (grey) encodes SC-
to-CoM distance. With increasing distance heaviness of tails increases. All
displacements are calculated for a time step of $\Delta t~{}k_{a}=0.2$
generations, combining data of multiple long simulation runs.

##  Appendix E Additional results for simulations of stem cells in tissues

###  E.1 Cell densities in confluent tissues

In confluent tissue simulations, we chose the SC number and $N_{TA}$ such that
an uninterrupted tissue evolves. Still, the TD cell density around SCs remains
very similar to the case of isolated SCs. The TD cells are located closer to
the SCs and form a more dense conglomerate. Exemplary cell densities for
$N_{TA}=6$ are shown in Fig. 12.

The total cell density of tissues with different $N_{TA}$ is controlled by the
number of SCs in the simulation, and chosen such that they are all very
similar (see Fig. 12).

A \begin{overpic}[scale={0.4}]{CellColonyDistribution_density.eps}
\put(-5.0,75.0){\captiontext*{}} \end{overpic}

A \begin{overpic}[scale={0.4}]{TissueDensity_NTA.eps}
\put(0.0,70.0){\captiontext*{}} \end{overpic}

Figure 12: (a) Cell density of TA (solid line) and all (dashed line) cells
around the SC for simulations of an isolated SC (blue) and SCs in tissue
(green). In tissue, distance is measured to the closest SC. (b) Average cell
density of tissue simulations as function of $N_{SC}$, the fixed number of SC,
$N_{TA}$ obtained from average cell number over box area.

###  E.2 Cluster analysis

To quantify the clustering of SCs in tissues, we calculate the cluster size
distribution function is given by

$\mathcal{N}(n)=\frac{1}{N_{SC}}n~{}p(n)$ (10)

and represents the fraction of SCs belonging to a cluster of size $n$, where
$p(n)$ is the number of clusters of size $n$. The distribution is normalized
such that $\sum_{i=1}^{N_{S}C}=1$. The average cluster size is given as

$\langle n\rangle=\frac{\sum_{n}np(n)}{\sum_{n}p(n)}$ (11)

Cells in a distance less than the characteristic TA cell distance estimated
from cell numbers belong to the same cluster. Results are shown in Fig. 13.

A \begin{overpic}[scale={0.4}]{Tissue-StemCell-ClusterAnalysis_CDF.eps}
\put(0.0,70.0){\captiontext*{}} \end{overpic}

B \begin{overpic}[scale={0.4}]{Tissue-StemCell-ClusterAnalysis_CDF_avg.eps}
\put(0.0,70.0){\captiontext*{}} \end{overpic}

Figure 13:  (a) Same as Fig. 3 (b) Average cluster size as function of
$N_{TA}$ remains clearly below two SC.

###  E.3 Effective diffusion coefficient and relaxation

The timescale for the thermal colloid simulation in the main text was adjusted
by matching the short term diffusion obtained from the MSD of SCs in the
tissue simulations. Here, we calculate the effective diffusion coefficient
$D_{eff}$ from the intermediate structure factor. For the relaxation of
structures it holds

$\frac{S(q,\Delta t)}{S(q,0)}=e^{-D_{eff}(q)~{}\Delta t~{}q^{2}\Delta t},$
(12)

which yields

$D_{eff}(q)=\log\left\\{\frac{S(q,\Delta t)}{S(q,0)}\right\\}/(q^{2}\Delta
t).$ (13)

Now, we can measure the effective diffusion for short time lags (here:
$0.02k_{a}$) and see that both models show the same short term diffusivity, as
should be (see Fig. 14).

The relaxation of structures is obtained from the half-life relaxation time of
$S(q,t)/S(q,0)$. We find that this time decreases with $q$, which means that
smaller structures relax faster. Further, in the tissue model structures relax
faster than in the thermal colloid system and the relaxation speeds up with
structure size (see Fig. 14).

A \begin{overpic}[scale={0.4}]{SCdynamics-structure-factor-
effectiveDiffusionCoefficient.eps} \put(-5.0,77.0){\captiontext*{}}
\end{overpic}

B \begin{overpic}[scale={0.4}]{Structurefactor_halflife_md_q.eps}
\put(0.0,70.0){\captiontext*{}} \end{overpic}

Figure 14:  (a) Effective diffusion coefficient $D_{eff}$ measured from
structure factor relaxation for $\Delta t=0.2k_{a}$. (b) Relaxation half life
time as function of $q$.

##  Appendix F Stochastic division rate model

The observed separation of stem cell in tissues is a generic effect arising
from short range mechanical interactions between pairs of cells. To highlight
this, we replace the size-threshold division (STD) mechanism by a stochastic
division with fixed rate $\tilde{k}_{d}$ (SDR). This turns off the homeostatic
pressure control on divisions. We determine the simulation parameter
$\tilde{k}_{d}$ for the SDR model from measurements of $k_{d}$ in the STD
model. Only the case of tissues with $N_{TA}=5$ is examined.

To implement, that cells do not divide at too small size and cause
computational flaws, we replace the growth force by a damped harmonic
oscillator

$\boldsymbol{F}_{i,j}^{Gharm.}=-k\boldsymbol{x}-\gamma_{k}\dot{\boldsymbol{x}},$
(14)

where $k=100$ is the spring constant, and $\gamma_{k}=20$ the damping
coefficient, chosen in a way to realize critical damping with relaxation time
faster than the average division time. The deflection
$\boldsymbol{x}=\left(|\boldsymbol{r}_{i}-\boldsymbol{r}_{j}|-R_{rest}\right)\hat{\boldsymbol{e}}_{ij}$
is given by the difference between cell size and rest length $R_{rest}$. We
set the rest length of proliferating cells, i.e. SC and TA cells, to
$R_{rest}^{SC,TA}=R_{ct}=0.8R_{PP}$ and for TD cells to
$R_{rest}^{TD}=0.1R_{PP}$.

Despite the fact that we tune the oscillator such that cells grow quickly to a
adequate size, some cells divide at very small size and cause computational
chaos. To prevent this, we introduce a minimum division size, much smaller
than the division size threshold used previously. With $R_{min}=0.15R_{PP}$ we
can ensure stable behavior of the computational model after division.

A \begin{overpic}[scale={0.4}]{StochsticDivisions_CDF_DivApoEvens.eps}
\put(-5.0,70.0){\captiontext*{}} \end{overpic}

B \begin{overpic}[scale={0.4}]{StochsticDivisions_CellsizeDistribution.eps}
\put(0.0,70.0){\captiontext*{}} \end{overpic}

Figure 15: (a) Complementary cumulative distribution function (cCDF) for
apoptosis (blue), division events including prevented divisions due to small
size (orange), and performed divisions (green). Black dotted line highlights
longest time for a prevented division event. Black dashed lines display median
of each cCDF. (b) Cellsize distribution for proliferating (SC and TA cells,
solid line) and TD cells (dotted lines) in the stochastic division model
(blue) and for size threshold divisions (red).

Figure 15 displays the complementary cumulative distribution function (cCDF)
of apoptosis, division, and non-permitted division events, i.e. division
events at too small cell size. Only a small fraction of events at very short
times is permitted and thus, we sufficiently removed the size-threshold
division mechanism. Further, the cell size distribution (see Fig. 15) shows,
that TD cells in both models match each other well, and for growing cells, we
find that the hard division size threshold changes to a blurred tail. The
stochastic division model has equal number of TA cells, both in total and for
each generation, and did produce only slightly more TD cells. However, cell
numbers fluctuate more, which can be to the detriment of stable tissue
homeostasis.

A \begin{overpic}[scale={0.4}]{StochsticDivisions_snapshotstoch.png}
\put(-15.0,85.0){\captiontext*{}} \end{overpic}

B \begin{overpic}[scale={0.4}]{StochsticDivisions_CellDistribution_SC.eps}
\put(-10.0,75.0){\captiontext*{}} \end{overpic}

C \begin{overpic}[scale={0.4}]{StochsticDivisions_RDF.eps}
\put(0.0,70.0){\captiontext*{}} \end{overpic}

D \begin{overpic}[scale={0.4}]{StochsticDivisions_MSD.eps}
\put(0.0,70.0){\captiontext*{}} \end{overpic}

Figure 16:  (a) Snapshot for stochastic division model, (b) TA and TD cell
distribution around nearest neighbor SC in tissue, (c) pair correlation
function, and (d) Mean squared displacement for the respective models. For
comparison, stochastic division model is shown by blue lines together with
results for model with size threshold division for $N_{TA}=5$.

When comparing both models, we find that SCs are still surrounded by their TA
cell cloud. In the stochastic division model the TA cell arrangement around
the stem cells looks more fringed compared to the deterministic model but no
large differences can be found in the measured TA cell distribution (see Fig.
16, 16). Still, for both models, the pair correlation function agree very well
(see Fig. 16) and the activity of SC in the tissue remains almost unchanged
(see Fig. 16).

We conclude that the results discussed in the main text are not a consequence
of homeostatic pressure growth control, but result from the constant cell
outflux generated by the SC. Yet, the pressure feedback on divisions is
beneficial to balance cell numbers in the tissue to ensure homeostasis.

## References

  * [1] P. A. Hall and F. M. Watt,  _Stem cells: the generation and maintenance of cellular diversity_ ,  Development 106(4), 619 (1989),  10.1242/dev.106.4.619.
  * [2] G. Mannino, C. Russo, G. Maugeri, G. Musumeci, N. Vicario, D. Tibullo, R. Giuffrida, R. Parenti and D. Lo Furno,  _Adult stem cell niches for tissue homeostasis_ ,  Journal of Cellular Physiology 237(1), 239 (2022),  10.1002/jcp.30562.
  * [3] P. Greulich and B. D. Simons,  _Dynamic heterogeneity as a strategy of stem cell self-renewal_ ,  Proceedings of the National Academy of Sciences 113(27), 7509 (2016),  10.1073/pnas.1602779113.
  * [4] D. J. Jörg, Y. Kitadate, S. Yoshida and B. D. Simons,  _Competition for Stem Cell Fate Determinants as a Mechanism for Tissue Homeostasis_ ,  Tech. Rep. arXiv:1901.03903, arXiv,  10.48550/arXiv.1901.03903 (2019).
  * [5] M. Basan, T. Risler, J. Joanny, X. Sastre-Garau and J. Prost,  _Homeostatic competition drives tumor growth and metastasis nucleation_ ,  HFSP Journal 3(4), 265 (2009),  10.2976/1.3086732.
  * [6] J. Ranft, M. Basan, J. Elgeti, J.-F. Joanny, J. Prost and F. Julicher,  _Fluidization of tissues by cell division and apoptosis_ ,  Proceedings of the National Academy of Sciences 107(49), 20863 (2010),  10.1073/pnas.1011086107.
  * [7] N. Podewitz, M. Delarue and J. Elgeti,  _Tissue homeostasis: A tensile state_ ,  EPL (Europhysics Letters) 109(5), 58005 (2015),  10.1209/0295-5075/109/58005.
  * [8] N. Podewitz, F. Jülicher, G. Gompper and J. Elgeti,  _Interface dynamics of competing tissues_ ,  New Journal of Physics 18(8), 083020 (2016),  10.1088/1367-2630/18/8/083020.
  * [9] N. Ganai, T. Büscher, G. Gompper and J. Elgeti,  _Mechanics of tissue competition: interfaces stabilize coexistence_ ,  New Journal of Physics 21(6), 063017 (2019),  10.1088/1367-2630/ab2475.
  * [10] T. Büscher, A. L. Diez, G. Gompper and J. Elgeti,  _Instability and fingering of interfaces in growing tissue_ ,  New Journal of Physics 22(8), 083005 (2020),  10.1088/1367-2630/ab9e88.
  * [11] T. Büscher, N. Ganai, G. Gompper and J. Elgeti,  _Tissue evolution: mechanical interplay of adhesion, pressure, and heterogeneity_ ,  New Journal of Physics 22(3), 033048 (2020),  10.1088/1367-2630/ab74a5.
  * [12] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts and P. Walter,  _Epidermis and Its Renewal by Stem Cells_ ,  Molecular Biology of the Cell. 4th edition (2002),  Publisher: Garland Science.
  * [13] P. H. Jones, B. D. Simons and F. M. Watt,  _Sic Transit Gloria: Farewell to the Epidermal Transit Amplifying Cell?_ ,  Cell Stem Cell 1(4), 371 (2007),  10.1016/j.stem.2007.09.014.
  * [14] S. W. Lane, D. A. Williams and F. M. Watt,  _Modulating the stem cell niche for tissue regeneration_ ,  Nature Biotechnology 32(8), 795 (2014),  10.1038/nbt.2978.
  * [15] M.-K. Hayward, J. M. Muncie and V. M. Weaver,  _Tissue mechanics in stem cell fate, development, and cancer_ ,  Developmental Cell 56(13), 1833 (2021),  10.1016/j.devcel.2021.05.011.
  * [16] S. Yoshida,  _Open niche regulation of mouse spermatogenic stem cells_ ,  Development, Growth & Differentiation 60(9), 542 (2018),  10.1111/dgd.12574.
  * [17] M. Basan, J. Prost, J.-F. Joanny and J. Elgeti,  _Dissipative particle dynamics simulations for biological tissues: rheology and competition_ ,  Physical Biology 8(2), 026014 (2011),  10.1088/1478-3975/8/2/026014.
  * [18] L. R. Strachan and R. Ghadially,  _Tiers of Clonal Organization in the Epidermis: The Epidermal Proliferation Unit Revisited_ ,  Stem Cell Reviews 4(3), 149 (2008),  10.1007/s12015-008-9020-6.
  * [19] P. G. Bolhuis and A. A. Louis,  _How To Derive and Parameterize Effective Potentials in Colloid-Polymer Mixtures_ ,  Macromolecules 35(5), 1860 (2002),  10.1021/ma010888r.
  * [20] H. Wang, F. H. Stillinger and S. Torquato,  _Sensitivity of pair statistics on pair potentials in many-body systems_ ,  The Journal of Chemical Physics 153(12), 124106 (2020),  10.1063/5.0021475.
  * [21] J. Decayeux, J. Fries, V. Dahirel, M. Jardat and P. Illien,  _Isotropic active colloids: explicit vs. implicit descriptions of propulsion mechanisms_ ,  10.48550/arXiv.2309.08455 (2023).
  * [22] A. Sánchez-Danés, E. Hannezo, J.-C. Larsimont, M. Liagre, K. K. Youssef, B. D. Simons and C. Blanpain,  _Defining the clonal dynamics leading to mouse skin tumour initiation_ ,  Nature 536(7616), 298 (2016),  10.1038/nature19069.
  * [23] M. Aragona, S. Dekoninck, S. Rulands, S. Lenglez, G. Mascré, B. D. Simons and C. Blanpain,  _Defining stem cell dynamics and migration during wound healing in mouse skin epidermis_ ,  Nature Communications 8(1), 14684 (2017),  10.1038/ncomms14684.
  * [24] K. Hara, T. Nakagawa, H. Enomoto, M. Suzuki, M. Yamamoto, B. D. Simons and S. Yoshida,  _Mouse Spermatogenic Stem Cells Continually Interconvert between Equipotent Singly Isolated and Syncytial States_ ,  Cell Stem Cell 14(5), 658 (2014),  10.1016/j.stem.2014.01.019.
  * [25] Y. Kobayashi, Y. Yasugahira, H. Kitahata, M. Watanabe, K. Natsuga and M. Nagayama,  _Interplay between epidermal stem cell dynamics and dermal deformation_ ,  npj Computational Materials 4(1), 45 (2018),  10.1038/s41524-018-0101-z.
  * [26] M. D. Prater, V. Petit, I. Alasdair Russell, R. R. Giraddi, M. Shehata, S. Menon, R. Schulte, I. Kalajzic, N. Rath, M. F. Olson, D. Metzger, M. M. Faraldo _et al._ ,  _Mammary stem cells have myoepithelial cell properties_ ,  Nature Cell Biology 16(10), 942 (2014),  10.1038/ncb3025.
  * [27] A. Van Keymeulen, A. S. Rocha, M. Ousset, B. Beck, G. Bouvencourt, J. Rock, N. Sharma, S. Dekoninck and C. Blanpain,  _Distinct stem cells contribute to mammary gland development and maintenance_ ,  Nature 479(7372), 189 (2011),  10.1038/nature10573.
  * [28] G. Mascré, S. Dekoninck, B. Drogat, K. K. Youssef, S. Brohée, P. A. Sotiropoulou, B. D. Simons and C. Blanpain,  _Distinct contribution of stem and progenitor cells to epidermal maintenance_ ,  Nature 489(7415), 257 (2012),  10.1038/nature11393.
  * [29] C. Gomez, W. Chua, A. Miremadi, S. Quist, D. Headon and F. Watt,  _The Interfollicular Epidermis of Adult Mouse Tail Comprises Two Distinct Cell Lineages that Are Differentially Regulated by Wnt, Edaradd, and Lrig1_ ,  Stem Cell Reports 1(1), 19 (2013),  10.1016/j.stemcr.2013.04.001.
  * [30] A. D. Lander,  _The 'stem cell' concept: is it holding us back?_ ,  Journal of Biology 8(8), 70 (2009),  10.1186/jbiol177.
  * [31] B. Corominas-Murtra, C. L. G. J. Scheele, K. Kishi, S. I. J. Ellenbroek, B. D. Simons, J. van Rheenen and E. Hannezo,  _Stem cell lineage survival as a noisy competition for niche access_ ,  Proceedings of the National Academy of Sciences 117(29), 16969 (2020),  10.1073/pnas.1921205117.
  * [32] M. Azkanaz, B. Corominas-Murtra, S. I. J. Ellenbroek, L. Bruens, A. T. Webb, D. Laskaris, K. C. Oost, S. J. A. Lafirenze, K. Annusver, H. A. Messal, S. Iqbal, D. J. Flanagan _et al._ ,  _Retrograde movements determine effective stem cell numbers in the intestine_ ,  Nature 607(7919), 548 (2022),  10.1038/s41586-022-04962-0.
  * [33] M. Ito, Y. Liu, Z. Yang, J. Nguyen, F. Liang, R. J. Morris and G. Cotsarelis,  _Stem cells in the hair follicle bulge contribute to wound repair but not to homeostasis of the epidermis_ ,  Nature Medicine 11(12), 1351 (2005),  10.1038/nm1328.
  * [34] V. Jaks, N. Barker, M. Kasper, J. H. van Es, H. J. Snippert, H. Clevers and R. Toftgård,  _Lgr5 marks cycling, yet long-lived, hair follicle stem cells_ ,  Nature Genetics 40(11), 1291 (2008),  10.1038/ng.239.
  * [35] M. E. Page, P. Lombard, F. Ng, B. Göttgens and K. B. Jensen,  _The Epidermis Comprises Autonomous Compartments Maintained by Distinct Stem Cell Populations_ ,  Cell Stem Cell 13(4), 471 (2013),  10.1016/j.stem.2013.07.010.
  * [36] J. Li, S. K. Schnyder, M. S. Turner and R. Yamamoto,  _Role of the Cell Cycle in Collective Cell Dynamics_ ,  Physical Review X 11(3), 031025 (2021),  10.1103/PhysRevX.11.031025.
  * [37] J. Li, S. K. Schnyder, M. S. Turner and R. Yamamoto,  _Competition between cell types under cell cycle regulation with apoptosis_ ,  Physical Review Research 4(3), 033156 (2022),  10.1103/PhysRevResearch.4.033156.
  * [38] P. J. Hoogerbrugge and J. M. V. A. Koelman,  _Simulating Microscopic Hydrodynamic Phenomena with Dissipative Particle Dynamics_ ,  Europhysics Letters (EPL) 19(3), 155 (1992),  10.1209/0295-5075/19/3/001.
  * [39] P. Español and P. Warren,  _Statistical Mechanics of Dissipative Particle Dynamics_ ,  Europhysics Letters (EPL) 30(4), 191 (1995),  10.1209/0295-5075/30/4/001.
  * [40] I. Pagonabarraga, M. H. J. Hagen and D. Frenkel,  _Self-consistent dissipative particle dynamics algorithm_ ,  Europhysics Letters (EPL) 42(4), 377 (1998),  10.1209/epl/i1998-00258-6.
  * [41] J. R. Howse, R. A. L. Jones, A. J. Ryan, T. Gough, R. Vafabakhsh and R. Golestanian,  _Self-Motile Colloidal Particles: From Directed Propulsion to Random Walk_ ,  Physical Review Letters 99(4), 048102 (2007),  10.1103/PhysRevLett.99.048102.
  * [42] J. Elgeti, R. G. Winkler and G. Gompper,  _Physics of microswimmers—single particle motion and collective behavior: a review_ ,  Reports on Progress in Physics 78(5), 056601 (2015),  10.1088/0034-4885/78/5/056601.
  * [43] M. E. Cates and J. Tailleur,  _When are active Brownian particles and run-and-tumble particles equivalent? Consequences for motility-induced phase separation_ ,  EPL (Europhysics Letters) 101(2), 20010 (2013),  10.1209/0295-5075/101/20010.
  * [44] A. P. Solon, M. E. Cates and J. Tailleur,  _Active brownian particles and run-and-tumble particles: A comparative study_ ,  The European Physical Journal Special Topics 224(7), 1231 (2015),  10.1140/epjst/e2015-02457-0.


# An Adaptive Contrastive Learning Model for Spike Sorting

Lang Qian , Cheng Yang
Tsinghua Shenzhen International Graduate School
Tsinghua University
Shenzhen, China
<EMAIL_ADDRESS><EMAIL_ADDRESS>
Shengjie Zheng
Shenzhen Institute of Advanced Technology, Chinese Academy of Sciences
University of Chinese Academy of Sciences
Shenzhen, China
<EMAIL_ADDRESS>
Chunshan Deng, Xiaojian Li
Brain Cognition and Brain Disease Institute (BCBDI)
Shenzhen Institute of Advanced Technology, Chinese Academy of Sciences
Shenzhen, China
<EMAIL_ADDRESS><EMAIL_ADDRESS>
Lang Qian and Shengjie Zheng contributed equally to this workCheng Yang and
Xiaojian Li are co-corresponding authors

###### Abstract

Brain-computer interfaces (BCIs), is ways for electronic devices to
communicate directly with the brain. For most medical-type brain-computer
interface tasks, the activity of multiple units of neurons or local field
potentials is sufficient for decoding. But for BCIs used in neuroscience
research, it is important to separate out the activity of individual neurons.
With the development of large-scale silicon technology and the increasing
number of probe channels, artificially interpreting and labeling spikes is
becoming increasingly impractical. In this paper, we propose a novel modeling
framework: Adaptive Contrastive Learning Model that learns representations
from spikes through contrastive learning based on the maximizing mutual
information loss function as a theoretical basis. Based on the fact that data
with similar features share the same labels whether they are multi-classified
or binary-classified. With this theoretical support, we simplify the multi-
classification problem into multiple binary-classification, improving both the
accuracy and the runtime efficiency. Moreover, we also introduce a series of
enhancements for the spikes, while solving the problem that the classification
effect is affected because of the overlapping spikes.

##  1 Introduction

Brain-computer interfaces (BCIs), are ways for electronic devices to
communicate directly with the brain. The most common brain-computer interface
applications are in medical care, typically in the form of collecting neural
signals from a subject’s brain to help the subject perform tasks such as
spelling wordsWillett et al. (2021), controlling a computer cursorSilversmith
et al. (2021), controlling a mechanical prosthesisEven-Chen et al. (2020), or
paralyzed muscles, and interacting with the outside world. On the other hand,
in the field of neuroscience research, brain-computer interfaces are also
powerful and unique tools for studying brain function. So far, the most
reliable way to record neural activity is to put microelectrodes into the
brain to record local electrical activity, from which we can observe the
action potential. We can assign it to the correct neuron by detecting the
channel where the action potential appears, the amplitude and waveform of
spike, etc. This process of neural information processing is called spike
sortingGibson et al. (2008)Li et al. (2019). However, microelectrodes are
still affected by different factors that make it hard to interpret and
challenge to assign this recorded waveform to the correct neuron, such as
microelectrode floating, spike overlapping, etc.

Nowadays, big data and high performance computing power have allowed deep
learning technology to make rapid development in applications, especially
supervised learning.Although supervised learning has been very successful on
CV and NLP, we cannot directly apply supervised learning to spike sorting.
Even if the model gets high accuracy, the model cannot transfer effectively
because same neurons from different subjects could fire different spike
waveforms. Moreover, the artificial labeling of spike is not only very
troublesome in processing, but also requires professionals to spend a lot of
timeLi et al. (2019)Saif-ur Rehman et al. (2021). At the same time, the use of
supervised learning can cause overfitting of the neural network when spike
overlapping is frequently encountered. Based on this, there is a slow shift
towards unsupervised but non-contrastive learningEom et al. (2021)Seong et al.
(2021) in dealing with spike sorting, which does not use labels, but the
purpose of such models is only to reconstruct the input and extract the
features of the middle layer for clustering, without really extracting the
latent representation of spikes fired by different neurons effectively.

With the rapid development of contrast learning, the ability of contrast
learning to extract the latent representation has caught up with supervised
learning in some aspects. Based on this situation, we propose a new framework
based on contrast learning that (1) simplifies a multi-classification task
into multiple binary-classification tasks, (2) proposes a series of augmented
approaches for spikes, and (3) introduces contrastive learning into spike
sorting. Our approach concerns extracting features of spike data from all
channels at the same time interval, rather than considering a single channel
spike data simultaneously. Because we believe that the multi-channel spike
data carries more valid information and can allow us to achieve better
classification effect of spikes.

In order to better achieve the latent representation of spikes, we use a
method named invariant information clustering which maximizes the mutual
information of raw data and augmented dataJi et al. (2019). In addition, we
used a series of augmentation methods of the data, with the aim of generating
data sharing similar semantics and the same labels as the original data. In
the past two years, contrastive learning has been used both in CVChen et al.
(2020) or in NLPYan et al. (2021)Gao et al. (2021) have had great success on
it. But for spike data augmentation, which is not the same as image or text
augmentation, we propose a series of spike augmentation methods based on the
theoretical basis of electrophysiology. We also tested the accuracy of our
model on a simulated dataset generated by kilosortPachitariu et al. (2016) and
in-vivo data from rat hippocampusMizuseki et al. (2013).

##  2 Method

###  2.1 Data Preprocessing

Step 1: The signal recordings obtained through microelectrodes are not
directly usable due to Local Field Potential (LFP) and high frequency noise.
We need to band-pass filter the raw data from 300Hz to 3000Hz, the data in the
range from 1Hz to 100Hz, which is usually LFP, and the data in the range from
300Hz to 3000Hz, which is the spike we use normally, the data at higher
frequencies are usually high frequency noise Rey et al. (2015), As shown in
Fig.1.

Step 2: When we get the recorded bandpass filtered signal, we need to locate
the spikes, there are many methods for example, Absolute Value, Nonlinear
Energy Operator, Stationary Wavelet Transform Product Gibson et al. (2008),
the task of locating the spikes is not the focus of this study.

Figure 1: The neural information is recorded using an array of electrodes,
each measuring nearby neural activity. The raw neural signal is transformed
into a digital signal by the signal acquisition system. The signals can be
processed by band-pass filters and spike detecting by detector. The image
modified from Zheng et al. (2022).

###  2.2 Data Augmentation

The key element of contrast learning is to construct augmented data with the
same semantic meaning as the original sample. We can use contrast learning to
find the maximum similarity in similar samples and the maximum difference in
different samples and based on that, we can extract the latent representation
in the samples. For example, we can recolor the image or scale the audio
amplitude, but we can still clearly perceive that the augmented data share the
same semantics as the original data Mohsenvand et al. (2020). However, the
augmentation methods of spike data is not exactly the same as the augmentation
of image data, text data, and audio data. To solve the problem of spike data
augmentation, we want to propose a system on multi-channel spike data
augmentation that can guarantee that the augmented sample and the original
sample have the same semantics between them both.

We uphold the idea that problem solving starts from the problem itself. We
augmented the spikes fired by neurons, and we need to try to simulate the
different spike data fired by neurons in the brain due to various possible
influences. This is because the spike fired by the same neuron at different
moments is not exactly the same, but they share similar semantics. After
observing the data fired by neurons and consulting with relevant biological
background experts, we propose five multi-channel spike data augmentation
methods, as shown in Fig. 2.(1) Random noise, (2) DC shifting (3) Horizontal
shifting (left or right), (4) Amplitude scaling, and (5) Spike overlapping.

This paragraph focused on the details of these five data augmentation methods.
the operation is performed on normalized data, with 10 channels and 81
sampling points.(1) Random noise: the distribution of random noise values is a
random normal distribution with a variance of 0.1. (2) DC shifting: obeying a
uniform distribution of - 0.1 to 0.1. (3) Horizontal shifting: Obeying a
uniform distribution of - 20 to 20 and the value is an integer. (4) Amplitude
scaling: Following the uniform distribution of 0.9 to 1.1. (5) Spike
overlapping: randomly selecting a data (including raw data), performing
horizontal offset, then randomly reducing the amplitude and summing it with
the raw data.

Figure 2: Five enhancements of spike data are described, from top to bottom,
as random noise, DC shifting, Horizontal shifting, Amplitude scaling, Spike
overlapping.

###  2.3 Model Architecture

Classifying spike data is not the same as classifying pictures in the
traditional way, because we cannot determine the number of spike
classifications, we cannot treat classifying spike data as a multi-
classification task. In this case, we propose a new spike data classification
model: an adaptive contrast classification model, referred to as ACCM, as
shown in Fig. 3.

Figure 3: Describing the ACCM model architecture

In this model, to give users more choices, we design three encoder
architectures for Backbone Network, as shown in Fig. 3:

(1) A transformer encoder is abbreviated as ACCM-T, as shown in Fig.4(a). This
model not only focuses on the peak and waveform of spike data, but also pays
more attention to the channel-to-channel relationship. Based on this
situation, we chose the Transformer encoder block because of its unique Multi-
Head attention mechanism. Also, to focus more on the channel information, this
model still includes the position encoding operation, which is also suitable
for large batches of data.

(2) A 1D-CNN encoder is abbreviated as ACCM-C, as shown in Fig.4(b). With aim
of extracting the different latent features of the current channel spike, the
data have three 1D-convolution operations. The kernel size of three
1D-convolution operations is all 1. Meanwhile, the output dimensions are
5,10,15 repectively. The reason for this design is to pay more attention to
the characteristics of the spike itself in terms of waveform, amplitude and
other relevant factors. After three 1D-convolution operations, we performed
the concatenate and flatten operations respectively. We also abandoned the use
of L2-Regularization with considering that spike data is very sensitive to
amplitude. Finally, we reduced the dimensionality of the data by performing
several fully connected operations. Meanwhile, we abandon the use of 2D-CNN in
this model, because the relationship between the microelectrode channels is
not really linear, but more like a graph structure. This model encoder is more
suitable for spike data that needs to be sensitive to peaks and patterns.

(3) A recurrent neural network encoder is abbreviated as ACCM-R, as shown in
Fig.4(c). This model is more concerned with channel-to-channel relationships,
when the number of channels of neural electrodes reaches hundreds or
thousands, not all channels can detect spike signals, and the neural
information contained in multiple channels is sufficient to meet the
requirements of spike sorting. This model encoder is more suitable for
channel-sensitive spike data.

At the same time, when we train the model, in order to avoid falling into the
local optimal point, we iteratively trained the two-class and four-class
classifiers, which can help us find the global optimal point. We set up a
controller under each binary-classification branch to decide whether the
classified data should be classified again.

Figure 4: Encoder architectures: (a) Transformer encoder block, (b) 1D-CNN
encoder block, (c) RNN encoder block.

###  2.4 Contrastive loss

In this model, we choosed to maximize the mutual information for contrastive
lossJi et al. (2019). Because we are dealing with a binary-classification
task, so the joint probability distribution was given by the 2x2 matrix
$\mathrm{M}$, as shown in Equation 1.

$\mathrm{M}=\frac{1}{n}\sum_{i=1}^{n}\Phi\left(\mathbf{x}_{i}\right)\cdot\Phi\left(\mathbf{x}_{i}^{\prime}\right)^{\mathrm{T}}$
(1)

The output $\Phi(\mathbf{x})\in[0,1]^{2}$ can be interpreted as the
distribution of a discrete random variable over $2$ classes. As we generally
consider symmetric problems, we usually redefine $\mathrm{M}$ as
$\left(\mathrm{M}+\mathrm{M}^{\mathrm{T}}\right)/2$.

$I(\mathrm{M})=\sum_{c=1}^{2}\sum_{c^{\prime}=1}^{2}\mathrm{M}_{cc^{\prime}}\cdot\ln\frac{\mathrm{M}_{cc^{\prime}}}{\mathrm{M}_{c}\cdot\mathrm{M}_{c}^{\prime}}$
(2)

where maximizing $I(\mathrm{M})$ trades-off minimizing the difference between
the raw data and the augmented data, as shown in Equation 2.

###  2.5 Dataset

In-vivo dataset: We used the hc3 dataset sourced from CRCNSMizuseki et al.
(2013), the data set contains recordings made from multiple hippocampal areas
in Long-Evans rats. The raw (broadband) data was recorded at 20KHz, and we
selected a segment of data with 8 channels, each with 32 sampling points, and
a total of three neurons firing spikes. Total number of spike data is 6000,
each neuron fired a total of 2,000 times.

Simulated data: The simulated dataset we used was derived from kilosort,
having 10 channels, each with 81 sampling points, and a total of eleven
neurons firing spikes. Total number of spike data is 22000, each neuron fired
a total of 2,000 timesPachitariu et al. (2016).

##  3 Experiment

In this section, we compared our results with one unsupervised model AE Eom et
al. (2021) which reconstructed the model own input and two contrastive
learning models, SimCLR Chen et al. (2020), SimCSE Gao et al. (2021). The
initialization weights of the model parameters and the selection of the
optimizer as well as the size of the learning rate are also presented.

###  3.1 Model Training

To train our model, we experimented with three architectures Fig. 4. In each
backbone network, the choice of batch sizes with training data is not exactly
the same. The epochs of each backbone network node are determined by the
number of training data, and usually fluctuate between 15 and 60. After
several model training sessions, we found that initializing the model
parameters to a random normal distribution with variance of 0.1 allowed the
model to converge faster. At the same time, we use the Adam optimizer instead
of the SGD optimizer in order to avoid getting trapped in the local optimal
point, and the learning rate was initialized to 0.005. We performed the
experiments on one GTX 2080ti GPU.

##  4 Results

We tested three different versions of the model, ACCM-R, ACCM-C, and ACCM-T,
and obtained the results of these three versions in terms of accuracy and
running time, respectively, for comparison with an unsupervised model AE Eom
et al. (2021), two contrastive learning models SimCLRChen et al. (2020) and
SimCSEGao et al. (2021).

Table 1: Simulated Model | ACC | Running Time
---|---|---
AE | 83.2% | 6min13s
SimCLR | 87.8% | 8min17s
SimCSE | 61.1% | 3min3s
ACCM-R | 95.3% | 21min3s
ACCM-C | 96.5% | 5min41s
ACCM-T | 98.4% | 11min4s
Table 2: In-vivo Model | ACC | Running Time
---|---|---
AE | 71.2% | 1min24s
SimCLR | 97.9% | 1min46s
SimCSE | 82.1% | 1min21s
ACCM-R | 98.7% | 1min9s
ACCM-C | 99.9% | 1min6s
ACCM-T | 99.9% | 58s

The above table 1 and table 2 shows the results of our experiments, the second
column represents the accuracy rate, and the third column is the model running
time to reach the highest accuracy rate. Our model is provided in three
different versions for different situations. If we need fast running time and
relatively high accuracy, we can choose ACCM-C version. For high efficiency,
we can choose ACCM-T. Also, ACCM-T achieves very high accuracy on simulated
data.

Through the above tables,we observed that the accuracy of model SimCSE is much
less than that of model SimCLR for both in-vivo and simulated datasets. Model
SimCSE augmented spike data only by dropout, while model SimCLR augmented
spike data through our model’s data augmentation methods. This can illustrate
the superiority of our data augmentation methods. Meanwhile, we observe that
model SimCLR is much more accurate when the number of clusters is relatively
small than when the number of clusters is relatively large. This illustrated
that in spike sorting, multiple binary-classification tasks could be easier to
get a higher accuracy rate than the multi-classification task. Again, we
concluded from the above two tables that our model could obtain high accuracy
in both in-vivo and simulated datasets, which is a good indication of the
superiority of our model in handling spike sorting.

##  5 Conclusion

In the field of neuroscience research, brain-computer interfaces provide a way
to study the direct causal relationship between small populations of neurons
and specific external outputs. Identification of the spiking waveform of each
neuron is the basis for distinguishing the activity of individual neurons in a
small population of neurons. In this paper we introduced ACCM, a self-
supervised framework for learning representation for spike data. To be able to
accomplish this work, we proposed five ways of augmenting the spike data. We
proposed three model variants, ACCM-C, ACCM-R, and ACCM-T, because different
professionals have different preferences for accuracy and runtime when
classifying spikes, and for the first time we introduced contrastive learning
into spike sorting, achieving state-of-the-art accuracy without the need for
electrode maps. And, since we cannot determine the number of clusters while
performing spike sorting, we also proposed an adaptive classification
algorithm to transform the multi-classification task into multiple binary-
classification tasks.

## References

  * Chen et al. [2020] Ting Chen, Simon Kornblith, Mohammad Norouzi, and Geoffrey Hinton.  A simple framework for contrastive learning of visual representations.  In _International conference on machine learning_ , pages 1597–1607. PMLR, 2020.
  * Eom et al. [2021] Junsik Eom, In Yong Park, Sewon Kim, Hanbyol Jang, Sanggeon Park, Yeowool Huh, and Dosik Hwang.  Deep-learned spike representations and sorting via an ensemble of auto-encoders.  _Neural Networks_ , 134:131–142, 2021.
  * Even-Chen et al. [2020] Nir Even-Chen, Dante G Muratore, Sergey D Stavisky, Leigh R Hochberg, Jaimie M Henderson, Boris Murmann, and Krishna V Shenoy.  Power-saving design opportunities for wireless intracortical brain–computer interfaces.  _Nature biomedical engineering_ , 4(10):984–996, 2020.
  * Gao et al. [2021] Tianyu Gao, Xingcheng Yao, and Danqi Chen.  Simcse: Simple contrastive learning of sentence embeddings.  _arXiv preprint arXiv:2104.08821_ , 2021.
  * Gibson et al. [2008] Sarah Gibson, Jack W Judy, and Dejan Markovic.  Comparison of spike-sorting algorithms for future hardware implementation.  In _2008 30th Annual International Conference of the IEEE Engineering in Medicine and Biology Society_ , pages 5015–5020. IEEE, 2008.
  * Ji et al. [2019] Xu Ji, Joao F Henriques, and Andrea Vedaldi.  Invariant information clustering for unsupervised image classification and segmentation.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 9865–9874, 2019.
  * Li et al. [2019] Yinghao Li, Shuai Tang, and Virginia R de Sa.  Supervised spike sorting using deep convolutional siamese network and hierarchical clustering.  _unpublished thesis_ , 2019.
  * Mizuseki et al. [2013] Kenji Mizuseki, Anton Sirota, Eva Pastalkova, Kamran Diba, and György Buzsáki.  Multiple single unit recordings from different rat hippocampal and entorhinal regions while the animals were performing multiple behavioral tasks.  _CRCNS org_ , 2013.
  * Mohsenvand et al. [2020] Mostafa Neo Mohsenvand, Mohammad Rasool Izadi, and Pattie Maes.  Contrastive representation learning for electroencephalogram classification.  In _Machine Learning for Health_ , pages 238–253. PMLR, 2020.
  * Pachitariu et al. [2016] Marius Pachitariu, Nicholas A Steinmetz, Shabnam N Kadir, Matteo Carandini, and Kenneth D Harris.  Fast and accurate spike sorting of high-channel count probes with kilosort.  _Advances in neural information processing systems_ , 29, 2016.
  * Rey et al. [2015] Hernan Gonzalo Rey, Carlos Pedreira, and Rodrigo Quian Quiroga.  Past, present and future of spike sorting techniques.  _Brain research bulletin_ , 119:106–117, 2015.
  * Saif-ur Rehman et al. [2021] Muhammad Saif-ur Rehman, Omair Ali, Susanne Dyck, Robin Lienkämper, Marita Metzler, Yaroslav Parpaley, Jörg Wellmer, Charles Liu, Brian Lee, Spencer Kellis, et al.  Spikedeep-classifier: A deep-learning based fully automatic offline spike sorting algorithm.  _Journal of Neural Engineering_ , 18(1):016009, 2021.
  * Seong et al. [2021] Changyu Seong, Wonjae Lee, and Dongsuk Jeon.  A multi-channel spike sorting processor with accurate clustering algorithm using convolutional autoencoder.  _IEEE Transactions on Biomedical Circuits and Systems_ , 15(6):1441–1453, 2021.
  * Silversmith et al. [2021] Daniel B Silversmith, Reza Abiri, Nicholas F Hardy, Nikhilesh Natraj, Adelyn Tu-Chan, Edward F Chang, and Karunesh Ganguly.  Plug-and-play control of a brain–computer interface through neural map stabilization.  _Nature Biotechnology_ , 39(3):326–335, 2021\.
  * Willett et al. [2021] Francis R Willett, Donald T Avansino, Leigh R Hochberg, Jaimie M Henderson, and Krishna V Shenoy.  High-performance brain-to-text communication via handwriting.  _Nature_ , 593(7858):249–254, 2021.
  * Yan et al. [2021] Yuanmeng Yan, Rumei Li, Sirui Wang, Fuzheng Zhang, Wei Wu, and Weiran Xu.  Consert: A contrastive framework for self-supervised sentence representation transfer.  _arXiv preprint arXiv:2105.11741_ , 2021.
  * Zheng et al. [2022] Shengjie Zheng, Wenyi Li, Lang Qian, Chenggang He, and Xiaojian Li.  A spiking neural network based on neural manifold for augmenting intracortical brain-computer interface data.  _arXiv preprint arXiv:2204.05132_ , 2022.