# Emergent properties of collective gene expression patterns in multicellular
systems

Matthew Smart<EMAIL_ADDRESS>Department of Physics, University of
Toronto, Canada  Anton Zilman<EMAIL_ADDRESS>Department of
Physics, University of Toronto, Canada  Institute for Biomedical Engineering,
University of Toronto, Canada

###### Abstract

Multicellular organisms comprise a diverse collection of stable tissues built
from different cell types. It remains unclear how large numbers of interacting
cells can precisely coordinate their gene expression during tissue self-
organization. We develop a generalized model of multicellular gene expression
that includes intracellular and intercellular gene interactions in tissue-like
collectives. We show that tuning the intercellular signaling strength results
in a cascade of transitions from single-cell autonomy towards different self-
organized collective states. Despite an enormous number of possible tissue
states, signaling tends to stabilize a small number of compositionally and
spatially simple tissue types even for disordered interaction networks.
Statistical properties of the stable phenotypes are preserved for different
interaction networks and initial conditions. These results provide a
theoretical framework to study how collections of cells in distinct organisms
robustly self-organize into relatively simple tissues even for complex
interaction networks mediated by large numbers of different molecules, and how
different stable tissues are maintained in homeostasis and disease. The close
alignment between this theoretical model of tissue self-organization and
modern sequencing techniques, particularly spatial transcriptomics, will
enable future applications in broad biological contexts.

††preprint: APS/123-QED

##  I Introduction

Tissue self-organization, enabled by the formation of distinct patterns of
collective gene expression, is a fundamental aspect of multicellular life. It
remains unclear how very large numbers of interacting cells can precisely
coordinate their phenotypes to ensure proper development, homeostasis, and
response to environmental challenges. A key turning point in our understanding
of gene regulation within single cells and tissues was the discovery that
cells can be experimentally “reprogrammed” to an undifferentiated state by
turning specific genes on or off [1]. It is now widely recognized that
cellular phenotypes are highly plastic not only in vitro [1, 2] but also in
vivo. Cell type transitions are critical for proper tissue homeostasis and
wound healing [3, 4, 5, 6], leading to a very dynamic picture of tissue self-
organization. Accordingly, disruption of the self-organized tissue state
through gene regulatory perturbations is implicated in the pathogenesis of
many diseases, most notably cancer [7, 8, 9, 10].

A major challenge for understanding tissue self-organization, in light of this
plasticity, is the complexity of gene regulation. Vast networks of interacting
regulatory molecules underpin the establishment and maintenance of the various
tissue states that compose developed organisms. These genetically encoded
networks are typically quite different from one organism to the next,
emphasizing the need for generalizable models which can provide the synthesis
of this diversity towards general principles. Motivated by the regularity of
cell fate specification during development, Waddington presciently put forward
the concept of an “epigenetic landscape” in the 1950s [11]. This seminal
concept has been mathematically expanded as a result of experimental advances
[12, 13, 14, 15, 16, 17]. In particular, modern experimental techniques such
as single-cell RNA sequencing (scRNA-seq) have generated enormous amounts of
transcriptomic data [18, 19, 20, 21]. Clusters in the data represent stable
regions in single-cell gene expression space, and can therefore be thought of
as dynamical attractors – the valleys of the epigenetic landscape [22]. The
expression signatures of certain clusters tend to align with known sets of
marker genes expressed by in vivo cell types [18, 22, 19]; this has led to the
Human Cell Atlas (HCA) project [19, 20, 21] and related data-driven efforts to
characterize the cell types composition of human tissue from development to
adulthood and in disease. The large number of observed clusters in the scRNA-
seq data indicates that the dynamical system governing single-cell state is
nonlinear and high-dimensional.

It remains an open question how tissues are able to organize into the numerous
and diverse architectures that make up multicellular organisms from the
building blocks of genetically identical cells. Due to significant limitations
in experimental knowledge preceding single-cell techniques, previous
approaches to study tissue self-organization focused on either low-dimensional
theoretical models, which fall into several categories, or on large systems
biology models for specific tissues in specific organisms.

One important category of the low-dimensional models, which includes classical
work [23] and more recent variants [24], explains how spatiotemporal patterns
can arise from reaction-diffusion dynamics of morphogens in a cellular
“continuum”. This type of approach captures certain spatial features of tissue
self-organization, but does not normally describe gene regulatory effects
within single cells, such as cellular autonomy and heterogeneity in the
absence of morphogens. This limitation has been partially addressed by a
different class of approaches that use cellular automata to integrate both
single-cell and tissue level effects [25, 26, 27]. However, due to their
relatively ad-hoc nature it is not clear how to systematically extend them to
larger gene networks associated with scRNA-seq data. To sidestep this
fundamental challenge, there is growing interest in a third class of models,
so-called “gene-free” approaches, which implicitly describe complex gene
regulation within interacting single cells by working with an abstract
phenotype space inspired by the Waddington landscape concept [28, 29, 30, 31,
15, 16]. However, it remains unclear how to relate these abstract variables to
the underlying – and experimentally measurable – gene activity, and how to
incorporate additional cell types or signaling pathways once the model is
constructed.

In contrast to the various types of low-dimensional theoretical frameworks,
systems biology approaches have been widely used to model gene regulation
within cells in different organisms [32], but parameterization difficulties
inherent to many such models make it difficult to draw generalizations from
their applications to specific systems. More comprehensive computational
approaches which additionally consider cell proliferation and other factors
have recently been reviewed [33]. Overall, previous works have commonly
focused on analytically studying cases with a small number of interacting
genes, or use large systems biology models fine-tuned to study specific
experimental systems.

Despite significant advances, it remains unclear how the diverse collections
of tissues that make up developed organisms can be generated and maintained
using genetically pre-determined rules. The key puzzle is the tendency of
multicellular systems to self-organize into stable configurations which are
spatially and compositionally simple, containing only a small fraction of the
possible phenotypic richness suggested by the combinatorics of collective gene
expression. To address this, it is critical to integrate the high-
dimensionality of gene expression at the single-cell level with cell-cell
interactions which can channel this high-dimensionality into relatively low-
complexity tissue states. As a representative example, mammalian cells exhibit
over a hundred stable phenotypic states (cell types), and have on the order of
$10^{4}$ protein-coding genes, with ${\sim}10^{3}$ identified as transcription
factors (genes that regulate the expression of other genes) [34].
Incorporating the high-dimensionality of transcriptional regulation – which
can now be interrogated experimentally using modern sequencing techniques –
provides a path towards addressing questions surrounding the assembly of
diverse tissues as well as cell type plasticity across different tissue
microenvironments.

To tackle this key problem, we use models from statistical physics, which are
naturally suited to describe how interactions between many microscopic degrees
of freedom within and between the cells lead to the formation of stable
macroscopic tissue states. To align our work with the conventional definitions
of cell type introduced above, we focus on minimal models which can encode a
large set of high-dimensional, binarized gene expression patterns as dynamical
attractors. These attractors will morph in the presence of cell-cell
signaling, reminiscent of cellular plasticity in vivo, thereby facilitating
the self-organization of diverse tissue configurations. We specifically employ
a type of spin glass [35, 36] which has been used to describe single-cell
reprogramming experiments [37, 38], generalizing it here to describe
multistability in cellular collectives.

In this paper, we present a model of multicellular gene expression that
couples the transcriptomic states of interacting cells in a systematic,
tunable manner. The model allows us to investigate the interplay between the
intra- and intercellular gene-gene interactions in forming stable collective
tissue states in a framework that is amenable to experimental input and
verification by single-cell transcriptomics data. We demonstrate how different
choices of cell-cell interactions can cause the multicellular system to self-
organize into a broad range of collective spatial patterns. To examine how
cell-cell interactions control the formation of such patterns, we consider an
ensemble of tissues composed of non-interacting cells and tune the strength of
signaling. This reveals a rich sequence of transitions in the space of tissue
gene expression. In the strong signaling regime, we characterize the
distribution of stable tissues and show that it may be partitioned into a
relatively small number of tissue types, offering insight into the assembly of
diverse tissue configurations from the coordination of genetically identical
but phenotypically plastic cells. Our results suggest that several statistical
properties of this distribution are invariant under different realizations of
the random signaling rules, while also shedding light on the puzzling
prevalence of compositionally simple tissues in nature despite the complex
connectivity of the “hairball” of signaling pathways that govern them [39]. We
conclude with a discussion of these results, which have implications for
understanding the self-organization of diverse tissues as well as the
plasticity of the cell types which constitute them. Our analysis may also
inform our understanding of diseases that disrupt the phenotypic composition
of tissues, such as autoimmunity and cancer.

##  II Model

We begin by detailing the spin glass model of multicellular gene expression
illustrated in Fig. 1 which describes the expression state of $M$ cells each
with $N$ genes. First, we introduce the fundamental unit of the model, the
single cell. We then present the multicell model wherein cells interact on a
graph describing their spatial couplings. Finally, we explain how the cell-
cell interactions are incorporated and how the state of each cell is updated.

###  II.1 Single-cell model

We assume that the phenotype of a single cell is defined by its gene
expression pattern, where each gene is in a binary “on” or “off” state. This
simplifying assumption follows Kauffman’s classical work [40]. The state of
the cell is then given by an $N$-dimensional vector
$\mathbf{s}\in\\{+1,-1\\}^{N}$ where $N$ is the number of genes and $s_{i}$
denotes the state of gene $i$.

Stable single-cell phenotypes (cell types) are represented as attractors of
the biological gene regulatory network [41]. The transcriptome of a given cell
type $\mu$ is denoted by $\bm{\xi}^{\mu}$ (as an $N$-dimensional binary
vector), and ongoing sequencing efforts have identified the stable
transcriptomic signatures of a large set of $p$ cell types
$\\{\bm{\xi}^{\mu}\\}_{\mu=1}^{p}$. The stability of each observed cell type
$\bm{\xi}^{\mu}$ is a key constraint for a candidate model of gene expression
dynamics.

Hopfield networks (HN) [35, 42, 36] provide a tractable minimal model to
encode such attractors. Mathematically, an HN is a form of Ising spin glass
and is defined by the Hamiltonian (defining a pseudo-energy landscape)

