Model-Based Prediction of an Effective Adhesion Parameter Guiding Multi-Type Cell Segregation

Abstract

The process of cell-sorting is essential for development and maintenance of tissues. Mathematical modeling can provide the means to analyze the consequences of different hypotheses about the underlying mechanisms. With the Differential Adhesion Hypothesis, Steinberg proposed that cell-sorting is determined by quantitative differences in cell-type-specific intercellular adhesion strengths. An implementation of the Differential Adhesion Hypothesis is the Differential Migration Model by Voss-Böhme and Deutsch. There, an effective adhesion parameter was derived analytically for systems with two cell types, which predicts the asymptotic sorting pattern. However, the existence and form of such a parameter for more than two cell types is unclear. Here, we generalize analytically the concept of an effective adhesion parameter to three and more cell types and demonstrate its existence numerically for three cell types based on in silico time-series data that is produced by a cellular-automaton implementation of the Differential Migration Model. Additionally, we classify the segregation behavior using statistical learning methods and show that the estimated effective adhesion parameter for three cell types matches our analytical prediction. Finally, we demonstrate that the effective adhesion parameter can resolve a recent dispute about the impact of interfacial adhesion, cortical tension and heterotypic repulsion on cell segregation.

1. Introduction

During the development of multicellular organisms, mechanical forces affect the internal states of cells as well as the interaction between cells, and they are, therefore, an integral part of all morphogenetic processes [1]. These forces are typically driven by molecular motors and transmitted via cytoskeleton elements and adhesion molecules within the cells and between them. Together, molecular motor driven movement and force transmission via adhesion complexes constitute two major self-organizing phenomena that drive tissue morphogenesis.

Since it constitutes a paradigmatic process for tissue morphogenesis, cell-sorting has found much attention both from experimental and theoretical research. It describes the tissue-scale segregation phenomenon that is observed when mixed heterotypic cell populations unmix into spatially confined homotypic cell clusters [2,3]. Over many decades this process has been studied in a number of in-vitro cell-sorting experiments for many different cell types [2,3,4,5,6,7]. To resolve which types of intercellular interaction guide individual mobile cells to find their homotypic neighbors, several theoretical hypotheses have been proposed and studied. Steinberg put forward the Differential Adhesion Hypothesis (DAH) which exploits the similarity of cell segregation with the unmixing of immiscible fluids such as water and oil. It states that cell-sorting is driven by the minimization of tissue surface tensions which result from quantitative differences in the strengths of cell-type specific intercellular adhesion [8]. However, Harris (1976) [9] remarked that tissue surface tension can result from several cellular mechanisms besides differential adhesion and suggested the differential contractility of cells as a major driver of cell-sorting. This is also acknowledged by the Differential Interfacial Tension Hypothesis (DITH) of Brodland and Chen (2000) [10]. There, the effects of intercellular adhesion and cellular contractility at cell–cell contacts are subsumed to the concept of differential interfacial tension at homo- and heterotypic cell contacts which determines the degree of mutual attachment between cells and thus drives cell-sorting. Canty et al. (2017) [11] demonstrated experimentally that repulsion at heterotypic cell–cell contacts with little to no contribution from adhesive or tensile differences between cell types can lead to cell-sorting and tissue separation. They propose the High Heterotypic Interfacial Tension Hypothesis (HIT) which asserts that heterotypic repulsion creates a situation where tension is strongly increased at heterotypic contacts compared to homotypic contacts and constitutes the major driver of cell-sorting and tissue separation.

The competing explanatory models DAH, DITH and HIT are unified within the Differential Migration Model (DMM) of Voss–Böhme and Deutsch (2010) [12]. They established that all the above hypotheses have in common that cell–cell contacts affect type-specifically the migration properties of the involved cells and parametrized this effect, independent of its specific nature - adhesion or repulsion, cortical tension or even signaling cascades. For the segregation of two cell types, the existence and form of an effective adhesion parameter (EAP) has been derived mathematically, which predicts the asymptotic sorting behavior [12].

However, while most previous studies focus on the simplest case of two cell types [7,13,14], the case of three and more cell types is more relevant in real biological systems. Studies of cell segregation and pattern formation with three or more cell types are still rare [15]. Although the DMM can be directly generalized to more than two cell types, it has only been studied for the minimal case of two cell types. In particular, the existence and form of the EAP for systems with three or more cell types remains an open question.

We study numerically the impact of the EAP and related parameters on the degree of tissue separation in terms of clear tissue boundaries as well as the separation time scale for mixtures of two cell types. We present new arguments to analytically predict the form of the EAP, thus generalizing the concept of the EAP to mixtures of arbitrary many cell types. We show numerically for three cell types that the analytically predicted EAP guides the asymptotic sorting behavior in the DMM, i.e., the EAP determines the asymptotic value of the normalized sum of heterotypic contacts. The form of the EAP on the asymptotic segregation behavior is independently quantified via two statistical learning methods, Support Vector Machines (SVM) and Logistic Regression model (Logit) [16]. The form estimated by the statistical learning methods matches our analytical prediction both for two and three cell types. Finally, we demonstrate that the EAP can explain different segregation behaviors in two cell- type Cellular Potts models which was previously attributed to the competing hypotheses DAH/DITH and HIT [11].

2. Materials and Methods

2.1. Differential Migration Model

The DMM is a special probabilistic cellular automaton (PCA) based on Voss-Böhme and Deutsch [12], which has been analyzed for two cell types. Here, we present the key features of this model for the case of arbitrary many cell types:

• Lattice: The PCA works on a squared lattice $S:={\{0,1,\cdots ,L\}}^2$ with periodic boundaries and we set $L=25$ throughout our numerical analysis.
• Cells: To emulate biological cells, every lattice site $x\in S$ is assigned a cell type $w\in W$ of all possible cell types $W:=\{0,1,\cdots \}$ in the system. In this way each cell regardless of its type, has approximately the same size and occupies exactly one lattice site. The number of cells of each type remains constant during the sorting process.
• Configurations: A configuration ${\eta }_t$ of the lattice S represents the state of the model at time t. It belongs to the set of all possible configurations $X:=W^S=\{\eta :S\to W\}$.
• Migration: A configuration $\eta$ is changed by a cell position switch $\eta \to {\eta }_{xy}$ involving two adjacent lattice sites x and y with $x,y\in S$:

  $${\eta }_{xy}\left(z\right):=\left\{\begin{array}{c}\eta \left(z\right),z\ne x\mathrm{and}z\ne y\hfill \\ \eta \left(y\right),z=x\hfill \\ \eta \left(x\right),z=y.\hfill \end{array}\right.$$

  (1)
  Position switches between two cells of the same type do not change the configuration and are therefore neglected. Thus, neighboring cells situated at lattice sites x and y only switch their positions if they are heterotypic, i.e., $\eta \left(x\right)\ne \eta \left(y\right)$.
• Differential adhesion: We assume that stronger bonds to neighboring cells hinder cell motility. Accordingly, a cell switch $\eta \to {\eta }_{xy}$ occurs with rate $\mathfrak{c}(x,y,\eta )$. The rate $\mathfrak{c}(x,y,\eta )$ depends on the parameters ${\beta }_{\eta \left(x\right)\eta \left(z\right)}$ and ${\beta }_{\eta \left(y\right)\eta \left(z\right)}$, which represent the binding strengths between cells at lattice sites x and y to positions from the von-Neumann-1 neighborhoods $\mathfrak{N}\left(x\right)$ and $\mathfrak{N}\left(y\right)$, see Figure 1. The von-Neumann-1 neighborhood of a lattice site corresponds to all neighboring lattice sites with Manhattan distance one. The cell switch rate of the two adjacent lattice sites $x,y\in S$ with $x\in \mathfrak{N}\left(y\right)$ and $y\in \mathfrak{N}\left(x\right)$ is as follows:

  $$\mathfrak{c}(x,y,\eta ):=\left\{\begin{array}{cc}{\displaystyle exp\left(-\sum _{z\in \mathfrak{N}\left(x\right)}{\beta }_{\eta \left(x\right)\eta \left(z\right)}-\sum _{z\in \mathfrak{N}\left(y\right)}{\beta }_{\eta \left(y\right)\eta \left(z\right)}\right),}\hfill & \eta \left(x\right)\ne \eta \left(y\right)\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

  (2)

Notice that the particular mechanism which affects the cellular motility on intercellular contact is not specified. The parameters ${\beta }_{ij}$ can be interpreted as repulsion, if ${\beta }_{ij}<0$, enhancing cellular motility, and it can be interpreted as binding strength resulting from the interplay of adhesion and relaxed cortical tension at cell–cell contact, if ${\beta }_{ij}>0$, inhibiting cellular motility. The details on the numerical implementation of this PCA are elaborated in the Appendix A. Depending on how many cell types are considered the number of intercellular adhesion parameters varies. For this, the vector ${\mathit{\beta }}_{\left(N\right)}$ of all intercellular adhesion parameters occurring in a system with N cell types is introduced, for instance

$$\begin{array}{c}{\mathit{\beta }}_{(N=2)}:={({\beta }_{00},{\beta }_{11},{\beta }_{01})}^T\mathrm{for}2-\mathrm{cell}-\mathrm{type}\mathrm{systems}\mathrm{and}\hfill \\ {\mathit{\beta }}_{(N=3)}:={({\beta }_{00},{\beta }_{11},{\beta }_{22},{\beta }_{01},{\beta }_{02},{\beta }_{12})}^T\mathrm{for}3-\mathrm{cell}-\mathrm{type}\mathrm{systems}.\hfill \end{array}$$

For a 2-cell-type model, an effective adhesion parameter (EAP)

$${\beta }_{(N=2)}^{*}:={\beta }_{00}+{\beta }_{11}-2{\beta }_{01}$$

was analytically predicted in Voss-Böhme and Deutsch [12]. The EAP determines the asymptotic cell segregation behavior of systems with $N=2$ cell types regardless of the value combinations of ${\mathit{\beta }}_{\left(2\right)}$. The higher the EAP the more segregated the cell types become.

For simulation of the DMM model a Gillespie-related algorithm is used where independent cell switches, see Equation (1), are chosen and the waiting times between these events are calculated based on the cell switch rates, see Equation (2), within lattice configuration $\eta$. For details about this algorithm, see Appendix A, Algorithm A1. Each numerical simulation yields a time-series of lattice configurations ${\eta }_t$, since every cell switch, i.e., whenever two neighboring cells switch their position, see Equation (1), results in a new configuration. The variable t refers to the number of cell switches performed.

2.2. Order Indicator

To quantify the level of cell segregation an order indicator is introduced that quantifies the level of cell segregation of a lattice configuration $\eta$. The order indicator $\omega \left(\eta \right)$ is the normalized sum of homotypic cell–cell contacts within a configuration $\eta$. This indicator increases the more cells of the same type are clustered together. In detail, the set of all homotypic von-Neumann-1 neighbors of a cell of type $\eta \left(x\right)$ at position $x\in S$ is given by

$$H(x,\eta ):=\{z\in S\mid \eta (x)=\eta (z),z\in \mathfrak{N}(x\left)\right\}.$$

With $n_{=}(x,\eta ):=\left|H(x,\eta )\right|$, the total sum of homotypic connections is defined as

$$d\left(\eta \right):={\displaystyle \frac{{\displaystyle \sum _{x\in S}{n}_{=}(x,\eta )}}{2}}.$$

(3)

The highest amount $d_{max}$ and the lowest amount $d_{min}$ of homotypic lattice site connections are reached in the case of a fully sorted and chessboard-patterned configuration, respectively. The exact values for 2- and 3-cell-type systems are listed in Appendix B. With $d_{max}$ and $d_{min}$, the order indicator value $\omega \left(\eta \right)$ for configuration $\eta$ is defined as the normalization of $d\left(\eta \right)$:

$$\omega \left(\eta \right):={\displaystyle \frac{d\left(\eta \right)-d_{min}}{d_{max}-d_{min}}}$$

Thus, $\omega \left(\eta \right)$ takes values in $[0,1]$, where configurations with almost perfect sorted patterns, such as in Figure 2I, yield $\omega \left(\eta \right)$ values above $0.95$. In the case of chessboard configurations as in Figure 2III, the order indicator values are below $0.05$. For a random lattice configuration, the order indicator value is $\omega \left(\eta \right)\approx 0.5$, in the case of a 2-cell-type system, and it is $\omega \left(\eta \right)\approx 0.33$ for a 3-cell-type system.

With the order indicator $\omega \left(\eta \right)$, a time-series of configurations ${\eta }_t$ can be converted into a time-series of order indicator values $\omega \left({\eta }_t\right)$, in short $\omega \left(t\right)$. Furthermore, the asymptotic cell segregation level for a simulation is estimated as the average over the last $10\%$ of the time-series $\omega \left(t\right)$ and is denoted $\overline{\omega }$, see Figure 2 for an illustration.

3. Results

3.1. Cell System Parameters for Two Cell Types

We present a new simple argument for the form of the effective adhesion parameter ${\beta }_{\left(2\right)}^{*}$ guiding the asymptotic behavior in the DMM for two cell types [12] to generalize it to systems with arbitrary many cell types. For this, we introduce two additional system parameters ${\beta }_{\left(2\right)}^s$ and ${\beta }_{\left(2\right)}^{\Delta }$ besides ${\beta }_{\left(2\right)}^{*}$, which determine how changes to the intercellular adhesion parameters ${\mathit{\beta }}_{\left(2\right)}$ affect the asymptotic cell segregation behavior as well as the dynamics. Although the parameters ${\beta }_{\left(2\right)}^s$ and ${\beta }_{\left(2\right)}^{*}$ set the temporal scale of the waiting times and the asymptotic level of segregation, respectively, we argue that the parameter ${\beta }_{\left(2\right)}^{\Delta }$ determines the number of cell switches required to reach the asymptotic level of segregation.

Indeed, the addition of a constant $\theta$ to all intercellular adhesion parameters simultaneously like

$${\mathit{\beta }}_{\left(2\right)}^{{}^{\prime }}={\mathit{\beta }}_{\left(2\right)}+\theta {\mathit{a}}_s\mathrm{with}{\mathit{a}}_s={(1,1,1)}^T$$

(4)

rescales the cell switch rates ${\mathfrak{c}}_{\mathit{\beta }}(x,y,\eta )$ in the DMM by a factor $exp(-8\theta )$, since

$$\begin{array}{cc}\hfill {\mathfrak{c}}_{{\mathit{\beta }}_{\left(2\right)}^{{}^{\prime }}}(x,y,\eta )& {\displaystyle =exp\left(-\sum _{z\in \mathfrak{N}\left(x\right)}({\beta }_{\eta \left(x\right)\eta \left(z\right)}+\theta )-\sum _{z\in \mathfrak{N}\left(y\right)}({\beta }_{\eta \left(y\right)\eta \left(z\right)}+\theta )\right)}\hfill \\ \hfill & =exp(-8\theta )·{\mathfrak{c}}_{{\mathit{\beta }}_{\left(2\right)}}(x,y,\eta ),\hfill \end{array}$$

(5)

independent of the actual configuration $\eta \in X$ and the considered lattice sites $x,y\in S$. This factor does not alter the relation between the cell switch rates, but only rescales all waiting times between cell switches by $exp(-8\theta )$. Thus, the spatio-temporal order of the cell sorting dynamics and the asymptotic behavior are not affected. In order to parametrize this rescaling, a scaling parameter ${\beta }_{\left(2\right)}^s$

$${\beta }_{\left(2\right)}^s={\beta }_{00}+{\beta }_{11}+{\beta }_{01}=\langle {\mathit{\beta }}_{\left(2\right)},{\mathit{a}}_s\rangle$$

(6)

is introduced, where the functional, which maps ${\mathit{\beta }}_{\left(2\right)}\to {\beta }_{\left(2\right)}^s$, is given as a scalar product with the vector ${\mathit{a}}_s={(1,1,1)}^T$.

Secondly, we notice that the choice which cell type is denoted by type $w=0$ and which one by $w=1$ is arbitrary and the dynamics in the model is invariant against relabeling of the cell types. This invariance is reflected in the cell switch rates $\mathfrak{c}(x,y,\eta )$ as well as in the order indicator $\omega$, which quantifies the level of segregation of a configuration $\eta$ as sum of homotypic lattice site connections. Accordingly, the functional which maps the intercellular adhesion parameters ${\mathit{\beta }}_{\left(2\right)}$ to the effective adhesion parameter ${\beta }_{\left(2\right)}^{*}$ must reflect this invariance as well. This implies that an effective adhesion parameter ${\beta }_{\left(2\right)}^{*}$, which controls the asymptotic sorting behavior, assuming it exists and is determined by a linear functional on the intercellular adhesion parameters, has the form

$$\begin{array}{cc}\hfill {\beta }_{\left(2\right)}^{*}& =\mathfrak{a}\left({\beta }_{00}+{\beta }_{11}\right)+\mathfrak{b}{\beta }_{01}=\langle {\mathit{\beta }}_{\left(2\right)},{\mathit{a}}_{*}\rangle \hfill \end{array}$$

(7)

with two real-valued constants $\mathfrak{a}$, $\mathfrak{b}$ and a vector ${\mathit{a}}_{*}={(\mathfrak{a},\mathfrak{a},\mathfrak{b})}^T$ introduced analogously to Equation (6) to express the functional. The constants $\mathfrak{a}$, $\mathfrak{b}$ are set by the condition that the asymptotic level of segregation and thus the differential adhesion parameter ${\beta }_{\left(2\right)}^{*}$ must not be affected by the temporal scaling parameter ${\beta }_{\left(2\right)}^s$, which implies

$$\langle {\mathit{a}}_{*},{\mathit{a}}_s\rangle =0\Rightarrow 2\mathfrak{a}=-\mathfrak{b}.$$

(8)

Without loss of generality, we can neglect an additional offset and factor in Equation (7) and choose $\mathfrak{a}=1$ for which follows

$$\begin{array}{cc}\hfill {\beta }_{\left(2\right)}^{*}& ={\beta }_{00}+{\beta }_{11}-2{\beta }_{01}=\langle {\mathit{\beta }}_{\left(2\right)},{\mathit{a}}_{*}\rangle \hfill \\ \hfill {\mathit{a}}_{*}& ={(1,1,-2)}^T.\hfill \end{array}$$

(9)

We present an additional heuristic argument to make the assumption of a linear dependency in Equation (7) plausible and to estimate that the critical value of the effective adhesion parameter, at which the system remains randomly mixed, should be zero: Consider two limit scenarios of a cell switch, first two cells at a straight interface between clusters and secondly a single cell of one type moving inside a cluster of opposite type. As illustrated in Figure 1a, two cells of different type at the interface switch their position and thus cause a less segregated configuration with a rate ${\mathfrak{c}}_{\mathrm{mix}}=exp\left[-3\left({\beta }_{00}+{\beta }_{11}\right)-2{\beta }_{01}\right]$. In contrast, a cell of type $w=0$ inside a cluster of type $w=1$ switches position with a given neighbor with rate ${\mathfrak{c}}_{\mathrm{unmix},0}=exp\left[-3{\beta }_{11}-5{\beta }_{01}\right]$, see Figure 1b. Analogously, for a single cell of type $w=1$ in a cluster of type $w=1$ the rate is ${\mathfrak{c}}_{\mathrm{unmix},1}=exp\left[-3{\beta }_{00}-5{\beta }_{01}\right]$. Both types of switches are necessary to move a single cell of opposite type out of a cluster. Since the former switch with rate ${\mathfrak{c}}_{\mathrm{mix}}$ reduces the level of segregation while the latter two switches with rates ${\mathfrak{c}}_{\mathrm{unmix},0}$, ${\mathfrak{c}}_{\mathrm{unmix},1}$ are required to increase segregation, the ratio between these types of switches ${\mathfrak{c}}_{\mathrm{unmix},0}/{\mathfrak{c}}_{\mathrm{mix}}·{\mathfrak{c}}_{\mathrm{unmix},1}/{\mathfrak{c}}_{\mathrm{mix}}=exp\left[3\left({\beta }_{00}+{\beta }_{11}-2{\beta }_{01}\right)\right]$ should determine the level of the asymptotic segregation. In particular, the point at which the system remains mixed is expected where mixing and unmixing rates are equal, i.e., where the ratio of the rates becomes 1. Since the exponential function is monotonous, it is sufficient to focus on the exponent, which is a multiple of the effective adhesion parameter ${\beta }_{\left(2\right)}^{*}={\beta }_{00}+{\beta }_{11}-2{\beta }_{01}$. Thus, the critical value 1 of the ratio of rates translates into a critical parameter ${\beta }_{\left(2\right)}^{*}=0$ for which a two-cell-type system is expected to remain mixed.

After establishing the system parameters ${\beta }_{\left(2\right)}^s$ and ${\beta }_{\left(2\right)}^{*}$, which are defined via scalar products with the corresponding vectors ${\mathit{a}}_s$, ${\mathit{a}}_{*}$, we point out that there is a third system parameter ${\beta }_{\left(2\right)}^{\Delta }$ defined by the vector perpendicular to ${\mathit{a}}_s$ and ${\mathit{a}}_{*}$

$$\begin{array}{cc}\hfill {\beta }_{\left(2\right)}^{\Delta }& =|{\beta }_{00}-{\beta }_{11}|=|\langle {\mathit{\beta }}_{\left(2\right)},{\mathit{a}}_{\Delta }\rangle |\hfill \\ \hfill {\mathit{a}}_{\Delta }& ={(1,-1,0)}^T\sim {\mathit{a}}_s\times {\mathit{a}}_{*}.\hfill \end{array}$$

(10)

The parameter ${\beta }_{\left(2\right)}^{\Delta }$ quantifies the difference between the homotypic adhesion strengths of the two cell types. It becomes zero if both types have the same strength of homotypic adhesion. We use the absolute value in the definition of the parameter ${\beta }_{\left(2\right)}^{\Delta }$ as we are not interested in which cell type has the stronger homotypic adhesion. We find numerically that the parameter ${\beta }_{\left(2\right)}^{\Delta }$ affects how many cell switches are required to reach the asymptotic level of segregation, see Figure 3. For a fixed effective adhesion parameter ${\beta }_{\left(2\right)}^{*}=3$, which means for the same asymptotic level of cell segregation, Figure 3 shows the number of cell switches required to reach certain thresholds of cell segregation, here $\omega \ge 0.7$, ≥0.8 and ≥0.9. This number increases with increasing ${\beta }_{\left(2\right)}^{\Delta }$ independently of the chosen thresholds.

This numerical observation is supported by a heuristic argument: Consider again the limit scenarios of a cell switch, i.e., two cells at a straight interface between clusters which switches with a particular neighbor with rate ${\mathfrak{c}}_{\mathrm{mix}}=exp\left[-3\left({\beta }_{00}+{\beta }_{11}\right)-2{\beta }_{01}\right]$ and a single cell of one type within a cluster of opposite type which switches with rate ${\mathfrak{c}}_{\mathrm{unmix},0}=exp\left[-3{\beta }_{11}-5{\beta }_{01}\right]$, or ${\mathfrak{c}}_{\mathrm{unmix},1}=exp\left[-3{\beta }_{00}-5{\beta }_{01}\right]$ respectively. In the symmetric case of equal homotypic adhesion parameters ${\beta }_{00}={\beta }_{11}$, which implies ${\beta }_{\left(2\right)}^{\Delta }=0$ and ${\mathfrak{c}}_{\mathrm{unmix},0}={\mathfrak{c}}_{\mathrm{unmix},1}$, cell switches of the two single cells in either cluster are equally frequent and these switches occur in a certain ratio with cell switches at the interface. If the symmetry is broken by increasing ${\beta }_{\left(2\right)}^{\Delta }$ at constant ${\beta }_{\left(2\right)}^{*}$ and ${\beta }_{\left(2\right)}^s$, for instance ${\beta }_{00}>{\beta }_{11}$, implying that ${\mathfrak{c}}_{\mathrm{unmix},0}>{\mathfrak{c}}_{\mathrm{unmix},1}$, the single cell of type $w=0$ is favored to leave the cluster of opposite type in comparison to the symmetric case and in particular compared to the other single cell of type $w=1$. Thus, segregation is now limited by the slower process of the single cell of type $w=1$ moving out of clusters of type $w=0$. Moreover, the rate ${\mathfrak{c}}_{\mathrm{mix}}$ of a cell switch at the interface is constant if only the parameter ${\beta }_{\left(2\right)}^{\Delta }$ is altered. Thus, in between two switches of a single cell of type $w=1$ (rate ${\mathfrak{c}}_{\mathrm{unmix},1}$) out of a type $w=0$ cluster there are more cell switches at the interface (rate ${\mathfrak{c}}_{\mathrm{mix}}$) than in the symmetric case. Please note that these back-and-forth switches at the interface do not progress the segregation but increase the number of total cell switches. This scenario illustrates how an increase of the parameter ${\beta }_{\left(2\right)}^{\Delta }$ increases the number of cell switches required to reach the asymptotic level of segregation. Thus, we denote ${\beta }_{\left(2\right)}^{\Delta }$ as convergence speed parameter.

Please note that the convergence speed parameter ${\beta }_{\left(2\right)}^{\Delta }$ in Equation (10) is proportional to the Euclidean distance between ${\mathit{\beta }}_{\left(2\right)}$ and the related symmetric adhesion parameters ${\mathit{\beta }}_{\left(2\right)}^{\mathrm{sym}}$ with the same values of ${\beta }_{\left(2\right)}^{*}$ and ${\beta }_{\left(2\right)}^s$. Indeed, for a given vector of adhesion parameters ${\mathit{\beta }}_{\left(2\right)}={({\beta }_{00},{\beta }_{11},{\beta }_{01})}^T$ with system parameters ${\beta }_{\left(2\right)}^s$, ${\beta }_{\left(2\right)}^{*}$, and ${\beta }_{\left(2\right)}^{\Delta }$, the corresponding symmetric adhesion parameters ${\mathit{\beta }}_{\left(2\right)}^{\mathrm{sym}}=\left(\frac{1}{2}[{\beta }_{00}+{\beta }_{11}],\frac{1}{2}[{\beta }_{00}+{\beta }_{11}],{\beta }_{01}\right)$ have the same system parameters ${\beta }_{\left(2\right)}^s$ and ${\beta }_{\left(2\right)}^{*}$ but a convergence speed parameter of zero. The Euclid distance of the adhesion parameters ${\mathit{\beta }}_{\left(2\right)}$ to this symmetric case ${\mathit{\beta }}_{\left(2\right)}^{\mathrm{sym}}$ is

$$∥{\mathit{\beta }}_{\left(2\right)}-{\mathit{\beta }}_{\left(2\right)}^{\mathrm{sym}}∥=\sqrt{\frac{1}{2}}\left|{\beta }_{00}-{\beta }_{11}\right|\sim {\beta }_{\left(2\right)}^{\Delta }.$$

(11)

Since the convergence speed parameter is zero, the symmetric case ${\mathit{\beta }}_{\left(2\right)}^{\mathrm{sym}}$ leads to the fastest convergence to the asymptotic level of segregation for the given system parameters ${\beta }_{\left(2\right)}^s$, ${\beta }_{\left(2\right)}^{*}$. In contrast to the definition of the convergence speed parameter ${\beta }_{\left(2\right)}^{\Delta }$ in Equation (10), the definition of ${\beta }_{\left(2\right)}^{\Delta }$ via Equation (11) can be generalized to N cell types, see below.

3.2. The Effective Adhesion Parameter for Arbitrary Number of Cell Types

For a mix of more than two cell types neither the form nor the existence of an effective adhesion parameter are known. Assuming that for an arbitrary number N of cell types an effective adhesion parameter ${\beta }_{\left(N\right)}^{*}$ exists, we postulate its form based on the generalization of our arguments for the case of two cell types. We also postulate the form of the convergence speed parameter ${\beta }_{\left(N\right)}^{\Delta }$ as introduced for two cell types.

The scaling parameter is directly generalized from Equation (6) to N cell types

$$\begin{array}{cc}\hfill {\beta }_{\left(N\right)}^s& =\sum _{i=0}^{N-1}{\beta }_{ii}+\sum _{i<j}^{N-1}{\beta }_{ij}=\langle {\mathit{\beta }}_{\left(N\right)},{\mathit{a}}_s\rangle \hfill \\ \hfill {\mathit{a}}_s& ={(1,\dots ,1)}^T.\hfill \end{array}$$

(12)

and only scales all waiting times between cell switches as in the case of two cell types.

The invariance against relabeling of the cell types applies as well to N cell types which implies, analogously to Equation (7),

$$\begin{array}{cc}\hfill {\beta }_{\left(N\right)}^{*}& =\mathfrak{a}\sum _{i=0}^{N-1}{\beta }_{ii}+\mathfrak{b}\sum _{i<j}^{N-1}{\beta }_{ij}=\langle {\mathit{\beta }}_{\left(N\right)},{\mathit{a}}_{*}\rangle \hfill \end{array}$$

(13)

where the constants $\mathfrak{a}$, $\mathfrak{b}$ are set again by the condition

$$\langle {\mathit{a}}_{*},{\mathit{a}}_s\rangle =0\Rightarrow N\mathfrak{a}=-\frac{N(N-1)}{2}\mathfrak{b}.$$

(14)

Choosing again $\mathfrak{a}=1$, the effective adhesion parameter has the form

$${\displaystyle {\beta }_{\left(N\right)}^{*}=\sum _{i=0}^{N-1}{\beta }_{ii}-{\displaystyle \frac{2}{(N-1)}}\sum _{i<j}^{N-1}{\beta }_{ij}}$$

(15)

with corresponding vector ${\mathit{a}}_{*}$. For $N=2$ this equals the result for two cell types from Equation (9) and for $N=3$ this takes the form

$${\beta }_{\left(3\right)}^{*}={\beta }_{00}+{\beta }_{11}+{\beta }_{22}-\left({\beta }_{01}+{\beta }_{02}+{\beta }_{12}\right).$$

(16)

In contrast to the case of two cell types, the space perpendicular to the vectors ${\mathit{a}}_s$ and ${\mathit{a}}_{*}$ for $N>2$ is not one- but $[N+N(N-1)/2-2]$-dimensional, i.e., two-dimensional for $N=3$. Thus, the definition of a convergence speed parameter ${\beta }_{\left(N\right)}^{\Delta }$, analogous to Equation (10), is ambiguous. One option is to consider all pairs of cell types and their corresponding convergence speed parameters, e.g., $|{\beta }_{00}-{\beta }_{11}|$, $|{\beta }_{00}-{\beta }_{22}|$, and $|{\beta }_{11}-{\beta }_{22}|$ for $N=3$. Since the segregation of all cell types requires the segregation of each pair of cell types, these pair-wise parameters should predict the convergence speed on subsets of intercellular adhesion parameters analogously to Figure 3. For instance, when two vectors of intercellular adhesion parameters are equal except they differ in one of these pair-wise convergence speed parameters, the one with the smaller parameter leads to a faster convergence. However, these pair-wise parameters do not allow comparison if several of them differ between two vectors of adhesion parameters and the influence of the heterotypic parameters is not even considered. Thus, the convergence speed parameter ${\beta }_{\left(N\right)}^{\Delta }$ is instead defined by the generalization of Equation (11)

$$\begin{array}{cc}\hfill {\beta }_{\left(N\right)}^{\Delta }& =∥{\mathit{\beta }}_{\left(N\right)}-{\mathit{\beta }}_{\left(N\right)}^{\mathrm{sym}}∥\hfill \\ \hfill {\mathit{\beta }}_{\left(N\right)}^{\mathrm{sym}}& =\left({\overline{\beta }}_{\mathrm{hom}},\dots ,{\overline{\beta }}_{\mathrm{hom}},{\overline{\beta }}_{\mathrm{het}},\dots ,{\overline{\beta }}_{\mathrm{het}}\right)\hfill \\ \hfill {\overline{\beta }}_{\mathrm{hom}}& =\frac{1}{N}\sum _{i=0}^{N-1}{\beta }_{ii},{\overline{\beta }}_{\mathrm{het}}=\frac{2}{N(N-1)}\sum _{i<j}^{N-1}{\beta }_{ij}.\hfill \end{array}$$

(17)

We postulate, in analogy to the case of two cell types, that the parameter ${\beta }_{\left(N\right)}^{\Delta }$ determines the number of cell switches required to reach the asymptotic level of segregation for a given effective adhesion parameter ${\beta }_{\left(N\right)}^{*}$.

3.3. Numerical Evidence for the Effective Adhesion Parameter

We test our analytical prediction for the effective adhesion parameter for more than two cell types, by comparison with numerical simulations of the DMM model for two and three cell types. We find in both cases that the proposed effective adhesion parameter determines the asymptotic level of cell segregation. However, there are deviations due to the limited simulation time, from which the asymptotic level of the order indicator $\overline{\omega }$ is extrapolated. This is revealed by additionally linking the simulation results to the convergence speed parameter ${\beta }^{\Delta }$, introduced above. For two cell types, we observe that for small values of ${\beta }_{\left(2\right)}^{\Delta }$, which refers to simulations which converge quickly to their asymptotic value, the dependency between ${\beta }_{\left(2\right)}^{*}$ and $\overline{\omega }$ agrees almost perfectly with the analytically proven prediction. Deviations from the analytical prediction are the most pronounced the larger the convergence speed parameter is. For three cell types, we observe an analogous relation between deviations from the postulated dependency on the effective adhesion parameter ${\beta }_{\left(3\right)}^{*}$ and a convergence speed parameter ${\beta }_{\left(3\right)}^{\Delta }$ as introduced above.

3.3.1. The Asymptotic Level of Cell Segregation Depends on the Effective Adhesion Parameter

At first the dependence of the level of asymptotic cell segregation $\overline{\omega }$ on the effective adhesion parameter ${\beta }^{*}$ is tested. For the 2-cell-type, the effective adhesion parameter clearly determines the asymptotic sorting state, see dark blue points in Figure 4a, i.e., $\overline{\omega }$ increases with increasing ${\beta }_{\left(2\right)}^{*}$, as predicted by Voss-Böhme and Deutsch [12], and all configurations stay randomly mixed at ${\beta }_{\left(2\right)}^{*}=0$, i.e., $\overline{\omega }({\beta }^{*}=0)=0.5$ (see black dashed lines), as predicted by the heuristic argument above. The monotonic dependence of the asymptotic order parameter $\overline{\omega }$ on the proposed effective parameter ${\beta }_{\left(2\right)}^{*}$ is supported by a value of $0.81$ of the spearman rank correlation coefficient (Kendall rank correlation coefficient of $0.58$). Additionally, Figure 4a shows the asserted influence of the convergence speed parameter ${\beta }_{\left(2\right)}^{\Delta }$ on level of cell segregation dynamics. More precisely, the higher the ${\beta }_{\left(2\right)}^{\Delta }$ value is the slower is the sorting process and the closer the respective $\overline{\omega }$ value remains to $\omega (t=0)\approx 0.5$ until the end of the simulation. For these intercellular adhesion parameters ${\mathit{\beta }}_{\left(2\right)}$, the asymptotic value of the order indicator has not yet been reached, which is why they are scattered above and below the increasing line, indicated by the dark blue points in Figure 4a. Note, that also the effective adhesion parameter ${\beta }^{*}$ shows impact on the amount of cell switches needed for a system to reach its asymptotic cell segregation $\overline{\omega }$. This is further investigated in Appendix D.

Analogous results for the case of three cell types with the postulated effective adhesion parameter are shown in Figure 4b. Qualitatively the same functional relation between the hypothesized effective adhesion parameter ${\beta }_{\left(3\right)}^{*}$ and $\overline{\omega }$ is visible as in the case of two cell types, see dark blue points in Figure 4b. This includes the prediction that all configurations stay randomly mixed at ${\beta }_{\left(3\right)}^{*}=0$, i.e., $\overline{\omega }({\beta }^{*}=0)=0.33$ (see black dashed lines). The monotonic dependence of the asymptotic order parameter $\overline{\omega }$ on the proposed effective parameter ${\beta }_{\left(3\right)}^{*}$ is supported by a value of $0.79$ of the Spearman rank correlation coefficient (Kendall rank correlation coefficient of $0.53$), which is of similar quality as in the case of two cell types. Based on the analogy of the numerical results for two and three cell types, we conclude that ${\beta }_{\left(3\right)}^{*}$ predicts the level of sorting reached asymptotically such as ${\beta }_{\left(2\right)}^{*}$ does in the 2-cell-type case. Correspondingly, the asymptotic cell segregation exhibits the postulated influence of the convergence speed parameter ${\beta }_{\left(3\right)}^{\Delta }$. Note that the restriction in Figure 4b to simulations with an ${\beta }_{\left(3\right)}^{\Delta }<3.0$ ensures a sufficient sampling of data points with smaller ${\beta }_{\left(3\right)}^{\Delta }$ parameter values. Smaller ${\beta }^{\Delta }$ values result in more simulations reaching their asymptotic cell segregation.

For the 2-cell-type system, we can additionally visualize the effective adhesion parameter ${\beta }_{\left(2\right)}^{*}$ as a plane in 3D, see Figure 5, with $E:({\beta }_{00}+{\beta }_{11})-2{\beta }_{01}=0$ for ${\beta }_{\left(2\right)}^{*}=0$ and its normal vector ${\mathit{a}}_{*}=(1,1,-2)$. This plane separates simulated data points based on their level of asymptotic cell segregation $\overline{\omega }$. This means, the direction ${\mathit{a}}_{*}$ indeed describes the direction of order progression, see color transition in Figure 5.

We present analogous analysis for 4- and 5-cell-type systems, see Appendix C, Figure A1 and find a strong support of the analytically predicted effective adhesion parameter as well.

3.3.2. Estimating the Effective Adhesion Parameter Using Statistical Learning Methods

We confirmed qualitatively the compliance between our prediction and the numerical data in the previous section. Now we confirm additionally our prediction by comparing the analytical effective adhesion parameters for two and three cell types to the corresponding estimates obtained via two statistical learning methods, Support Vector Machines (SVM) and Logistic Regression (Logit-Model) [16]. More precisely, for the general case of N cell types, the learning methods provide a hyperplane equation of the effective adhesion parameter ${\beta }_{\left(N\right)}^{*}$. Therefore, we quantify statistically the deviation from our analytical prediction. This is done for ${\beta }_{\left(2\right)}^{*}$ and ${\beta }_{\left(3\right)}^{*}$ based on the numerical data $\overline{\omega }\left(\mathit{\beta }\right)$ referred to in the previous section and displayed in Figure 4a,b, respectively. Both methods are suited for this task because they fit the linear dependency of the parameter ${\beta }^{*}$ on the intercellular adhesion parameters $\mathit{\beta }$, are well documented and directly available from libraries such as scikit-learn for Python [17], and, in contrast to some other approaches like deep learning, both quantify the factors that lead to classification.

The learning methods provide a hyperplane equation as a linear predictor to separate classified data points. Thus, the numerical data points $\overline{\omega }\left(\mathit{\beta }\right)$ are first classified based on their asymptotic order indicator $\overline{\omega }$ and a classification threshold. As thresholds, we chose the random configurations, i.e. $\overline{\omega }\ge 0.5$ for two cell types and $\overline{\omega }\ge 0.33$ for three cell types. Due to the symmetry of the effective adhesion parameter ${\beta }^{*}$, see Equations (7) and (13), the number of coefficients, which must be estimated, can be reduced to the constant $\mathfrak{a}$ for the homotypic and the constant $\mathfrak{b}$ for the heterotypic adhesion parameters. The resulting hyperplane equations for both 2-cell-type and 3-cell-type cases are presented in the headlines of

Table 1. The effective adhesion parameter is confirmed by statistical learning methods. (a) For the estimation of the effective adhesion parameter ${\beta }_{\left(2\right)}^{*}$, 2000 simulations are conducted under the same conditions as described in Figure 4a. After simulation, each data point is classified: $\overline{\omega }>0.5$ as class 1, else class $-1$. The classified data points are used by the SVM and Logit-Model to estimate the homotypic $\mathfrak{a}$ and the heterotypic $\mathfrak{b}$ coefficient as well as the intercept $\mathfrak{i}$ according to the table head. The model accuracy is tested via 5-fold cross-validation. (b) Analogous for the 3-cell-type systems the simulation conditions are described in Figure 4b. For data point classification applies: $\overline{\omega }>0.33$ as class 1, else class $-1$. For details on the estimation process as well as the model accuracy test, see Appendix E.

(a) 2-cell-type systems	$0=\mathfrak{a}({\mathit{\beta }}_{00}+{\mathit{\beta }}_{11})+\mathfrak{b}{\mathit{\beta }}_{01}+\mathfrak{i}$			Model
	$\mathfrak{a}/\mathfrak{a}$	$\mathfrak{b}/\mathfrak{a}$	$\mathfrak{i}/\mathfrak{a}$	Accuracy
SVM	1.0000	−1.9764	−0.0374	0.9965
Logit	1.0000	−1.9763	−0.0455	0.9965
Theoretical prediction	1	−2	0
(b) 3-cell-type systems	$\mathbf{0}=\mathfrak{a}({\mathit{\beta }}_{\mathbf{00}}+{\mathit{\beta }}_{\mathbf{11}}+{\mathit{\beta }}_{\mathbf{22}})+\mathfrak{b}({\mathit{\beta }}_{\mathbf{01}}+{\mathit{\beta }}_{\mathbf{02}}+{\mathit{\beta }}_{\mathbf{12}})+\mathfrak{i}$			Model
	$\mathfrak{a}/\mathfrak{a}$	$\mathfrak{b}/\mathfrak{a}$	$\mathfrak{i}/\mathfrak{a}$	Accuracy
SVM	1.0000	−0.9994	0.2835	0.9895
Logit	1.0000	−0.9935	0.2624	0.9900
Theoretical prediction	1	−1	0

a,b.

For the 2-cell-type case in Table 1a, both relative coefficients, for homotypic $\mathfrak{a}/\mathfrak{a}\approx 1.00$ and heterotypic $\mathfrak{b}/\mathfrak{a}\approx -1.99$ intercellular adhesion parameters, are close to their predicted values 1 and $-2$, respectively. The relative intercept $\mathfrak{i}/\mathfrak{a}$, i.e. the value of ${\beta }_{\left(2\right)}^{*}$ which corresponds to the chosen threshold $\overline{\omega }=0.5$, is for both learning methods near zero, as predicted by our heuristic argument. These predictions of the SVM and the Logit-Model are of high quality as indicated by model accuracy values close to one. This quantitative result of the statistical learning methods in Table 1a support the analytical prediction for ${\beta }_{\left(2\right)}^{*}={\beta }_{00}+{\beta }_{11}-2{\beta }_{01}$.

Analogous for three cell types, the results of Table 1b support the analytical results by estimating ${\beta }_{\left(3\right)}^{*}={\beta }_{00}+{\beta }_{11}+{\beta }_{22}-{\beta }_{01}-{\beta }_{02}-{\beta }_{12}$ with convincing model accuracy of about $0.99$ for both learning methods. The relative intercept $\mathfrak{i}/\mathfrak{a}$, i.e., the value of ${\beta }_{\left(3\right)}^{*}$ which corresponds to the chosen threshold $\overline{\omega }=0.33$, is for both learning methods close to zero, as predicted by our heuristic argument. We attribute the deviation of the relative intercept $\mathfrak{i}/\mathfrak{a}$ from zero as well as the lower model accuracy compared to two cell types to insufficient convergence of the order indicators as predicted by the convergence speed ${\beta }_{\left(3\right)}^{\Delta }$.

We present analogous analysis for 4- and 5-cell-type systems, see Appendix C,

Table A2. The effective adhesion parameter is confirmed by statistical learning methods for 4 and 5 cell types. (a) For the estimation of the effective adhesion parameter ${\beta }_{\left(4\right)}^{*}$, 2000 simulations are conducted under the same conditions as described in Figure A1a. After simulation, each data point is classified: $\overline{\omega }>0.25$ as class 1, else class $-1$. The classified data points are used by the SVM and Logit-Model to estimate the homotypic $\mathfrak{a}$ and the heterotypic $\mathfrak{b}$ coefficient as well as the intercept $\mathfrak{i}$ according to the table head. The model accuracy is tested via 5-fold cross-validation. (b) Analogous for the 5-cell-type systems the simulation conditions are described in Figure A1b. For data point classification applies: $\overline{\omega }>0.20$ as class 1, else class $-1$. For details on the estimation process as well as the model accuracy test, see Appendix E.

(a) 4-cell-type systems	$0=\mathfrak{a}{\displaystyle \sum _{\mathit{i}=0}^3}{\mathit{\beta }}_{\mathbf{ii}}+\mathfrak{b}{\displaystyle \sum _{\mathit{i}<\mathit{j}}^3}{\mathit{\beta }}_{\mathbf{ij}}+\mathfrak{i}$			Model
	$\mathfrak{a}/\mathfrak{a}$	$\mathfrak{b}/\mathfrak{a}$	$\mathfrak{i}/\mathfrak{a}$	Accuracy
SVM	1.0000	−0.6553	0.8022	0.9615
Logit	1.0000	−0.6625	0.7330	0.9610
Theoretical prediction	1	−2/3	0
(b) 5-cell-type systems	$\mathbf{0}=\mathfrak{a}{\displaystyle \sum _{\mathbf{i}=\mathbf{0}}^{\mathbf{4}}}{\mathit{\beta }}_{\mathbf{ii}}+\mathfrak{b}{\displaystyle \sum _{\mathbf{i}<\mathbf{j}}^{\mathbf{4}}}{\mathit{\beta }}_{\mathbf{ij}}+\mathfrak{i}$			Model
	$\mathfrak{a}/\mathfrak{a}$	$\mathfrak{b}/\mathfrak{a}$	$\mathfrak{i}/\mathfrak{a}$	Accuracy
SVM	1.0000	−0.5030	0.5889	0.9555
Logit	1.0000	−0.4913	0.5992	0.9565
Theoretical prediction	1	−0.5	0

and again find a strong support of the analytically predicted parameters defining the effective adhesion parameter.

3.4. The Impact of Interfacial Tension, Adhesion or Repulsion on Cell Segregation

For two cell types, the EAP resolves previous discussions about the impact of interfacial tension on cell segregation compared to interfacial adhesion or repulsion [11]: We predict that a higher level of segregation is reached when the EAP is large, which can be achieved by both heterotypic repulsion or differential adhesion. In fact, transforming the relative contact tensions used in ref. [11] to the corresponding EAP by identifying high tensions $T_{ij}$ in the experiments and the Cellular Potts model (CPM) with low adhesion parameters ${\beta }_{ij}=-T_{ij}$ in the DMM, we can predict the experimental and numerical observations there, see Figure 6. In particular, we find that low lengths of heterotypic interfaces (LHI) are observed for high EAP values, analogous to the DMM. Figure 6 shows that the EAP is a better predictor for the asymptotic segregation behavior than the classification according to the DAH/DITH and HIT hypotheses.

Thus, we propose that rather than differentiating between intercellular adhesion and contact tension, the combined effect of intercellular contact on cellular motility, which is quantified by the EAP, should be focused on. Furthermore, the correct prediction of the asymptotic level of segregation in the CPM simulations of Ref. [11] using the effective adhesion parameter ${\beta }^{*}$ derived for the DMM suggests that our results may be directly applicable to CPM models of cell segregation. This is plausible due to the analogous structure of the exponent of the cell switch rates in the DMM and the energy functional in the CPM.

4. Discussion

We study analytically and numerically the differential migration model of Voss-Böhme and Deutsch (2010) [12], which is a cell-based model incorporating differential hypothesis [8], differential interfacial hypothesis [10], and High Heterotypic Interfacial Tension Hypothesis [11] within a unified framework. We generalize the existence and form of an effective adhesion parameter (EAP) guiding the asymptotic level of segregation, which is already known for systems with two cell types [12], to systems with an arbitrary number of cell types. We additionally predict the critical value of this effective adhesion parameter at which the system remains randomly mixed. For the case of two and three cell types, we confirm these theoretical predictions numerically and quantify the form of the effective adhesion parameter independently using statistical learning methods. The analogous results for 4- and 5-cell types suggest that our findings are valid in systems with even higher number of cell types. For two cell types, we show that the EAP resolves previous discussions about the impact of interfacial tension on cell segregation compared to interfacial adhesion or repulsion [11]. Thus, we propose that rather than differentiating between intercellular adhesion and contact tension, the combined effect of intercellular contact on cellular motility, which is quantified by the EAP, should be focused on.

Most previous studies focus on the simplest case of two cell types [7,13,14]. Although three and more cell types are more relevant in real biological systems, studies of cell segregation and pattern formation with three or more cell types are still rare [15]. By analytically and numerically demonstrating the existence of cell segregation for three cell types and characterizing its dynamics, we extent results from two cell types to the more general case.

The differential migration model has recently been applied successfully to reproduce experimental cell segregation data for two cell types [18]. In particular, the match between model and experimental observation demonstrated that contact tension and adhesion are sufficient to explain segregation data without additional mechanisms of collective motion. Analogously, our analytical and numerical results can be applied to cell segregation experiments with three or more cell types.