Emergence of division of labor in tissues through cell interactions and spatial cues

Abstract

Most cell types in multicellular organisms can perform multiple functions. However, not all functions can be optimally performed simultaneously by the same cells. Functions incompatible at the level of individual cells can be performed at the cell population level, where cells divide labor and specialize in different functions. Division of labor can arise due to instruction by tissue environment or through self-organization. Here, we develop a computational framework to investigate the contribution of these mechanisms to division of labor within a cell-type population. By optimizing collective cellular task performance under trade-offs, we find that distinguishable expression patterns can emerge from cell-cell interactions versus instructive signals. We propose a method to construct ligand-receptor networks between specialist cells and use it to infer division-of-labor mechanisms from single-cell RNA sequencing (RNA-seq) and spatial transcriptomics data of stromal, epithelial, and immune cells. Our framework can be used to characterize the complexity of cell interactions within tissues.

Keywords: CP: Developmental biology; Pareto optimality; division of labor; enterocytes; fibroblasts; lateral inhibition; macrophages; morphogens; self-organization.

Graphical abstract

Figure 1

Pareto optimality framework (A) Single-cell RNA sequencing data in expression space can reveal various low-dimensional structured patterns. (B) Mapping of the representation of individual cells across task performance, tissue, and expression space. (C) We consider two underlying mechanisms driving the allocation of tasks in the tissue: external signaling gradients and cell-cell interactions. (D) Mapping patterns in tissue space with those in expression space is an ill-posed inverse problem—there are multiple possible tissue compositions consistent with a given expression configuration.

Figure 2

A variety of tissue and expression patterns emerge from the Pareto optimality framework with a cell-cell communication mechanism of lateral inhibition (A) Theoretical framework of Pareto optimality with cell-cell interactions—cell $c$’s performance is affected by the performance of its neighboring cells in the tissue, $N_c$. The total performance function, $F$, is a product over the tasks, where the performance in each task is the sum over the contribution of all cells considering a self-component ($P_t$) and the effect of interaction with nearby cells ($H_t$). (B) Simulation results of lateral inhibition in expression and tissue space. Varying the range of cellular interactions produces diverse tissue patterns akin to patterns observed in real tissues.

Figure 3

Distinct spatial patterns emerge from external gradients and cell-cell interactions (A and B) The Pareto-optimal solution of cells that are affected by an external gradient across a 1D and a 2D tissue in gene expression (A) and tissue (B) space. (C) With external gradients, the pairwise expression distances (y axis) versus the pairwise physical distances (x axis) show high Pearson correlation (for 1D: corr = 0.88, p < 0.001, for 2D: corr = 0.72, p < 0.001). (D and E) The Pareto-optimal solution of cells that are affected by cell-cell interactions in a 1D and a 2D tissue in gene expression (D) and tissue (E) space. (F) With cell-cell interactions, the pairwise expression distances (y axis) versus the pairwise physical distances (x axis) are anticorrelated (in terms of Pearson correlation, for 1D, corr = −0.1, p = 0.002, for 2D, corr = −0.13, p < 0.001). (G) Schematics of the colon tissue including the intestinal villi and the muscularis layer. (H) Gene expression profile of enterocytes where the cells are colored by their task specializations (combinations of red, green, and blue colors representing the three tasks). (I) Arrangements of enterocytes along the crypt-to-villus axis. The distance of the cells from each archetype is plotted as a function of the distance from the crypt. (J) The pairwise distances in expression versus physical space of enterocytes show high Pearson correlation (corr = 0.67, p < 0.001). (K) Gene expression profile of colon fibroblasts from Slide-seq data (Avraham-Davidi et al.²⁴), where the cells are colored by their task specializations. At the vertices (H,K) are ellipses that indicate STD of vertex position from bootstrapping. (L) Spatial arrangements of the fibroblasts in the colon tissue. (M) The pairwise distances in expression versus physical space of fibroblasts show negative Pearson correlation (corr = −0.1, p < 0.001).

Figure 4

Inferring archetype crosstalk networks of colon fibroblasts based on ligand-receptor interactions (A) Interactions between specialist cells are inferred from enrichments of ligands and their corresponding receptors. We use a directed edge to represent a pair of a ligand and its corresponding receptor that are enriched near each of the archetypes it connects. (B) A projection of gene expression profiles of single-cell RNA sequencing (scRNA-seq) colon fibroblasts on the first three PCs. Fibroblasts fill in a 5-vertex polytope (p < $10^{-4}$). (C) A table showing examples of ligands and their respective receptors that are enriched near the archetypes. We plot the complete archetype crosstalk network inferred for the colon fibroblasts where the thickness of each edge corresponds to the number of ligand-receptor pairs enriched between its vertex archetypes. (D and E) Using TACCO, we compute a mapping from fibroblast cells assayed by scRNA-seq and Slide-seq beads based on their expression agreement. Using the mapping, we (D) infer the scRNA-seq task components for each bead (depicted in pie charts per bead) and (E) compute the corresponding correlation of pairwise task distances versus physical distances (corr = −0.05, p < ${10}^{-3}$). (F) Projection of the scRNA-seq expression profiles onto the Slide-seq beads. To view the expression of Delta and Notch, we image their log ratio. Beads whose inferred Delta expression is greater/less than their Notch expression lean toward a turquoise/brown shade, respectively.

Figure 5

Archetype crosstalk networks can help estimate the role of cell-cell interactions in shaping task allocation even without spatial information (A) Human lung fibroblasts (from Adams et al.³³) fit in a 5-vertex polytope in expression space (p = 0.004) and show 5 archetypes that correspond to regulation of immune response, ECM degradation, protein biosynthesis and metabolism, ECM production, and myofibroblast specializations. (B) Archetype crosstalk network of lung fibroblasts. (C) Human lung macrophages (from Adams et al.³³) fit in a tetrahedron in expression space (p < 10⁻⁴) and show 4 archetypes that correspond to ECM degradation, phagocytosis, pro-inflammatory response, and metabolism specializations. (D) Archetype crosstalk network of lung macrophages. (E) Expression profile of intestine enterocytes (from Moor et al.¹²). (F) Archetype crosstalk network of intestine enterocytes.