$H_{0}(\mathbf{s})=-\frac{1}{2}\mathbf{s}^{\mathsf{T}}\mathbf{J}\mathbf{s}-\mathbf{h}^{\mathsf{T}}\mathbf{s}$
(1)

where $\mathbf{h}$ is an applied field on each gene and the gene-gene
interactions $\mathbf{J}$ are chosen so that each pattern $\bm{\xi}^{\mu}$ is
a global minimum of Eq. (1) in the absence of an external field. Choosing
$\mathbf{J}$ based on a rule to encode a set of patterns as minima is commonly
called pattern storage in the HN literature.

The pseudo-energy landscape of the spin model reflects the Waddington
landscape of the gene-gene interaction network, with local minima representing
the stable configurations of gene expression corresponding to the biologically
observed single-cell phenotypes.

Because the gene expression vectors associated with cell types are typically
correlated (i.e. non-orthogonal), we use the projection rule [43, 44] for
pattern storage as in ref. [37]. Given an $N\times p$ matrix of cell type
transcriptomes
$\bm{\xi}=[\bm{\xi}^{1}\>\>\,\bm{\xi}^{2}\>\ldots\>\bm{\xi}^{p}]$, the
projection rule for storing the $p$ patterns is

$\mathbf{J}^{\textrm{proj}}=\bm{\xi}(\bm{\xi}^{\mathsf{T}}\bm{\xi})^{-1}\bm{\xi}^{\mathsf{T}},$
(2)

which specifies an $N\times N$ matrix of gene-gene interactions. As in [43] we
set the diagonal elements $J_{ii}$ to zero.

Minima of the Hamiltonian Eq. (1) correspond to stable single-cell gene
expression states. In the absence of noise or external signals $\mathbf{h}$, a
given cell state $\mathbf{s}\in\\{+1,-1\\}^{N}$ will decrease its energy
$H_{0}(\mathbf{s})$ until it reaches a local minimum $\mathbf{s}^{*}$ of
$H_{0}$. The projection rule Eq. (2) ensures each cell type $\bm{\xi}^{\mu}$
is a global minimum of $H_{0}$ [43, 44].

This framework has been shown to recapitulate aspects of in vitro cellular
reprogramming in single cells [37, 38]. HNs have also been used in a variety
of other biological contexts [45, 46, 47, 48, 49, 50]. However, these works
did not account for the cell-cell interactions within tissue, which play an
essential role in regulating tissue structure and composition. To describe
multicellular systems, we next consider how cells may influence one another’s
gene expression.

###  II.2 Multicellular model

To treat a multicellular system of $M$ cells, we extend Eq. (1) by adding
intercellular interaction terms. The Hamiltonian for the collection of $M$
cells, as a set of $N$-dimensional gene expression vectors
$\\{\mathbf{s}^{a}\\}_{a=1}^{M}$, is

$H(\\{\mathbf{s}^{a}\\})=\sum_{a}H_{0}(\mathbf{s}^{a})+\gamma\sum_{a}\sum_{b}A^{ab}f(\mathbf{s}^{a},\mathbf{s}^{b})$
(3)

where the summation is over individual cells. The first sum describes the $M$
individual cells, while the second sum contains the cell-cell couplings with a
symmetric functional form $f(\mathbf{x},\mathbf{y})=f(\mathbf{y},\mathbf{x})$.
Biologically, these cell-cell interactions can be mediated by a variety of
factors, and we detail a particular form for the coupling inspired by ligand-
receptor signaling in the following subsection. The spatial adjacency matrix
$\mathbf{A}$ defines which cells are interacting, with $A^{ab}=1$ if $a$ and
$b$ are neighbors (interact) and $0$ otherwise (Fig. 1(a)).

The overall interaction strength is quantified by the global parameter
$\gamma\geq 0$. As $\gamma\rightarrow 0$ the tissue acts as a collection of
$M$ independent, non-interacting cells. Beyond a certain threshold $\gamma>0$
the system may exhibit emergent multicellular behavior, such as signaling
dependent cell types or collective spatial patterns.

Figure 1: Multicellular model. (a) Biological tissue consisting of $M$
interacting cells is represented by a graph. Graph edges are defined by the
adjacency matrix $\mathbf{A}$ with $A^{ab}=1$ if cells $a,b$ interact, $0$
otherwise. (b) Single-cell state is defined by the expression status of $N$
genes, $s_{i}\in\\{+1,-1\\},i=1,\ldots,N$. We concisely denote the state of
cell $a$ by the $N$-dimensional vector $\mathbf{s}^{a}$. Every cell has the
same intracellular gene regulatory rules specified by an $N\times N$ matrix
$J_{ij}$. Adjacent cells interact according to a second $N\times N$ matrix
$W_{ij}$ which describes the effect that gene $j$ in a “sender” cell has on
gene $i$ in a “recipient” cell. The cell-cell interaction strength $\gamma$
tunes the magnitude of $\mathbf{W}$ relative to $\mathbf{J}$. (c) The
interaction matrix for tissue gene expression, denoted by
$\bm{\mathcal{J}}(\gamma)$ in Eq. (2), defines interactions between all $NM$
genes in the multicellular system. To illustrate its construction, we show
$\bm{\mathcal{J}}(\gamma=0.5)$ for a particular choice of
$\mathbf{J}$,$\mathbf{W}$, and $\mathbf{A}$.

For simplicity, we use an adjacency matrix corresponding to next-nearest-
neighbor interactions on a square lattice. This choice is biologically
inspired by cells interacting over relatively short distances (e.g. Notch
signaling [31]). While we focus on the two-dimensional case for ease of
visualization, an arbitrary cell-cell interaction matrix can be chosen in
principle. Likewise, the interaction strength $\gamma$ could also be distance-
dependent instead of constant.

###  II.3 Cell-cell interactions

Many genes participate in cell-cell interactions by generating, sensing, or
transducing signals that are sent between cells. For example, a cell may
express and secrete signaling molecules (ligands) into the surrounding
environment, which in turn affects the gene expression of neighboring cells
that sense and respond to these signals. Each ligand may influence the
expression of many genes in a target cell resulting in an intercellular gene-
gene interaction network. We represent this network of sender-recipient
signaling interactions via an $N\times N$ matrix $\mathbf{W}$, where $W_{ij}$
represents the effect gene $j$ (in a sender cell) has on gene $i$ (in a
recipient cell). A complete empirical $\mathbf{W}$ is not yet available from
the biological data and so we will focus on the case of randomly sampled
$\mathbf{W}$.

The cell-cell signaling matrix $\mathbf{W}$ effectively couples the gene
regulatory networks of neighboring cells, depicted graphically in Fig. 1(b).
We denote the “signaling field” which cell $b$ exerts on a neighboring cell
$a$ by $\mathbf{h}^{ab}=\gamma\mathbf{W}\mathbf{s}^{b}$, which acts as an
applied field on the single-cell Hamiltonian for cell $a$,
$H_{0}(\mathbf{s}^{a})$ (Eq. (1)). Summing over all neighbors gives the total
applied field that the tissue collectively exerts on cell $a$,

$\mathbf{h}^{a}=\gamma\mathbf{W}\sum_{b}A^{ab}\mathbf{s}^{b}.$ (4)

Note that this collective applied field is not static; it changes with the
state of the tissue. This aspect is essential for the self-organizing
properties of the model.

The collective applied field Eq. (4) we have considered for our initial
investigation depends only on the state of the sender cells. In principle,
this function could also depend on the state of the recipient cell. For
example, a cell that is not expressing certain receptors will not be able to
sense and respond to the corresponding ligands. This more complex signaling
form will be explored in future work.

When both $\mathbf{J}$ and $\mathbf{W}$ are symmetric, the collective behavior
can be studied from a statistical mechanics perspective through the
multicellular Hamiltonian introduced above. The choice of collective applied
field in Eq. (4) corresponds to choosing
$f(\mathbf{s}^{a},\mathbf{s}^{b})=-\frac{1}{2}{\mathbf{s}^{a}}^{\mathsf{T}}\mathbf{W}\mathbf{s}^{b}$
for the interaction terms in Eq. (3). The parameterized multicellular
Hamiltonian is

$H(\\{\mathbf{s}^{a}\\})=-\frac{1}{2}\sum_{a}{\mathbf{s}^{a}}^{\mathsf{T}}\mathbf{J}\mathbf{s}^{a}-\frac{\gamma}{2}\sum_{a}\sum_{b}A^{ab}{\mathbf{s}^{a}}^{\mathsf{T}}\mathbf{W}\mathbf{s}^{b}.$
(5)

Equation (5) is compactly expressed as
$H(\mathbf{x},\gamma)=-\frac{1}{2}\mathbf{x}^{\mathsf{T}}\bm{\mathcal{J}}(\gamma)\mathbf{x}$,
where $x$ is an $NM$ vector of each cell’s transcriptome concatenated, and the
$NM\times NM$ interaction matrix $\bm{\mathcal{J}}(\gamma)$ has a block form,

$\bm{\mathcal{J}}(\gamma)=\mathbf{I}_{M}\otimes\mathbf{J}+\mathbf{A}\otimes\gamma\mathbf{W},$
(6)

with $\otimes$ denoting the Kronecker product for constructing block matrices
and $\mathbf{I}_{M}$ the $M\times M$ identity matrix. Each term represents a
distinct layer of gene regulation: the diagonal blocks correspond to the
intracellular gene interactions $\mathbf{J}$, whereas the off-diagonal blocks
correspond to the intercellular signaling matrix $\gamma\mathbf{W}$ tiled
according to the adjacency matrix $\mathbf{A}$ (see Fig. 1(c)). The cell-cell
interaction strength $\gamma$ tunes the magnitude of $\mathbf{W}$ relative to
$\mathbf{J}$ (Fig. 1(b)).

Various types of block Ising models have been studied in different contexts
[51, 52, 53, 54], but they have largely been restricted to either two cells or
to uniform off-diagonal interactions (i.e. constant $W_{ij}$ or simple
adjacency $A^{ab}=1-\delta^{ab}$). We are focused here on much more general
off-diagonal interactions (randomly sampled $W_{ij}$ and structured adjacency
matrices; see Results). Of note, refs. [55, 56, 57] studied disordered “multi-
species” block systems with an arbitrary number of cells, but they consider
alternative couplings between cells and do not focus on the deterministic
limit.

Although in this work we consider discrete gene expression states, we note
that related continuous state models known as coupled map lattices [58] have
been used to describe lattices of interacting cells [59, 60, 61]. These works
focused on a few genes or underlying cell types, whereas the approach we
outline is inherently scalable to many genes and cell types.

###  II.4 Gene expression dynamics

The minima of Eq. (5) correspond to stable configurations of the tissue
(collective gene expression patterns). To identify them, we use a discrete
analog of gradient descent on Eq. (5). Our results are also relevant in the
case of mild gene expression noise. We present the full stochastic update rule
here for completeness, then introduce the deterministic limit.

We use Glauber dynamics [35, 62] as an asynchronous update rule for the
single-cell spin glass Eq. (1) and its multicell extension Eq. (5). We
emphasize that we are not focused on the dynamics itself, but rather in using
it as a tool to sample the steady states which arise for interacting cells.
For a given cell, a gene is selected at random and updated according to

$p(s_{i}(t+1)\rightarrow 1)=\frac{1}{1+\exp{(-2\beta
h_{i}^{\textrm{total}})}}$ (7)

where $h_{i}^{\textrm{total}}=\sum_{j}J_{ij}s_{j}(t)+h_{i}$, and $\beta^{-1}$
represents the strength of the gene regulatory noise arising from various
sources (analogous to thermal noise). The timestep $t$ is expressed in units
of single gene updates. Note that the mean spin update is $\langle
s_{i}(t+1)\rangle=\tanh{(\beta h_{i}^{\textrm{total}})}$.

In the deterministic limit ($\beta\rightarrow\infty$) of Eq. (7), the update
rule for a single cell becomes
$\mathbf{s}(t+N)=\textrm{sgn}(\mathbf{J}\mathbf{s}(t)+\mathbf{h})$, where
$\textrm{sgn}(\cdot)$ is applied element-wise and the $N$ genes are updated in
a fixed sequential order. We fix the sequence of updates in order to ensure a
well-defined mapping from an initial condition to a resulting fixed point. A
state is a fixed point of the update rule when
$\mathbf{s}=\textrm{sgn}(\mathbf{Js+h})$. The sequential update order does not
impact whether a state is a fixed point. For $\mathbf{h}=\mathbf{0}$, it can
be verified that the encoded cell types $\bm{\xi}^{\mu}$ are fixed points when
$\mathbf{J}$ is defined via the projection rule Eq. (2).

For the multicellular model, the deterministic update rule for each cell is

$\mathbf{s}^{a}(t+N)=\textrm{sgn}(\mathbf{J}\mathbf{s}^{a}(t)+\gamma\mathbf{W}\sum_{b}A^{ab}\mathbf{s}^{b}(t))$
(8)

where $a\in\\{1,\ldots,M\\}$. Tissue level updates can be expressed compactly
using Eq. (6) as
$\mathbf{x}(t+NM)=\textrm{sgn}(\bm{\mathcal{J}}\mathbf{x}(t))$.

Very similar systems of equations have been used as continuous-time dynamical
systems, most notably in Hopfield’s classical work on associative memory [63].
They take the form
$\frac{\mathrm{d}u_{i}}{\mathrm{d}t}=-\frac{u_{i}}{\tau_{i}}+\sum_{j}\tanh{(\beta
v_{j})}$ with $v_{i}\equiv\sum_{j}J_{ij}u_{j}+h_{i}$, $\tau_{i}>0$. An
analogous system has been applied to scRNA-seq data, where it generated
experimentally validated predictions in the context of differentiation [22].
Interestingly, such systems have also been used as recurrent neural networks
which may be trained to reproduce time series from other dynamical systems
[64].

When either of $\mathbf{J}$, $\mathbf{W}$ are asymmetric, the stochastic
dynamics no longer satisfy detailed balance. In this case, the model is known
as an asymmetric kinetic Ising system [65, 66, 67]. In addition to fixed point
attractors, such systems can exhibit oscillatory behavior which is necessary
to describe phenomena such as the cell cycle or spatiotemporal patterns. This
out-of-equilibrium dynamics will be investigated in future work.

###  II.5 Low-dimensional system with three cell types and nine genes

For simplicity and to facilitate visualization, we consider a low-dimensional
system with $N=9$ genes and $p=3$ encoded single-cell types (shown in Fig.
2(a)). Three cell types is the minimal non-trivial encoding, and nine genes
gives a large but tractable space of $512$ transcriptomic states for each
cell. The set of cell type gene expression vectors $\\{\bm{\xi}^{\mu}\\}$
determines the intracellular gene regulatory interactions $\mathbf{J}$ through
the projection rule Eq. (2).

The single-cell energy landscape Eq. (1) is depicted in Fig. 2(b) for the
single-cell types from Fig. 2(a). In addition to the three encoded minima,
there are five “spurious” minima consisting of the negation of each cell type,
$-\bm{\xi}^{\mu}$ (due to spin-flip symmetry of $H(\mathbf{s})$), and the sum
of the three types,
$\pm\mathbf{S}=\pm\frac{1}{3}(\bm{\xi}^{1}+\bm{\xi}^{2}+\bm{\xi}^{3})$.
Although some spurious minima may be biologically significant [37], others are
likely not (e.g. the negations $-\bm{\xi}^{\mu}$). We assume throughout this
paper that $\mathbf{h}=\mathbf{0}$, but we note that certain minima can be
eliminated using a global external field. Briefly, the embedding in Fig. 2(b)
is generated by computing the Hamming distance between four minima
$\\{\bm{\xi}^{1}$, $\bm{\xi}^{2}$, $\bm{\xi}^{3}$, $\mathbf{S}\\}$ and each of
the $2^{N}$ possible gene expression states. This results in a matrix
$\mathbf{X}\in\mathbb{R}^{4\times 2^{N}}$. Principle component analysis (PCA)
is then applied to $\mathbf{X}$ to obtain a two-dimensional embedding
$\tilde{\mathbf{X}}\in\mathbb{R}^{2\times 2^{N}}$. The columns of
$\tilde{\mathbf{X}}$ correspond to embeddings of each gene expression state
$\mathbf{s}$, which we plot in three dimensions using $H_{0}(\mathbf{s})$ as
the vertical axis.

In the absence of signaling ($\gamma=0$), the stable multicellular states are
composed of prescribed single-cell types $\\{\bm{\xi}^{\mu}\\}$. To
characterize the stable tissues which may self-organize in a given
multicellular system for $\gamma>0$, we sample local minima of Eq. (5) from
arbitrary initial conditions for a given choice of intracellular interactions
$\mathbf{J}$, cell-cell interactions $\mathbf{W}$, and adjacency matrix
$\mathbf{A}$.

As mentioned previously, the cell-cell interaction matrix $\mathbf{W}$ is not
yet fully characterized experimentally. In addition to the many genes that
have direct signaling roles (e.g. ligands and receptors), even more genes
participate indirectly in signaling networks (e.g. as co-factors, downstream
sensing molecules, activatable transcription factors, etc.) or as cargo for
extracellular vesicles such as exosomes which can shuttle RNA between cells
[68]. In lieu of definitive data constraining $\mathbf{W}$ (in contrast to the
single-cell transcriptomics data constraining $\mathbf{J}$), we consider dense
symmetric matrices, sampling the upper triangular elements as $W_{ij}\sim
U[-1,1]$ (see e.g. Fig. 1(c)). We consider structured alternatives to this
simple choice in the discussion.

##  III Results

###  III.1 Different intercellular signaling rules stabilize qualitatively
distinct tissue types

We are first interested in assessing the range of possible tissue states which
self-organize under different signaling rules. In particular, we fix all
aspects of the model except for the signaling matrix $\mathbf{W}$ (i.e.
$\mathbf{J}$, $\mathbf{A}$, and $\gamma$). As described above, $\mathbf{J}$ is
set by the choice of encoded single-cell types $\\{\bm{\xi}^{\mu}\\}$ which
are stable in the absence of signaling (Figs. 2(a) and 2(b)), and $\mathbf{A}$
represents a next-nearest-neighbor square lattice. We then sample different
realizations of $W_{ij}\sim U[-1,1]$ and identify the tissues which self-
organize starting from a fixed initial condition of the tissue gene
expression. See the preceding Section II.5 for details of the simulated
system.

Figure 2: Stable tissue patterns for varying signaling rules. (a) An example
system with $N=9$ genes and $p=3$ single-cell types
$\\{\bm{\xi}^{\mu}\\}_{\mu=1}^{3}$ defined by their gene expression patterns
(arrow up/down denotes on/off). Each cell state is represented by a grid of up
to $N$ white dots where the presence (absence) of a dot indicates the
associated gene is on (off). Additionally, each unique cell state
$\mathbf{s}\in 2^{N}$ is assigned a unique color. The choice of patterns fixes
the intracellular rules $\mathbf{J}$ through the projection rule Eq. (2). (b)
Two-dimensional embedding of the single-cell energy landscape
$H_{0}(\mathbf{s})$ Eq. (1). (c) A $400$ cell ($20\times 20$ grid) lattice
initial condition is prepared with $\bm{\xi}^{2}$ cells on the left half and
$\bm{\xi}^{1}$ cells on the right half. (d) Stable tissue patterns reached
from the initial condition in panel (c) are found for eight random symmetric
signaling rules $W_{ij}\sim U[-1,1]$ with $\gamma=1$.

In Fig. 2(c) we prepare an initial condition of $M=400$ cells arranged on a
$20\times 20$ lattice. This choice mimics a 2D sheet consisting of two
different cell types: the left half is composed of cells in state
$\bm{\xi}^{2}$ and the right half is composed of cells in state
$\bm{\xi}^{1}$. This tissue state evolves according to the regulatory rules,
eventually reaching a local minimum of Eq. (5). Each square in the grid
describes a cell at a particular location. The gene expression state of each
cell is visualized both quantitatively (through the presence/absence of $N=9$
dots corresponding to “on” genes inside the square as in Fig. 2(a)) and
qualitatively (each state $\mathbf{s}\in 2^{N}$ is assigned a unique color).

We identify fixed points reached from the initial condition in Fig. 2(c) for
different random interaction matrices W when the signaling strength $\gamma$
is set to $1$. Representative examples are displayed in Fig. 2(d).
Sufficiently strong interactions can destabilize the encoded single-cell types
$\\{\bm{\xi}^{\mu}\\}$. The stable configurations reached by different choices
of W are quite diverse and can be divided into several “tissue types”. These
include homogeneous (all cells are in the same state), ordered layers, and
labyrinthine (a few states which are spatially interwoven). More heterogeneous
patterns are possible, as are multi-phasic patterns where different spatial
regions of the tissue exhibit different patterns (such as homogeneous in one
region and labyrinthine in another). The spatial features of the observed
patterns also scale with the signaling range (Fig. 3).

Figure 3: Tissue patterns scale with signal radius. We fix the cell-cell
signaling rules $\mathbf{W}$ which stabilized the upper-left tissue state in
Fig. 2(d). Details of the visualization approach are provided in Fig. 2. Each
tissue is generated from the same dual initial condition as in Fig. 2(c) for
$\gamma=1$. The interaction radius $R$, which defines a corresponding
adjacency matrix $\mathbf{A}_{R}$, is indicated above each tissue. The white
shaded region denotes the neighborhood which exchanges signals with the
indicated cell.

Overall, different signaling rules $\mathbf{W}$ can cause an arbitrary
configuration of single-cells to self-organize into a variety of stable tissue
configurations. These patterns are maintained through coherent signaling
between the cells within their local neighborhoods (Eq. (4)). This tissue
self-organization process is important for biological functions such as
homeostasis in different organisms but is difficult to characterize
experimentally due to the high-dimensionality of gene expression. This
theoretical approach describes the formation and stabilization of collective
gene expression patterns for particular choices of the gene regulatory
parameters ($\mathbf{J}$,$\mathbf{W}$,$\gamma$) and spatial organization of
cells ($\mathbf{A}$), which may aid understanding of how particular stable
states arise both in vivo and in vitro.

###  III.2 Tuning the signaling strength destabilizes single-cell types,
leading to a rich sequence of tissue transitions

Complex organisms exhibit a broad range of stable tissue configurations. It
remains unclear how such diversity arises in the context of genetically
predefined intracellular and intercellular regulatory rules (i.e.
$\mathbf{J}$, $\mathbf{W}$ are fixed by the genome). In a given organism there
are numerous tissue types that are distinguished by their cellular composition
and spatial patterning. To carry out different functional roles, such tissue
types are often strikingly distinct (e.g. adipose tissue and hepatic lobules).
In addition to the macroscopic variation across tissue types, there can also
be fine-scale variation in a given type. In this case, the tissues are
constructed from the same set of cell states, but there is minor variation in
cell number or spatial arrangement (e.g. the fingerprints of identical twins).

Here we characterize the diversity of stable tissue configurations in the
multicellular model when the regulatory parameters $\mathbf{J}$, $\mathbf{W}$,
and $\mathbf{A}$ are fixed. We analyze how this diversity emerges as the cell-
cell signaling strength $\gamma$ is tuned, as might occur during development
or as part of homeostasis and adaptation to environmental pressures.

As an illustrative example, in Fig. 4(a) we reconsider the initial condition
from Fig. 2(c) and gradually increase $\gamma$ for a particular choice of
cell-cell signaling rules $\mathbf{W}$. We keep the same $\mathbf{J}$,
$\mathbf{A}$ as in Fig. 2 but lower the size of the system to a $10\times 10$
lattice ($M=100$ cells) to facilitate computation and visualization. The
initial configuration (Fig. 4(a), left) is composed of single-cell types that
are stable in isolation (i.e. in the absence of signaling), and it therefore
remains stable for very mild levels of signaling. However, once $\gamma$
passes a certain threshold, the tissue is destabilized and self-organizes into
a different stable configuration. These multicellular patterns are
characterized by signaling-dependent single-cell states (i.e. they are
maintained by cell-cell interactions). As $\gamma$ increases the tissue
undergoes a cascade of such transitions which are punctuated by intervals of
stability. At $\gamma\gtrsim 1$ the tissue reaches a limiting configuration
that no longer changes, which we term the strong signaling regime.

Fig. 4(a) displays only a few of the many tissue transitions which are
observed from this particular initial condition. While the transition sequence
appears quite complex overall, the earliest transitions at low $\gamma$ may be
anticipated. By definition, a tissue state is stable if all cells
$\\{\mathbf{s}^{a}\\}$ present in the tissue are stable, as defined by Eq.
(8). Because the tissue states at $\gamma=0$ consist of only stable single-
cell types $\bm{\xi}^{\mu}$, we can readily enumerate the possible
neighborhoods (which are defined by the adjacency matrix $\mathbf{A}$). These
neighborhoods determine the collective signaling field experienced by a given
cell, Eq. (4). The simplest neighborhood present in the initial condition in
Fig. 4(a) consists of a cell in state $\bm{\xi}^{\mu}$ surrounded by $z\equiv
8$ cells also in state $\bm{\xi}^{\mu}$ ($z$ is the coordination number for
$\mathbf{A}$). We therefore ask: at what $\gamma$ does the fixed point
condition $\bm{\xi}^{\mu}=\textrm{sgn}(\bm{\xi}^{\mu}+\gamma
z\mathbf{W}\bm{\xi}^{\mu})$ no longer hold?

In this example, the right half of the lattice is composed of $\bm{\xi}^{1}$
while the left half is composed of $\bm{\xi}^{2}$. For $\mu=1$ this criterion
gives $\gamma_{1}^{*}\approx 0.0231$, whereas for $\mu=2$ it gives
$\gamma_{2}^{*}\approx 0.0559$. This analysis is reflected in Fig. 4(a), which
shows that by $\gamma=0.024$ the right half is destabilized but not the left,
and by $\gamma=0.056$ the left half becomes destabilized as well. A general
consequence of this analysis is that the encoded single-cell types
$\\{\bm{\xi}^{\mu}\\}_{\mu=1}^{p}$ are destabilized in a sequential fashion,
from which we identify at least $\sim p$ low-$\gamma$ transitions. This
heuristic is less useful for describing the many transitions at intermediate
levels of $\gamma$ because it requires enumerating all possible neighborhoods
and the combinatorics become prohibitive.

###  III.3 Nonlinear dimension reduction reveals the emergence of self-
organized tissue types

To generalize beyond a single initial condition, we numerically investigate
the distribution of stable gene expression states (of the entire tissue)
reached by an ensemble of initial conditions. Specifically, we choose
$k=10^{4}$ random initial conditions $\\{\mathbf{x}_{i}^{0}\\}_{i=1}^{k}$ and
for each we compute the stable tissue $\mathbf{x}_{i}$ which self-organizes at
a particular value of $\gamma$. This mapping generates a large matrix
$\mathbf{X}_{\gamma}\in\\{+1,-1\\}^{NM\times k}$. Each column of
$\mathbf{X}_{\gamma}$, denoted $\mathbf{x}_{i}$, is a local minimum of Eq.
(5).

Classically, low-dimensional “order parameters” (e.g., the magnetization
$m=N^{-1}\sum_{i}s_{i}$ in the case of a ferromagnet) have served as the
quantification of the degree of disorder and a way to categorize different
macroscopic phases. However, the rich space of distinct minima in our systems
makes direct identification of order parameters impractical. Instead we use an
unsupervised nonlinear dimension reduction technique, Uniform Manifold
Approximation and Projection (UMAP) [69], to embed and identify the tissue
states that belong to the same class based on their phenotypic composition and
spatial gene expression patterns. Unsupervised learning techniques are
increasingly being leveraged to study complex physical systems [70, 71, 72],
and a similar approach to the one taken here has recently been used to
visualize high-dimensional energy landscapes in materials science [73].

Despite the high-dimensionality of the data, UMAP is able to generate
informative two-dimensional embeddings. In Fig. 4 (b) and (c) we display the
nonlinear embedding of $\mathbf{X}_{\gamma}$ for several representative values
of $\gamma$. Each point in a given panel corresponds to a stable tissue
$\mathbf{x}_{i}$. In Fig. 4(b) we color these points by $n(\mathbf{x}_{i})$,
which denotes the number of unique single-cell states that are present in the
tissue, whereas Fig. 4(c) shows the same points colored according to their
relative energy, $H(\mathbf{x}_{i})$.

Figure 4: Emergence of clusters of stable tissue patterns through cell-cell
signaling. We fix $\mathbf{J}$ as in Fig. 2 and sample a particular set of
cell-cell signaling rules $W_{ij}\sim U[-1,1]$. We use the same type of
adjacency matrix $\mathbf{A}$ as in Fig. 2 but with a $10\times 10$ lattice
($M=100$ cells). (a) From a specific tissue initial condition, we increase
$\gamma$ and identify the resulting stable configuration (the same initial
condition is used for each $\gamma$). Numerous tissue transitions are observed
as $\gamma$ is tuned (see Fig. 2 for visualization details). (b) We fix
$k=10^{4}$ random initial conditions and compute the corresponding stable
tissue states $\mathbf{x}_{i}\in\\{+1,-1\\}^{NM}$ ($N$ genes, $M$ cells) for
eight different values of $\gamma$. For each $\gamma$, this generates a matrix
$\mathbf{X}\in\\{+1,-1\\}^{NM\times k}$. We embed $\mathbf{X}$ in two
dimensions using nonlinear dimension reduction (aligned UMAP [69]). Each point
$\mathbf{x}_{i}$ is colored according to the number of unique single-cell
states present in the tissue, $n(\mathbf{x}_{i})$. (c) Each point in (b) is
colored instead by its multicellular energy $H(\mathbf{x}_{i})$.

At $\gamma=0$, the points are embedded in a relatively unstructured manner.
This is expected and suggests that the way we are sampling local minima is not
overtly biased. A heuristic argument for this observation is as follows: in
the absence of interactions ($\gamma=0$), a tissue state is stable whenever
each cell is stable. The number of stable single-cell states $n$ includes the
$p$ encoded single-cell types $\\{\bm{\xi}^{\mu}\\}$ but can also include
additional spurious stable states. For the system with $p=3$ encoded cell
types considered here, there are n=8 stable states (Fig. 2(b)). This gives
$n^{M}$ distinct stable tissue states. When this very large space of $n^{M}$
minima is sub-sampled at random, one expects little structure in the nearest-
neighbor graph. This disorder is ultimately reflected in the embedding.

As $\gamma$ increases this unstructured fine-scale diversity of non-
interacting cells in arbitrary spatial arrangements gradually coalesces into a
relatively small number of clusters. Interestingly, these clusters have a
clear signature in terms of $n(\mathbf{x})$ – the unique single-cell states
that are present in the tissue – as well as the energy $H(\mathbf{x})$ Eq. (5)
(as shown in Fig. 4 (b) and (c)). This indicates that the embedding is
clustering distinct “types” of tissue (in analogy to cell types) which
maintain some intra-cluster variation. For instance, for $\gamma\geq 0.15$
there are two large clusters which contain tissues composed of just two cell
states.

Overall, tuning the interaction strength $\gamma$ promotes the self-
organization of a wide array of stable tissue configurations. At $\gamma=0$,
the model is simply a collection of $M$ non-interacting single-cell types. As
$\gamma$ increases, multicellular structures emerge through collective
interactions between the cells. We have shown that the low-$\gamma$
transitions are associated with sequential destabilization of the encoded
single-cell types, whereas the transitions for higher $\gamma$ are more
nuanced. In the next subsection, we show that the unsupervised approach we
take here is especially useful for identifying and distinguishing tissue types
in the strong signaling regime.

###  III.4 Strong signaling causes the tissue to self-organize into a
relatively small number of types

In the preceding subsection, we applied nonlinear dimension reduction to the
stable gene expression patterns $\mathbf{X}_{\gamma}$ reached from an ensemble
of random initial conditions. We showed that the self-organized tissue states
form several clusters in the low-dimensional space as $\gamma$ is increased.
Here we focus on $\mathbf{X}_{\gamma}$ in the limiting regime of strong
signaling ($\gamma=1$). Fig. 5(a) shows the UMAP embedding of
$\mathbf{X}_{\gamma=1}$ colored by the number of unique single-cell states
present in the tissue, $n(\mathbf{x}_{i})$. Several representative points from
each cluster are annotated and visualized in Fig. 5(b).

Figure 5: Distribution of stable tissue configurations for strong signaling.
(a) The $\gamma=1$ sub-panel of Fig. 4(b) is shown with (b) twelve
representative examples from the observed clusters (see Fig. 2 for
visualization details). (c) For each point $\mathbf{x}_{i}$ from panel (a),
the number of unique single-cell states present in the tissue,
$n(\mathbf{x}_{i})$, is plotted against the energy $E=H(\mathbf{x}_{i})$ (Eq.
(5)). Histograms are shown for both distributions. (d) The embedding from
panel (a) is colored according to the multicell Hamiltonian
$H(\mathbf{x}_{i})$.

By comparing the annotated points within and between the clusters, we conclude
that this unsupervised approach is producing intuitive clusters which group
very similar tissue states (e.g. Fig. 5(b), examples $4,5,6$) while separating
very different ones. In analogy to “cell types”, which are defined based on
clusters of scRNA-seq data (and specified by gene expression signatures
$\bm{\xi}^{\mu}$), we refer to these clusters as “tissue types”. In
particular, we observe clusters of labyrinthine and stripe-like patterns
wherein each cell is in one of two specific cell states (Fig. 5(b), examples
$4-6,11,12$), clusters where the tissue gene expression is homogeneous (Fig.
5(b), examples $1$ and $9$), and clusters where two of the aforementioned
tissue types are separated by an interface (Fig. 5(b), examples $2,7,8,10$).

We also note that within each cluster described above, there can be extensive
fine-scale variation (for instance, in the spatial arrangement of cell states
within the labyrinthine clusters). This suggests a hierarchical picture of the
multicellular energy landscape Eq. (5). At a high level, the landscape is
partitioned into several basins of attraction corresponding to the different
tissue types. In more detail, each of these basins may be locally very rugged
(i.e. contain many local minima in close proximity), reflecting the large
number of similar but distinct stable configurations we observe in certain
clusters.

Additionally, these clusters appear in symmetric pairs. Inspection of the
elements of each pair reveals that they have opposite gene expression
patterns: $\\{\mathbf{s}^{a}\\}\rightarrow\\{-\mathbf{s}^{a}\\}$, compare e.g.
Fig. 5(b) points $6$ and $12$. This is a reflection of the spin-flip symmetry
present in Eq. (5), and indicates that the way the local minima are being
sampled (i.e. the ensemble size and dynamical update rule) is capturing
expected aspects of the energy landscape in a relatively unbiased manner. As
an aside, the antisymmetric minima can be eliminated by applying an external
field to the gene expression of each cell (e.g. by biasing certain
“housekeeping” genes to remain on), and this will be investigated in further
work.

Separately, we also report the distribution of both $n(\mathbf{x}_{i})$ and
the multicellular energy $E=H(\mathbf{x}_{i})$ over all sampled minima
$\\{\mathbf{x}_{i}\\}_{i=1}^{k}$ of $H(\mathbf{x})$. We display the data in
Fig. 5(c), which contains a scatter plot as well as the two marginal
distributions (i.e. $p(n)$ and $p(E)$). In Fig. 5(d) we provide a version of
Fig. 5(a) colored instead by $E=H(\mathbf{x}_{i})$.

The data exhibits three main features. First, there is a clear correlation
between $n(\mathbf{x}_{i})$ and $H(\mathbf{x}_{i})$ (note this is also
apparent when comparing Fig. 5 (a) and (d)). This means that the deepest
minima tend to be compositionally simple, that is, characterized by a low
number of unique single-cell states (low $n$ – e.g. point $1$ in Fig. 5(b)).
In contrast, the higher energy minima tend to be more complex, with large
numbers of distinct single-cell states $n$, illustrated by point $3$ in Fig.
5(b). Second, these low energy, simple minima are also the most probable (over
an ensemble of random initial conditions). This is reflected in the small $n$,
low energy peaks in $p(n)$ and $p(E)$. In terms of the energy landscape, this
suggests that the associated basins of attraction have large volume. And
third, looking specifically at the distribution of the number of single-cell
states within a tissue, $p(n)$, there is a smooth “bulk” of compositionally
complex minima with $5<n<35$ which appears bimodal.

In the following subsection we investigate how these properties of the local
minima of Eq. (5) depend on the choice of random cell-cell signaling rules
$\mathbf{W}$.

###  III.5 Distribution of stable tissues under different random signaling
rules displays universal characteristics

As above we denote the frequency of local minima with $n$ unique single-cell
states by $p(n)$, and the frequency of local minima with energy
$E=H(\mathbf{x})$ by $p(E)$. We select eight different signaling rules
$\\{\mathbf{W}_{\alpha}\\}_{\alpha=1}^{8}$ with symmetric elements
$W_{ij}=W_{ji}\sim U[-1,1]$. For each, we sample the local minima (self-
organized tissue configurations) reached by an ensemble of random initial
conditions. Each panel of Fig. 6 corresponds to a particular
$\mathbf{W}_{\alpha}$, and shows a scatter plot of $n(\mathbf{x}_{i})$,
$H(\mathbf{x}_{i})$ for the sampled minima $\\{\mathbf{x}_{i}\\}$ as well as
the marginal distributions $p(n)$, $p(E)$.

Figure 6: Distribution of stable tissues under different signaling rules.
Eight signaling matrices $\\{\mathbf{W}_{\alpha}\\}$ are sampled according to
$W_{ij}\sim U[-1,1]$. For each $\mathbf{W}_{\alpha}$, the space of stable
tissues (local minima of Eq. (5)) is sampled via $k=10^{4}$ random initial
conditions at $\gamma=1$. For each stable tissue $\mathbf{x}_{i}$ we compute
its energy $E=H(\mathbf{x}_{i})$ and the number of unique single-cell states
$n(\mathbf{x}_{i})$. Each panel visualizes the joint distribution of these
quantities as well as their marginal distributions $p(n)$, $p(E)$. The dashed
green (gold) box denotes the compositionally simple, low energy (complex, high
energy) segments of the distribution. The top-left panel $\mathbf{W}_{1}$
corresponds to the cell-cell signaling rules used in Figs. 4 and 5. All panels
use the same $\mathbf{J}$, $\mathbf{A}$ as in Figs. 4 and 5.

Despite some expected variation between the eight plots, the three features
identified in the preceding subsection are largely preserved. All plots
exhibit a positive correlation between the energy $E$ and the number of unique
single-cell states $n$. This indicates that the deepest minima in the energy
landscape Eq. (5) are also the simplest (i.e. composed of few types of cells),
which is denoted in Fig. 6 by the dashed green box. Likewise, the more complex
minima (large $n$) tend to have higher energies and form a “bulk” (indicated
by a dashed gold box in Fig. 6) which also has a positive slope. In all plots
we see that these simpler, deep minima appear with much higher probability
(note the logarithmic scale) than the more complex, shallow minima.

This data suggests that, in the strong interaction regime, the energy
landscape exhibits several universal features which are relatively invariant
under different realizations of the cell-cell interaction matrix $\mathbf{W}$.
Minor deviations arising from atypical $\mathbf{W}$ are observed, (e.g. the
lack of a strong single peak at low $n$ for $\mathbf{W}_{5}$). Most notably,
the lowest energy minima tend to also be the simplest minima, as quantified by
$n(\mathbf{x})$. Furthermore, these deep, compositionally simple minima appear
to have large basins of attraction when compared to the “bulk” of more
disordered, shallow minima.

##  IV Summary and Discussion

Our work is motivated by a question that is central to the organization of
multicellular life: how do tissues self-assemble into the numerous spatial
architectures with diverse phenotypic compositions that make up multicellular
organisms from the building blocks of genetically identical single cells? How
can the cellular phenotypes on the one hand be sufficiently plastic and
adaptable to generate this tissue diversity, while on the other remaining
sufficiently stable to the dynamic signals found within different tissue
microenvironments? While natural gene-gene interaction networks frequently
appear disordered, with many weakly interacting genes and cross-wired
signaling pathways, the resulting tissues are comparatively simple, often
being composed of a handful of tissue-specific cell types that represent a
small fraction of the possible phenotypic richness. The high-dimensionality of
the interaction networks, in addition to their variability across different
organisms, have posed longstanding challenges for studying questions
concerning tissue self-organization and identifying general principles.

Towards addressing these challenges, we have developed a generalized model of
multicellular gene expression that couples single-cell gene regulation with
cell-cell signaling in a manner that is scalable to large numbers of genes and
cell types while being highly aligned to modern data modalities including
scRNA-seq and spatial transcriptomics. Our model allows us to systematically
study how tuning the degree of cell-to-cell coupling can control self-
organization of multicellular collectives into different classes of stable
tissue types with distinct spatial and compositional architecture. The basic
unit of the model, the single cell, is parameterized by a set of gene
expression vectors $\\{\bm{\xi}^{\mu}\\}$ corresponding to single-cell types.
These cell types are encoded as stable fixed points of a Hopfield network
(HN), defining an effective gene-gene interaction matrix $\mathbf{J}$. In
multicellular tissues described by a spatial adjacency matrix $\mathbf{A}$,
genetically identical cells with the same internal rules $\mathbf{J}$ interact
according to a cell-cell signaling matrix $\gamma\mathbf{W}$ which couples the
gene expression of adjacent cells. For a given set of regulatory rules
$(\mathbf{J},\mathbf{W})$, collective behavior emerges when the cell-cell
signaling strength $\gamma$ is increased beyond a critical value $\gamma_{c}$,
characterized by a cascade of collective transitions between different self-
organizing tissue states which are generally unstable in the absence of cell-
cell signaling. The punctuated nature of these transitions explains how
flexible, phenotypically plastic cells can robustly assemble into a wealth of
tissue configurations to carry out different functions when the strength of
collective signaling is varied, as it does during development, wound healing,
and the immune response.

To emulate maximally disordered versions of the cross-wired signaling pathways
often found in nature, we considered random cell-cell signaling networks
$\mathbf{W}$. For strong signaling, where collective effects dominate, we
found that disordered interactions tend to stabilize a surprisingly small
number of tissue “types” defined by large basins of attraction which can be
locally rugged with many subtly distinct local minima. Interestingly, several
features of the energy landscape appear invariant to different realizations of
the random signaling rules. In particular, we find that minima depth and
simplicity (as characterized by number of unique single-cell states and
spatial structure) are positively correlated. Additionally, these deep, simple
minima tend to have very large basins of attraction compared to the more
complex, shallow minima. This contraction of the transcriptional landscape is
consistent with the observation that multicellular tissues in different
organisms generally contain only a small fraction of the possible phenotypic
richness, often displaying just a handful of specialized cell types in each
tissue, even when the microscopic interactions which drive the system –
defined here by the cell-cell signaling networks – being highly disordered.
These observations may have implications in other contexts, such as
hierarchical self-organization in neuroscience [74, 75] as well as non-
biological systems, for instance material self-assembly [76].

Another important motivation for our work is the increasingly observed
plasticity of cellular phenotypes across diverse tissue contexts. Certain
classes of mammalian cells such as fibroblasts and macrophages exhibit tissue-
specific variations in addition to shared subtypes found throughout the body
[77]. Tissue-specific cell type plasticity and abnormal changes in the cell
type composition of tissue also play important roles in disease pathogenesis,
including cancer initiation and maintenance [10, 7]. These observations of
cross-tissue cell type heterogeneity raise numerous questions: When a cell
type is ectopically placed in a “sister” tissue rather than its native tissue,
will it morph to its “sister” cell type variant, and what compositional
changes does one expect in the surrounding tissue? Such questions are
difficult to address in vivo or with existing modeling approaches but are
highly suited to the model introduced here, which integrates the concept of
cell type plasticity in diverse tissue contexts with tissue self-organization
and multistability.

Accordingly, our work has important implications for the definition of “cell
type”. While the general concept is central to ongoing efforts to understand
and ultimately treat disease by cataloging the gene expression of human
tissues [19, 20], the exact definition is still under debate [78]. A common
approach to define cell type is based on the identification of clusters of
similar cells in single-cell gene expression data. A more nuanced view is
suggested by our framework, which reveals how the fixed points of the
multicellular system (stable cell states within tissue) evolve with the
strength of the intercellular signaling $\gamma$. At low $\gamma$ the
prescribed single-cell types remain stable, but as $\gamma$ increases they can
morph as multicellular structures progressively arise. In vivo, cells are
continuously interacting and these interactions control the stability of
different cell states in tissue. This effect – which should be taken into
account by emerging definitions of cell type – is formalized within our model.
By disentangling the intra- and intercellular levels of gene regulation, our
work provides a quantitative picture of cell type plasticity in which gene
regulatory fixed points shift and bifurcate upon exposure to different tissue
microenvironments.

Recent gene-free approaches which model the effects of cell-cell signaling on
a low-dimensional phenotype space [28, 29, 30, 31, 15, 16] provide an elegant
method for describing cell fate transitions. Such approaches rely on
bifurcation analysis of proximal cell types, and thus far have been limited to
a handful of related cell types represented as phenotypic fixed points. Our
framework is able to accommodate many disparate cell and tissue types while
using a less abstract (but much more detailed) state space that represents the
expression state of individual genes in each cell. Towards bridging the two
modeling techniques, we note that the presented model is naturally equipped
with a low-dimensional phenotype space: the projection of the single-cell
states onto the embedded cell types [35, 43, 48, 37, 38]. This may be an
avenue to develop a theoretical description of cell fate transitions in our
framework while improving the interpretability of the low-dimensional
phenotype space.

The current theoretical work provides the basis for future applications to
specific experimental systems. Inspired by the past successes of other
modeling approaches [28, 30, 29], this will involve experimentally informed
(non-random) $\mathbf{W}$, tissue-specific choices of adjacency matrix
$\mathbf{A}$, and gene expression noise. More broadly, our framework is
designed with emerging sequencing techniques in mind, specifically spatial
transcriptomics [79]. By using many snapshots of gene expression in tissue,
$\mathbf{J}$, $\mathbf{W}$, $\mathbf{A}$, and $\gamma$ can be conditionally
inferred (or at least constrained). Towards overcoming potential inference
challenges related to undersampling, we note that the present model and future
generalizations may be represented as graph neural networks (GNNs) [80, 81],
and the rapidly progressing GNN literature may reveal efficient inference
techniques. Future work in this direction will serve as an important test of
the framework’s experimental predictions.

There are several potential next steps to deepen the theoretical understanding
of the presented framework and extend its applicability to broader biological
contexts. First, a theoretical explanation should be developed to understand
the rich sequence of tissue transitions at intermediate values of $\gamma$
(which will depend on the detailed structure of $\mathbf{J}$ and
$\mathbf{W}$), as well as the invariant statistical properties for the minima
distributions in the high-$\gamma$ limit. This will serve to generalize our
results which are based on finite systems with $N=9$ genes and $p=3$ encoded
cell types. Relatedly, it would be useful for engineering purposes to be able
to predict the types of tissues than can self-organize (e.g. labyrinthine or
homogeneous) given arbitrary $\mathbf{J}$, $\mathbf{W}$, $\mathbf{A}$, and
$\gamma$. Second, for the random cell-cell signaling we have considered, the
stable tissue phenotypes at high $\gamma$ do not necessarily contain the
encoded cell types; the case of non-random cell-cell interactions matrices
$\mathbf{W}$ which preserve the encoded single-cell types
$\\{\bm{\xi}^{\mu}\\}$ will be investigated to better understand cellular
plasticity within tissue. Finally, cellular division and apoptosis events
occur when particular gene expression patterns are reached; incorporating
these events in the model will provide a natural framework to study the
sequences of self-organized tissue transitions which characterize development
[82]. The intra- and intercellular regulatory rules $\mathbf{J}$, $\mathbf{W}$
will determine the fate of these unfolding trajectories.

###### Acknowledgements.

We acknowledge the support of the Natural Sciences and Engineering Research
Council of Canada (NSERC) through Discovery Grant RGPIN 402591 to A.Z.; CGS-D
Graduate Fellowship to M.S. We thank Duncan Kirby and Jeremy Rothschild for
useful discussions, and Tatyana Gavrilchenko and Pearson Miller for helpful
comments on the manuscript.

## References

  * Takahashi and Yamanaka [2006] K. Takahashi and S. Yamanaka, Induction of Pluripotent Stem Cells from Mouse Embryonic and Adult Fibroblast Cultures by Defined Factors, Cell 126, 663 (2006).
  * Xu  _et al._ [2015] J. Xu, Y. Du, and H. Deng, Direct lineage reprogramming: Strategies, mechanisms, and applications, Cell Stem Cell 16, 119 (2015).
  * de Sousa e Melo and de Sauvage [2019] F. de Sousa e Melo and F. J. de Sauvage, Cellular Plasticity in Intestinal Homeostasis and Disease, Cell Stem Cell 24, 54 (2019).
  * Guerrero-Juarez  _et al._ [2019] C. F. Guerrero-Juarez, P. H. Dedhia, S. Jin, R. Ruiz-Vega, D. Ma, Y. Liu, K. Yamaga, O. Shestova, D. L. Gay, Z. Yang, K. Kessenbrock, Q. Nie, W. S. Pear, G. Cotsarelis, and M. V. Plikus, Single-cell analysis reveals fibroblast heterogeneity and myeloid-derived adipocyte progenitors in murine skin wounds, Nature Communications 10, 10.1038/s41467-018-08247-x (2019).
  * Plikus  _et al._ [2017] M. V. Plikus, C. F. Guerrero-Juarez, M. Ito, Y. R. Li, P. H. Dedhia, Y. Zheng, M. Shao, D. L. Gay, R. Ramos, T.-C. Hsi, J. W. Oh, X. Wang, A. Ramirez, S. E. Konopelski, A. Elzein, A. Wang, R. J. Supapannachart, H.-L. Lee, C. H. Lim, A. Nace, A. Guo, E. Treffeisen, T. Andl, R. N. Ramirez, R. Murad, S. Offermanns, D. Metzger, P. Chambon, A. D. Widgerow, T.-L. Tuan, A. Mortazavi, R. K. Gupta, B. A. Hamilton, S. E. Millar, P. Seale, W. S. Pear, M. A. Lazar, and G. Cotsarelis, Regeneration of fat cells from myofibroblasts during wound healing, Science 355, 748 (2017), https://www.science.org/doi/pdf/10.1126/science.aai8792 .
  * Sinha  _et al._ [2018] M. Sinha, C. K. Sen, K. Singh, A. Das, S. Ghatak, B. Rhea, B. Blackstone, H. M. Powell, S. Khanna, and S. Roy, Direct conversion of injury-site myeloid cells to fibroblast-like cells of granulation tissue, Nature Communications 9, 1 (2018).
  * Giroux and Rustgi [2017] V. Giroux and A. K. Rustgi, Metaplasia: Tissue injury adaptation and a precursor to the dysplasia-cancer sequence, Nature Reviews Cancer 17, 594 (2017), arXiv:15334406 .
  * Means  _et al._ [2005] A. L. Means, I. M. Meszoely, K. Suzuki, Y. Miyamoto, A. K. Rustgi, R. J. Coffey, C. V. Wright, D. A. Stoffers, and S. D. Leach, Pancreatic epithelial plasticity mediated by acinar cell transdifferentiation and generation of nestin-positive intermediates, Development 132, 3767 (2005).
  * Smart  _et al._ [2021] M. Smart, S. Goyal, and A. Zilman, Roles of phenotypic heterogeneity and microenvironment feedback in early tumor development, Physical Review E 103, 032407 (2021).
  * Yuan  _et al._ [2019] S. Yuan, R. J. Norgard, and B. Z. Stanger, Cellular plasticity in cancer, Cancer Discovery 9, 837 (2019).
  * Waddington [1957] C. H. Waddington,  _The strategy of the genes. A discussion of some aspects of theoretical biology._ (Allen & Unwin, London, 1957) p. 274 pp.
  * Chang  _et al._ [2008] H. H. Chang, M. Hemberg, M. Barahona, D. E. Ingber, and S. Huang, Transcriptome-wide noise controls lineage choice in mammalian progenitor cells, Nature 453, 544 (2008).
  * Huang [2012] S. Huang, The molecular and mathematical basis of Waddington’s epigenetic landscape: A framework for post-Darwinian biology?, BioEssays 34, 149 (2012).
  * Huang  _et al._ [2005] S. Huang, G. Eichler, Y. Bar-Yam, and D. E. Ingber, Cell fates as high-dimensional attractor states of a complex gene regulatory network, Physical Review Letters 94, 1 (2005).
  * Rand  _et al._ [2021] D. A. Rand, A. Raju, M. Sáez, F. Corson, and E. D. Siggia, Geometry of gene regulatory dynamics, Proceedings of the National Academy of Sciences of the United States of America 118, 10.1073/pnas.2109729118 (2021), arXiv:2105.13722 .
  * Sáez  _et al._ [2021] M. Sáez, R. Blassberg, E. Camacho-Aguilar, E. D. Siggia, D. A. Rand, and J. Briscoe, Statistically derived geometrical landscapes capture principles of decision-making dynamics during cell fate transitions, Cell Systems 10.1016/j.cels.2021.08.013 (2021).
  * Wang  _et al._ [2011] J. Wang, K. Zhang, L. Xu, and E. Wang, Quantifying the Waddington landscape and biological paths for development and differentiation, Proceedings of the National Academy of Sciences 108, 8257 (2011).
  * Han  _et al._ [2018] X. Han, R. Wang, Y. Zhou, L. Fei, H. Sun, S. Lai, A. Saadatpour, Z. Zhou, H. Chen, F. Ye, D. Huang, Y. Xu, W. Huang, M. Jiang, X. Jiang, J. Mao, Y. Chen, C. Lu, J. Xie, Q. Fang, Y. Wang, R. Yue, T. Li, H. Huang, S. H. Orkin, G. C. Yuan, M. Chen, and G. Guo, Mapping the Mouse Cell Atlas by Microwell-Seq, Cell 172, 1091 (2018).
  * Karlsson  _et al._ [2021] M. Karlsson, C. Zhang, L. Méar, W. Zhong, A. Digre, B. Katona, E. Sjöstedt, L. Butler, J. Odeberg, P. Dusart, F. Edfors, P. Oksvold, K. von Feilitzen, M. Zwahlen, M. Arif, O. Altay, X. Li, M. Ozcan, A. Mardonoglu, L. Fagerberg, J. Mulder, Y. Luo, F. Ponten, M. Uhlén, and C. Lindskog, A single–cell type transcriptomics map of human tissues, Science Advances 7, 10.1126/sciadv.abh2169 (2021).
  * Regev  _et al._ [2017] A. Regev, S. Teichmann, E. Lander, I. Amit, C. Benoist, E. Birney, B. Bodenmiller, P. Campbell, P. Carninci, and W. Enard, Science Forum: The Human Cell Atlas, eLife , 1 (2017).
  * Rozenblatt-Rosen  _et al._ [2017] O. Rozenblatt-Rosen, M. J. Stubbington, A. Regev, and S. A. Teichmann, The Human Cell Atlas: from vision to reality, Nature 550, 451 (2017), arXiv:121202 [10.1101] .
  * Jang  _et al._ [2017] S. Jang, S. Choubey, L. Furchtgott, L. N. Zou, A. Doyle, V. Menon, E. B. Loew, A. R. Krostag, R. A. Martinez, L. Madisen, B. P. Levi, and S. Ramanathan, Dynamics of embryonic stem cell differentiation inferred from single-cell transcriptomics show a series of transitions through discrete cell states, eLife 6, 10.7554/eLife.20487 (2017).
  * Turing [1952] A. M. Turing, The chemical basis of morphogenesis, Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 237, 37 (1952).
  * Murray [2002] J. D. Murray,  _Mathematical Biology I: An Introduction_, Vol. 17 (Springer, New York, 2002) p. 551.
  * Dang  _et al._ [2020] Y. Dang, D. A. Grundel, and H. Youk, Cellular Dialogues: Cell-Cell Communication through Diffusible Molecules Yields Dynamic Spatial Patterns, Cell Systems 10, 82 (2020).
  * Maire and Youk [2015] T. Maire and H. Youk, Molecular-Level Tuning of Cellular Autonomy Controls the Collective Behaviors of Cell Populations, Cell Systems 1, 349 (2015).
  * Olimpio  _et al._ [2018] E. P. Olimpio, Y. Dang, and H. Youk, Statistical Dynamics of Spatial-Order Formation by Communicating Cells, iScience 2, 27 (2018), arXiv:1706.06481 .
  * Camacho-Aguilar  _et al._ [2021] E. Camacho-Aguilar, A. Warmflash, and D. A. Rand, Quantifying cell transitions in C. elegans with data-fitted landscape models, PLoS Computational Biology 17, 10.1371/journal.pcbi.1009034 (2021).
  * Corson and Siggia [2017] F. Corson and E. D. Siggia, Gene-free methodology for cell fate dynamics during development, eLife 6, 10.7554/eLife.30743 (2017).
  * Corson and Siggia [2012] F. Corson and E. D. Siggia, Geometry, epistasis, and developmental patterning, Proceedings of the National Academy of Sciences of the United States of America 109, 5568 (2012).
  * Corson  _et al._ [2017] F. Corson, L. Couturier, H. Rouault, K. Mazouni, and F. Schweisguth, Self-organized Notch dynamics generate stereotyped sensory organ patterns in Drosophila, Science 356, 10.1126/science.aai7407 (2017).
  * Klipp  _et al._ [2011] E. Klipp, W. Liebermeister, C. Wierling, A. Kowald, H. Lehrach, and R. Herwig,  _Systems Biology: A Textbook_ , 2nd ed. (Wiley-Blackwell, 2011).
  * Osborne  _et al._ [2017] J. M. Osborne, A. G. Fletcher, J. M. Pitt-Francis, P. K. Maini, and D. J. Gavaghan, Comparing individual-based approaches to modelling the self-organization of multicellular tissues, PLoS Computational Biology 13, 10.1371/journal.pcbi.1005387 (2017).
  * Babu  _et al._ [2004] M. M. Babu, N. M. Luscombe, L. Aravind, M. Gerstein, and S. A. Teichmann, Structure and evolution of transcriptional regulatory networks, Current Opinion in Structural Biology 14, 283 (2004).
  * Amit [1989] D. J. Amit,  _Modeling Brain Function: The World of Attractor Neural Networks_ (Cambridge University Press, 1989).
  * Hopfield [1982] J. J. Hopfield, Neural networks and physical systems with emergent collective computational abilities., Proceedings of the National Academy of Sciences of the United States of America 79, 2554 (1982).
  * Lang  _et al._ [2014] A. H. Lang, H. Li, J. J. Collins, and P. Mehta, Epigenetic Landscapes Explain Partially Reprogrammed Cells and Identify Key Reprogramming Genes, PLoS Computational Biology 10, 10.1371/journal.pcbi.1003734 (2014), arXiv:1211.3133 .
  * Pusuluri  _et al._ [2018] S. T. Pusuluri, A. H. Lang, P. Mehta, and H. E. Castillo, Cellular reprogramming dynamics follow a simple 1D reaction coordinate, Physical Biology 15, 10.1088/1478-3975/aa90e0 (2018).
  * Levchenko [2003] A. Levchenko, Dynamical and Integrative Cell Signaling: Challenges for the New Biology, Biotechnology and Bioengineering 84, 773 (2003).
  * Kauffman [1969] S. A. Kauffman, Metabolic stability and epigenesis in randomly constructed genetic nets, Journal of Theoretical Biology 22, 437 (1969).
  * Mojtahedi  _et al._ [2016] M. Mojtahedi, A. Skupin, J. Zhou, I. G. Castaño, R. Y. Leong-Quong, H. Chang, K. Trachana, A. Giuliani, and S. Huang, Cell Fate Decision as High-Dimensional Critical State Transition, PLoS Biology 14, 10.1371/journal.pbio.2000640 (2016).
  * Amit  _et al._ [1985] D. J. Amit, H. Gutfreund, and H. Sompolinsky, Spin-glass models of neural networks, Physical Review A 32, 1007 (1985).
  * Kanter and Sompolinsky [1987] I. Kanter and H. Sompolinsky, Associative recall of memory without errors, Physical Review A 35, 380 (1987).
  * Personnaz  _et al._ [1986] L. Personnaz, I. Guyon, and G. Dreyfus, Collective computational properties of neural networks: New learning mechanisms, Physical Review A 34, 4217 (1986).
  * Fard  _et al._ [2016] A. T. Fard, S. Srihari, J. C. Mar, and M. A. Ragan, Not just a colourful metaphor: Modelling the landscape of cellular development using hopfield networks, npj Systems Biology and Applications 2, 10.1038/npjsba.2016.1 (2016).
  * Fard and Ragan [2017] A. T. Fard and M. A. Ragan, Modeling the attractor landscape of disease progression: A network-based approach, Frontiers in Genetics 8, 1 (2017).
  * Guo and Zheng [2017] J. Guo and J. Zheng, HopLand: Single-cell pseudotime recovery using continuous Hopfield network-based modeling of Waddington’s epigenetic landscape, in  _Bioinformatics_, Vol. 33 (2017) pp. i102–i109.
  * Krämer [2021] A. Krämer, Master regulators as order parameters of gene expression states, Physical Review E 103, 10.1103/PhysRevE.103.012409 (2021), arXiv:2011.06478 .
  * Maetschke and Ragan [2014] S. R. Maetschke and M. A. Ragan, Characterizing cancer subtypes as attractors of Hopfield networks, Bioinformatics 30, 1273 (2014).
  * Szedlak  _et al._ [2017] A. Szedlak, S. Sims, N. Smith, G. Paternostro, and C. Piermarocchi, Cell cycle time series gene expression data encoded as cyclic attractors in Hopfield systems, PLoS Computational Biology 13, 1 (2017), arXiv:1707.07707 .
  * Barra  _et al._ [2011] A. Barra, G. Genovese, and F. Guerra, Equilibrium statistical mechanics of bipartite spin systems, Journal of Physics A: Mathematical and Theoretical 44, 1 (2011), arXiv:1012.1261 .
  * Fedele and Contucci [2011] M. Fedele and P. Contucci, Scaling Limits for Multi-species Statistical Mechanics Mean-Field Models, Journal of Statistical Physics 144, 1186 (2011).
  * Kirsch and Toth [2020] W. Kirsch and G. Toth, Two Groups in a Curie–Weiss Model with Heterogeneous Coupling, Journal of Theoretical Probability 33, 2001 (2020).
  * Knöpfel  _et al._ [2020] H. Knöpfel, M. Löwe, K. Schubert, and A. Sinulis, Fluctuation Results for General Block Spin Ising Models, Journal of Statistical Physics 178, 1175 (2020).
  * Agliari  _et al._ [2018] E. Agliari, D. Migliozzi, and D. Tantari, Non-convex Multi-species Hopfield Models, Journal of Statistical Physics 172, 1247 (2018), arXiv:1807.03609 .
  * Barra  _et al._ [2015] A. Barra, P. Contucci, E. Mingione, and D. Tantari, Multi-Species Mean Field Spin Glasses. Rigorous Results, Annales Henri Poincare 16, 691 (2015).
  * Panchenko [2015] D. Panchenko, The free energy in a multi-species Sherrington-Kirkpatrick model, Annals of Probability 43, 3494 (2015), arXiv:arXiv:1310.6679v2 .
  * Kaneko [1992] K. Kaneko, Overview of coupled map lattices, Chaos 2, 279 (1992).
  * Bignone [1993] F. A. Bignone, Cells-gene interactions simulation on a coupled map lattice, Journal of Theoretical Biology 161, 231 (1993).
  * García-Morales  _et al._ [2017] V. García-Morales, J. A. Manzanares, and S. Mafe, Weakly coupled map lattice models for multicellular patterning and collective normalization of abnormal single-cell states, Physical Review E 95, 1 (2017), arXiv:1612.05006 .
  * Klevecz [1998] R. R. Klevecz, Phenotypic heterogeneity and genotypic instability in coupled cellular arrays, Physica D: Nonlinear Phenomena 124, 1 (1998).
  * Glauber [1963] R. J. Glauber, Time-dependent statistics of the Ising model, Journal of Mathematical Physics 4, 294 (1963).
  * Hopfield [1984] J. J. Hopfield, Neurons with graded response have collective computational properties like those of two-state neurons, Biophysics 81, 3088 (1984).
  * ichi Funahashi and Nakamura [1993] K. ichi Funahashi and Y. Nakamura, Approximation of dynamical systems by continuous time recurrent neural networks, Neural Networks 6, 801 (1993).
  * Aguilera  _et al._ [2021] M. Aguilera, S. A. Moosavi, and H. Shimazaki, A unifying framework for mean-field theories of asymmetric kinetic Ising systems, Nature Communications 12, 10.1038/s41467-021-20890-5 (2021), arXiv:2002.04309 .
  * Mézard and Sakellariou [2011] M. Mézard and J. Sakellariou, Exact mean-field inference in asymmetric kinetic Ising systems, Journal of Statistical Mechanics: Theory and Experiment 2011, 10.1088/1742-5468/2011/07/L07001 (2011), arXiv:1103.3433 .
  * Roudi  _et al._ [2009] Y. Roudi, E. Aurell, and J. A. Hertz, Statistical physics of pairwise probability models, Frontiers in Computational Neuroscience 3, 10.3389/neuro.10.022.2009 (2009).
  * Lo Cicero  _et al._ [2015] A. Lo Cicero, P. D. Stahl, and G. Raposo, Extracellular vesicles shuffling intercellular messages: For good or for bad, Current Opinion in Cell Biology 35, 69 (2015).
  * McInnes  _et al._ [2018] L. McInnes, J. Healy, and J. Melville, UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction, arXiv:1802.03426v3 (2018), arXiv:1802.03426 .
  * Löpez  _et al._ [2019] C. A. Löpez, V. V. Vesselinov, S. Gnanakaran, and B. S. Alexandrov, Unsupervised Machine Learning for Analysis of Phase Separation in Ternary Lipid Mixture, Journal of Chemical Theory and Computation 15, 6343 (2019).
  * Mehta  _et al._ [2019] P. Mehta, M. Bukov, C. H. Wang, A. G. Day, C. Richardson, C. K. Fisher, and D. J. Schwab, A high-bias, low-variance introduction to Machine Learning for physicists, Physics Reports 810, 1 (2019), arXiv:1803.08823 .
  * Wang and Ferguson [2018] J. Wang and A. L. Ferguson, Nonlinear machine learning in simulations of soft and biological materials, Molecular Simulation 44, 1090 (2018).
  * Shires and Pickard [2021] B. W. Shires and C. J. Pickard, Visualizing Energy Landscapes through Manifold Learning, Physical Review X 11, 10.1103/PhysRevX.11.041026 (2021), arXiv:2111.07843 .
  * Naim  _et al._ [2018] M. Naim, V. Boboeva, C. J. Kang, and A. Treves, Reducing a cortical network to a Potts model yields storage capacity estimates, Journal of Statistical Mechanics: Theory and Experiment 2018, 10.1088/1742-5468/aab683 (2018), arXiv:1710.04897 .
  * O’Kane and Treves [1992] D. O’Kane and A. Treves, Short- and long-range connections in autoassociative memory, Journal of Physics A: Mathematical and General 25, 5055 (1992).
  * Murugan  _et al._ [2015] A. Murugan, Z. Zeravcic, M. P. Brenner, and S. Leibler, Multifarious assembly mixtures: Systems allowing retrieval of diverse stored structures, Proceedings of the National Academy of Sciences of the United States of America 112, 54 (2015), arXiv:1408.6893 .
  * Elmentaite  _et al._ [2022] R. Elmentaite, C. Domínguez Conde, L. Yang, and S. A. Teichmann, Single-cell atlases: shared and tissue-specific cell types across human organs, Nature Reviews Genetics 10.1038/s41576-022-00449-w (2022).
  * Cel [2017] What Is Your Conceptual Definition of “Cell Type” in the Context of a Mature Organism?, Cell Systems 4, 255 (2017).
  * Longo  _et al._ [2021] S. K. Longo, M. G. Guo, A. L. Ji, and P. A. Khavari, Integrating single-cell and spatial transcriptomics to elucidate intercellular tissue dynamics, Nature Reviews Genetics 22, 627 (2021).
  * Scarselli  _et al._ [2009] F. Scarselli, M. Gori, A. C. Tsoi, M. Hagenbuchner, and G. Monfardini, The graph neural network model, IEEE Transactions on Neural Networks 20, 61 (2009).
  * Zhou  _et al._ [2020] J. Zhou, G. Cui, S. Hu, Z. Zhang, C. Yang, Z. Liu, L. Wang, C. Li, and M. Sun, Graph neural networks: A review of methods and applications, AI Open 1, 57 (2020), arXiv:1812.08434 .
  * Stanoev  _et al._ [2021] A. Stanoev, C. Schröter, and A. Koseska, Robustness and timing of cellular differentiation through population-based symmetry breaking, Development (Cambridge) 148, 10.1242/dev.197608 (2021).