A Voß-Böhme

A Deutsch

The cellular basis of cell sorting kinetics

Keywords: differential adhesion hypothesis, tissue surface tension, phase segregation, stochastic lattice gas, interacting particle system 2000 MSC: 92C15, 92C17, 92C37

Cell sorting is a dynamical cooperative phenomenon that is fundamental for tissue morphogenesis and tissue homeostasis. According to Steinberg's differential adhesion hypothesis, the structure of sorted cell aggregates is determined by physical characteristics of the respective tissues, the tissue surface tensions. Steinberg postulated that tissue surface tensions result from quantitative differences in intercellular adhesion. Several experiments in cell cultures as well as in developing organisms support this hypothesis.

The question of how tissue surface tension might result from differential adhesion was addressed in some theoretical models. These models describe the cellular interdependence structure once the temporal evolution has stabilized. In general, these models are capable of reproducing sorted patterns. However the model dynamics at the cellular scale are defined implicitly and are not welljustified. The precise mechanism describing how differential adhesion generates the observed sorting kinetics at the tissue level is still unclear.

It is necessary to formulate the concepts of cell level kinetics explicitly. Only then it is possible to understand the temporal development at the cellular and tissue scales. Here we argue that individual cell mobility is reduced the more the cells stick to their neighbors. We translate this assumption into a precise mathematical model which belongs to the class of stochastic interacting particle systems. Analyzing this model, we are able to predict the emergent sorting behavior at the population level. We describe qualitatively the geometry of cell segregation depending on the intercellular adhesion parameters. Furthermore, we derive a functional relationship between intercellular adhesion and surface tension and highlight the role of cell mobility in the process of sorting. We show that the interaction between the cells and the boundary of a confining vessel has a major impact on the sorting geometry.

A c c e p t e d m a n u s c r i p t 1. Introduction

Cell sorting is a biological process where heterotypic cell populations in composite aggregates segregate into spatially confined homotypic cell clusters. In a series of cell-sorting experiments, Holtfreter [START_REF] Holtfreter | Gewebeaffinität, ein Mittel der embryonalen Formbildung[END_REF] and subsequently many other biologists [for instance [START_REF] Armstrong | Cell sorting out: the self-assembly of tissues in vitro[END_REF][START_REF] Beysens | Cell sorting is analogous to phase ordering in fluids[END_REF] showed that dissociated embryonic tissue in a Petri-dish sorts spontaneously. The sorting proceeds via the coalescence of small islands into larger ones until the cells are aggregated within cell clusters that are homogeneous with respect to the cell type. In general, cells of any pair of types segregate into distinct tissues, not only cells that are in contact during normal development [START_REF] Forgacs | Biological implications of tissue viscoelasticity[END_REF].

It is widely assumed that the observed bulk behavior of cell populations in vitro results from type-dependent differences in individual cell properties. In particular, cell-type dependent disparities in the expression of molecules that regulate intercellular adhesion are held to be responsible for cell sorting. In vivo, cell segregation phenomena play a role in a range of developmental processes, most importantly tissue morphogenesis and tissue homeostasis. One would like to understand to which extent the intracellular regulation of spatiotemporal expression patterns of binding molecules guides the development and stabilization of biological tissues. This cannot be accomplished until one understands the details and mechanisms of cell sorting.

Cell sorting is commonly regarded to be the result of biological processes at different spatial scales. At the tissue level (macroscopic scale) cells of different populations separate into immiscible cell aggregates where striking similarities with the phase-separation of fluids are observed phenomenologically. Intercellular adhesion, the supposed cause of cell sorting, as well as other characteristics of individual cells such as cell migration properties are attributes at the cellular level (mesoscopic scale). They in turn result from biochemical interactions at the level of proteins (microscopic scale). Considering time-scales, cell sorting, like the related morphogenetic processes, is observed over long timescales (minutes to hours) whereas processes at cellular and lower levels typically develop over short (seconds) to medium (minute) timescales. Accordingly, several experimental and theoretical studies aim to understand cell sorting and intercellular adhesion at various spatial and temporal scales [START_REF] Armstrong | A continuum approach to modelling cell-cell adhesion[END_REF][START_REF] Forgacs | Viscoelastic properties of living embryonic tissues: a quantitative study[END_REF][START_REF] Krieg | Tensile forces govern germ-layer organization in zebrafish[END_REF][START_REF] Lecuit | Cell surface mechanics and the control of cell shape, tissue patterns and morphogenesis[END_REF][START_REF] Niessen | Cadherin-mediated cell sorting not determined by binding or adhesion specificity[END_REF][START_REF] Shi | Biophysical properties of cadherin bonds do not predict cell sorting[END_REF].

At the population level, cells in a tissue behave remarkably similar to molecules in a liquid. This was observed and formulated by Steinberg in his famous differential adhesion hypothesis (DAH) [START_REF] Steinberg | Reconstruction of tissues by dissociated cells[END_REF]. Steinberg postulates that there exist measurable physical quantities of homotypic tissues, the tissue

A c c e p t e d m a n u s c r i p t

surface tensions, which determine the geometry of sorted cell patterns. He proposed further, in analogy to the linkage of liquid surface tension with cohesive interactions between molecules in a liquid, that cell-type dependent differences in tissue surface tensions are induced by quantitative differences in the strength of intercellular adhesion. Subsequently, Steinberg, Foty and colleagues developed a concept of tissue surface tension and conducted a series of experiments to show that this quantity behaves in many ways analogously to liquid surface tension. The phenomenological similarity between tissues and fluids could thus be successfully exploited and quantified [START_REF] Beysens | Cell sorting is analogous to phase ordering in fluids[END_REF][START_REF] Forgacs | Viscoelastic properties of living embryonic tissues: a quantitative study[END_REF][START_REF] Foty | Liquid properties of embryonic tissues: Measurement of interfacial tensions[END_REF][START_REF] Steinberg | Reconstruction of tissues by dissociated cells[END_REF].

Liquids are cohesive materials with mobile subunits, the molecules, that are subjected to thermal fluctuations. Utilizing the phenomenological analogy between liquids and tissues, the DAH therefore implicitly acknowledges that there are two key mechanisms at the cell level that contribute to the observed behavior at the tissue level, cell mobility and intercellular adhesion. Cell mobility is a prerequisite for active tissue rearrangement and is thereby essential for the process of sorting. It is widely agreed, that the cells in a tissues are able to change their spatial positions [START_REF] Upadhyaya | Anomalous diffusion and non-Gaussian velocity distribution of Hydra cells in cellular aggregates[END_REF]. Regarding intercellular adhesion, it is accepted [START_REF] Wolpert | Principles of development[END_REF], that there is a complex system of attraction and repulsion mechanisms operating between cells of various types. Holtfreter coined the expression 'cell affinity' for the resulting effective attractive (or repulsive) force. In this spirit the term 'intercellular adhesion' is used here to describe the affinity of two cells of possibly different types regardless of the particular molecular basis for their interaction.

By now it is widely accepted [START_REF] Green | Sophistications of cell sorting[END_REF], that sticky and mobile cells at the mesoscopic scale produce the cell sorting which is observed at the macroscopic scale. However the precise link that connects the two scales is still unknown. It was shown experimentally, that tissue surface tension, as measured by specific surface tensiometers [START_REF] Foty | Surface tensions of embryonic tissues predict their mutual envelopment behavior[END_REF], increases linearly with the expression level of certain adhesion molecules such as cadherins [START_REF] Foty | The differential adhesion hypothesis: a direct evaluation[END_REF]. However it is still open from the experimental as well as theoretical point of view, in what way cell migration is damped by intercellular adhesion and to what extent the dynamic behavior at the cell level determines the dynamic behavior at the tissue level. In short, what is the mechanism that couples intercellular adhesion and surface tension and which role does cell migration play in this process?

If one tries to transfer the idea that tissues behave like fluids from tissue to cell-level, then it is necessary to acknowledge that there are fundamental differences in the microscopic sources of mesoscopic behavior. Liquid molecules move passively in response to thermal fluctuations. In contrast, cells move passively as well as actively. Active cell displacements are powered by the cell metabolism and are controlled by cell-intrinsic regulation processes. Considering the cohesive forces between the subunits we also find major differences. Liquid molecules are held together by van der Waals forces, hydrogen bonds and ionic interactions. The molecular basis for intercellular adhesion are specific bind-

A c c e p t e d m a n u s c r i p t

ing molecules, particularly cadherins, that interact with molecules of adjoining cells. The expression of those binding molecules is controlled by cell-intrinsic regulation mechanisms. There is evidence that binding molecules function in other intracellular processes [START_REF] Halbleib | Cadherins in development: cell adhesion, sorting and tissue morphogenesis[END_REF][START_REF] Marsden | Integrin-ECM interactions regulate cadherin-dependent cell adhesion and are required for convergent extension in Xenopus[END_REF]. This means that the formation of adhesive bonds can alter the intracellular structure and directly affect cell mobility. In the context of cancer invasion, this was already acknowledged in theoretical models [START_REF] Gerisch | Mathematical modelling of cancer cell invasion of tissue: Local and non-local models and the effect of adhesion[END_REF][START_REF] Ramis-Conde | Modeling the influence of the E-Cadherin-β-Catenin pathway in cancer cell invasion: A multiscale approach[END_REF]. Thus the analogy between liquids and tissues cannot be simply transferred to cell level. Instead a mathematical model is necessary to formulate the concepts on the cell level kinetics explicitly. Then the analysis of this model reveals whether these ideas are consistent with experimental and theoretical findings at cellular and tissue scales.

The question that we address here is how the mesoscopic and macroscopic properties in the process of cell sorting are connected. We conjecture that cell movement is less the more a cell sticks to its neighbors. We translate this assumption into a precise mathematical model. Analyzing this model, we solely rely on stringent mathematical reasoning, since we believe that one does not fully understand the consequences of a model until one can prove rigorous theorems. We are able to predict the emergent behavior at the population level. We qualitatively describe the geometry of cell segregation depending on the intercellular adhesion parameters. Thereby we find a functional relationship between intercellular adhesion and surface tension and highlight the role of cell mobility in the process of sorting. We show that the interaction between the cells and the boundary of a confining vessel has a major impact on the sorting geometry.

Models of Sorting

There are already several mathematical models that address the question of how tissue surface tension is connected with individual cell properties, in particular with the strength of intercellular adhesion. We discuss here two influential models which are to some respects related to our approach, the cellular Potts model of [START_REF] Glazier | Simulation of the differential adhesion driven rearrangement of biological cells[END_REF] and an Ising-like model of Mouchizuki [START_REF] Mouchizuki | A stochastic model for cell sorting and its application[END_REF]. Both models are spatially discrete stochastic equilibrium models which are computationally analyzed with the help of Monte Carlo simulations. It is assumed that typical cell configurations which are observed once the temporal development has stabilized are just those which minimize some (surface) energy functional. This energy functional or 'Hamiltonian' represents the overall interdependence structure between the cells and is therefore dependent on the individual cell characteristics. Although the physical term 'energy' is used here, the constructed Hamiltonians do not represent physical energies but are idealized mathematical concepts used to describe at a macroscopic scale the state of a large number of interacting components. These models, which have their origin in statistical physics, do

A c c e p t e d m a n u s c r i p t

not describe dynamical behavior but characterize the system's state once the temporal evolution has stabilized.

Exploring the energy landscape corresponding to a given Hamiltonian can be hard to accomplish rigorously, therefore one often uses Markov chain Monte Carlo (MCMC) methods. That is, one constructs a Markov chain whose stationary measure is exactly the equilibrium state that one intends to study. The configurations that are visited by the Markov chain after some relaxation period are then typical configurations of the model system. There are various possible methods to construct this auxiliary Markov chain, and one can try to find one which converges exceptionally fast or one which can be readily computed. In any case, it is not clear to what extent this Markov chain resembles the actual or assumed dynamical behavior of the biological system at hand.

In both of the above cited studies, such MCMC methods were applied where the particular Markov transition mechanism allows a meaningful interpretation in terms of the original cell sorting problem. However neither study includes a thorough examination that clarifies whether there is at least a qualitative agreement of essential dynamical characteristics between model system and original cell sorting problem. Therefore it cannot be assumed that the reported kinetic properties emerge from the conjectured interconnection between mesoscopic and macroscopic scale. Instead the model's temporal development can be rather artificial since it is a byproduct of the more-or-less ad-hoc MCMC dynamics.

Glazier, Graner and coworkers [START_REF] Glazier | Simulation of the differential adhesion driven rearrangement of biological cells[END_REF][START_REF] Mombach | Quantitative comparison between differential adhesion models and cell sorting in the presence and absence of fluctuations[END_REF], were the first to address the question of how tissue surface energy might relate to the adhesive bonds between individual cells in a tissue. They developed a model, where cells of certain types 1, ..., K are placed on a two-or three-dimensional square lattice in such a way that each cell covers several lattice sites. The Hamiltonian consists of two terms. The first term weighs the deviation from a cell-type dependent target cell volume and is parameterized by an elasticity parameter λ and the cell-type dependent target cell volumes V 1 , ..., V K . The second term weights the cell surface length and is parameterized by cell-type dependent surface-energy terms (J ij ) i,j=1,...K . Here any particular piece of unit length of a type-i cell surface is given weight J ij , if the cell adjoins with this part of the surface to a cell of type j. Thus the Hamiltonian incorporates variations in cell volume as well as differential adhesion.

Using computer simulations, Glazier and Graner studied the equilibrium measure corresponding to the Hamiltonian, which is the Gibbs measure in mathematical terms. They implemented a mixture of Voter and Metropolis dynamics to sample typical configurations and observed sorted patterns for selected choices of the parameters. This supports their claim that the Hamiltonian can be interpreted to quantify the interfacial and surface free energies of a tissue in the sense of Steinberg's differential adhesion hypothesis.

In addition, they regarded the temporally discrete dynamical process that is given by the chosen Markov transition mechanism as being a dynamical model of the process of sorting. It switches the occupancy of one single lattice site There are some drawbacks of this model. First, spatial cell displacements are modeled implicitly via volume fluctuations and are therefore dependent on the interplay of several parameters that are not directly related to cell mobility nor to concepts on cell-intrinsic events (microscopic scale). Therefore it is not possible to assess to what extent theoretical and empirical ideas on cell migration are reflected correctly. Second, it is not clear whether the implemented Markov transition mechanism results in a correct MCMC algorithm. Since the authors modified the standard Metropolis algorithm (supposedly in order to obtain a biologically more realistic behavior), common results in MCMC theory do no longer apply. Hence it has to be shown explicitly that the stationary measure of the chosen Markov chain is indeed the one which is specified by the given Hamiltonian. Actually, there are indications that the system fixates in the longtime limit in a state where one cell occupies all lattice sites, since these states are absorbing [Brémaud, 1999, §4.6]. This would clearly be an artifact and motivates the question whether the authors really explore the energy landscape of the specified Hamiltonian. If not, the equilibrium behavior is not determined by the Hamiltonian alone which reopens the question about a suitable concept of surface tension and surface energy. And lastly, there is a question of spatial resolution. The decision to model cells in such a way that one cell covers several lattice sites brings about some technical difficulties. For instance, the elasticity parameter cannot be varied much since a single cell would otherwise decompose into several separated components. For the same reason, the standard Metropolis algorithm cannot be applied. This high spatial resolution would be justified if more biological detail was captured. However it is not apparent that the cellular Potts model realistically describes biological properties below cell level.

Mouchizuki and colleagues [START_REF] Mouchizuki | A stochastic model for cell sorting and its application[END_REF] claimed that it is necessary to 'model the stochastic or random movement of the cells explicitly, and model the effect of differential adhesion as a bias of cell exchange rate'. However they rather picked a standard MCMC model, the so-called spin-exchange Ising model, and interpreted it in cell sorting terms.

In detail, cells of two types B, W were placed on a two-dimensional square lattice in such a way that each cell occupies exactly one site and such that there are no empty lattice sites. The strength of homo-and heterotypic intercellular adhesion between two adjoint cells was quantified by parameters λ BB , λ W W and λ BW . A Hamiltonian, called 'total adhesion' was specified to be the weighted sum of all nearest-neighbor bonds. Then the standard Kawasaki dynamics was applied and interpreted as random exchange of cell locations between neighboring cells. This resulting rate for an exchange of spatial cell position is biased towards increasing the total adhesion of the spatial pattern.

Mouchizuki and his coworkers studied the temporal development of their Markov chain model with the help of computer simulations and developed sev- eral order parameters which helped to characterize the sorting process. They studied the equilibrium behavior by applying a heuristic computational routine and proposed an abrupt transition from mixed to sorted equilibrium behavior as some effective adhesion parameter is increased. This effective adhesion parameter depends on λ BB , λ W W and λ BW as well as on some migration parameter m.

Despite the author's claim to present a dynamical model for cell sorting we would like to classify Mouchizuki's Ising-like model as an equilibrium model. Kawasaki dynamics are standard MCMC transition mechanisms for lattice models when one wants to conserve the number of individuals, such as cells or particles, that are located on the lattice. Therefore the corresponding Markov chain can be interpreted to resemble nearest-neighbor cell migration. However we challenge whether this dynamics shows the correct qualitative behavior. For instance, for some situations the rate for a change of spatial positions is higher the stronger the homo-and heterotypic bindings to adjacent cells are. It is not apparent and was not discussed in the cited paper whether there are biological mechanisms that justify this model behavior.

We find that it is necessary to define a model of the dynamical features of cell sorting which describes our conception of individual cell behavior. Then we can explore the implications for the longtime behavior of large cell populations and can check whether our conclusions are consistent with experimental observations.

Often important qualitative characteristics are revealed only with the help of some temporal or spatial limit procedures. Then qualitatively different behavior in certain regions of the phase space becomes apparent and phase transitions can be detected. This cannot be accomplished by computer simulations. Therefore we formulate a model which is analytically tractable.

Interacting particle model

We propose a dynamic individual-based mathematical model for the process of sorting. In particular, we explicitly model the migration of individual cells based solely on spatial competition and the adhesive properties of the individual cells. We focus on pattern formation aspects of sorting which result from active rearrangement of cell configurations and exclude proliferation and apoptosis.

A population of cells of different types is considered, where the cell-types are denoted by numbers 0, 1, ..., n, n ≥ 1.

We assume that all cells regardless of type have similar volume and spheroid shape. It is not appropriate to treat the cells to be points in space, since migration needs to be modeled at the spatial scale of cell size. Instead, we incorporate volume exclusion. We place the cells on a d-dimensional square lattice, d = 2, 3, where each cell occupies exactly one lattice site. Then the lattice spacing corresponds to the cell diameter. Each arrangement of cells on the lattice is a configuration. Note that the discreteness of space makes the model simpler to handle and analyze.

A c c e p t e d m a n u s c r i p t

We aim to understand the dynamics of cell populations that consist of n different cell types and whose spatial arrangement can have medium in between. To this end, we interpret the cell type 0 to be an unoccupied lattice point and types 1, ..., n to be biological cells of certain type.

However there is another interpretation in case n = 1, that is W = {0, 1}. In this case, each arrangement of 0's and 1's on the lattice models a cell population of two different cell types without medium in between. We call this a fully occupied system.

We denote the lattice by S ⊂ Z d , d = 2, 3. Here we concentrate on the two cases that S = Z d or that S is bounded, that is |S| < ∞. We define T := {T ⊂ S : |T | < ∞}, where |T | is the cardinal number of a set T ⊂ S. Let be W = {0, 1, ..., n} the set of all possible cell types. Then the configuration space is given by

X = W S = {η : S → W }.
A configuration is an element of the configuration space,

η = {η(x)} x∈S ∈ X.
With only one cell of type j, cell migration is modeled by a simple nearestneighbor random walk on the lattice. The corresponding transition matrix is given by p j := α j p = (α j p(x, y)) x,y∈S with p(x, y) = 1, |x -y| = 1 0, otherwise.

(3.1)

The parameter α j ≥ 0 characterizes the migration speed of type-j cells, j ∈ W , and the expression |x -y| denotes the Manhattan distance between two lattice sites x = (x 1 , ..., x d ), y = (y 1 , ..., y d ) ∈ S, that is

|x -y| := |x 1 -y 1 | + ... + |x d -y d |.
We assume that α j > 0 for j = 0. If there is more than one cell on the lattice but no adhesive interplay, we define that the cells interact by virtue of the exclusion principle. A randomly chosen cell of type j attempts to jump on a given neighboring site with rate α j . If the target site is already occupied by a cell of type i, then both cells interchange their positions with rate α i + α j . Effectively, the configuration changes only for i = j. Given a configuration η ∈ X, the transition

η → η xy with η xy (z) := ⎧ ⎪ ⎨ ⎪ ⎩ η(y), z = x, η(x), z = y, η(z), z = x, z = y,
describes the change in the configuration as the cell at x jumps to an empty site y, or vice versa, or if both locations x, y are occupied the cells at x and y interchange their positions. In the case that both sites x, y are empty or occupied by cells of the same type, the configuration is not altered by the above transition. We specify the rates for a transition η → η xy to be

c 0 (x, y, η) = (α η(x) + α η(y) )p(x, y) if η(x) = η(y) 0 otherwise,
where x, y ∈ S, η ∈ X. For W = {0, 1} this model is a so-called symmetric nearest-neighbor exclusion process [START_REF] Liggett | Interacting Particle Systems[END_REF].

Adhesive interactions between cells bias the cell migration. If a cell strongly sticks to neighboring cells then the cell mobility is lessened. Thus the rate for the cell's transition to one of the neighboring sites is decreased. We assume that adhesion works only between cells in contact. Therefore we only consider bonds between cells on neighboring sites, that are sites at Manhattan distance one. The strength of the bond between two cells depends on the cell types. Thus the bond between two cells of type i and j, respectively, is weighted by a value β ij ∈ R. Naturally, we assume the matrix

β = (β ij ) i,j∈W
of adhesion parameters to be symmetric. For x, y ∈ S, η ∈ X, we model the rate for a transition η → η xy to be c(x, y, η) = c 0 (x, y, η) exp

⎧ ⎨ ⎩ - z:|z-x|=1 β η(x)η(z) - z:|z-y|=1 β η(y)η(z) ⎫ ⎬ ⎭ = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (α η(x) + α η(y) )p(x, y)× × exp -z:|z-x|=1 β η(x)η(z) -z:|z-y|=1 β η(y)η(z) , if η(x) = η(y), 0, otherwise. (3.2)
The interpretation is as follows. Suppose we have a cell of type i at position x and a cell of type j at y, then both cells exchange their positions with a rate that is proportional to exp

⎧ ⎨ ⎩ - z:|z-x|=1 β iη(z) - z:|z-y|=1 β jη(z) ⎫ ⎬ ⎭ .
This rate is smaller the more the cells adhere to neighboring cells. Positional exchanges of other than nearest-neighbor cells are suppressed. See Figure 1 for an example. Note that for W = {0, 1} the model falls in the class of so-called stochastic lattice gases with speed change [START_REF] Spohn | Large Scale Dynamics of Interacting Particles[END_REF].

The transition rates embody the local behavior of cell mobility. We do not assume that there are any chemical markers which guide cell behavior, therefore the cells have no ability to gain knowledge of their more distant or even global Type-1 cells (dark gray) and type-2 cells (light gray) are arranged on a twodimensional square lattice. Left: The dark cell located within the shadowed area can move freely since it has no nearest neighbors and is not subject to intercellular adhesion therefore. A transition to the empty site within the shadowed area happens with rate α 1 . Right: After the transition has taken place, the dark cell adheres to two neighboring type-2 cells. The rate for the reverse transition is decreased to α 1 exp{-2β 12 }.

A c c e p t e d m a n u s c r i p t

environments. The cell dynamics depends solely on the cell configuration in the direct neighborhood of the regarded cells.

We assume that cell jumps to neighboring lattice sites follow a stochastic rule. This reflects the observation that, besides adhesive bindings and spatial competition, a variety of minor influences determines the migration behavior of biological cells. Those influences can be the inner cell structure, which is not modeled explicitly here, inhomogeneities of the environment or fluctuations in cell shape and volume. All of them together give rise to the stochastic component of the model. We suppose further that the future evolution of the cell population depends solely on the present spatial arrangement of the cells, that is the present cell configuration. We therefore formulate a Markov process model.

We treat time as a continuous quantity. A discretization of time, as very often used for computer simulations, would not help us to analyze the model. From our analytical viewpoint, continuous-time models are easier to handle and more natural. See [START_REF] Klauss | Modelling and simulation by stochastic interacting particle systems[END_REF] for the derivation of a according time-discrete model.

The family c = (c(x, y, .)) x,y∈S of transition rates which are specified by (3.2) describes the dynamics at the cellular scale. It determines a Markov process on X which shall be called the cell sorting model. For bounded volumes, that is for |S| < ∞, this Markov process is a finite Markov chain, which is well-defined as soon as the effect of the boundary is described. This will be done in the following. The specifics of how to define the cell sorting model for S = Z d are presented at the end of this section.

If the lattice is finite, that is if |S| < ∞, the cell sorting model is uniquely determined as soon as the initial state and the boundary conditions are fixed. We consider the cases of periodic and fixed boundary conditions. In both cases, we alter the expression in (3.2) for those cells that are located next to the boundary of the considered volume in order to obtain a meaningful behavior near the boundary.

If

S = Λ N := [0, N -1] d ∩ Z d , N ∈ N, we can define periodic boundary conditions. For N ∈ N, let be N Z d := z ∈ Z d : z = (Nx 1 , ..., N x d ) for some x 1 , ..., x d ∈ Z and |x| Λ N := min{|y| : y ∈ x + N Z d }, x ∈ Λ N ,
where z+B

:= {y+z : y ∈ B}, z ∈ R d , B ⊂ R d .
A N p f (η) := x,y∈Λ N c N p (x, y, η)(f (η xy ) -f (η)), f ∈ C(X N ), η ∈ X N . (3.3)
The generator A N p is a bounded operator, since X N is a finite set. In the context of Markov chains one often works with the Q-matrix a N p instead of the generator, which is is given by the relation

a N p := a N p (η, ξ) η,ξ∈X N = A N p 1 {ξ} (η) η,ξ∈X N . (3.4)
Note that the correspondence between a N p and A N p is one-to-one. It is wellknown from standard Markov chain theory (see, for example, [START_REF] Brémaud | Markov chains. Gibbs fields, Monte Carlo simulation and queues[END_REF]) that for each probability measure ν 0 on X N the operator A N p uniquely determines a continuous-time Markov chain with initial distribution ν 0 and state space X N which has right continuous paths with left limits (cadlag paths).

Definition 3.1. Fix N ∈ N and choose a probability measure ν 0 on X N . Let c N p = (c N p (x, y, •))
x,y∈Λ N be the family of transition rates which is defined above. The corresponding continuous-time Markov chain on X N which has cadlag paths and initial distribution ν 0 is called cell sorting model in Λ N with periodic boundary conditions and with start in ν 0 .

If we want to define a cell sorting model with fixed boundary condition, we set S = Z d , fix a volume V ⊂ S with |V | < ∞ and choose a so-called external condition ζ ∈ X. During the temporal development, the external condition ζ shall remain fixed but it will influence the rates for cell migration next to the boundary. To accomplish this, we define,

for ζ ∈ X, v ∈ X V := W V , a configuration τ V (ζ, v) ∈ X by τ V (ζ, v)(z) := ζ(z), z ∈ S\V, v(z), z ∈ V. (3.5) We put c V ζ (x, y, η) := c(x, y, τ V (ζ, η)), x,y ∈ V, η ∈ X V , (3.6)
where c(x, y, .) is given by (3.2), and obtain a family

c V ζ = (c V ζ (x, y, •)) x,y∈V of transition rates on X V . As above, the family c V ζ determines a Markov generator A V ζ : C(X V ) → C(X V ) by A V ζ f (η) := x,y∈V c V ζ (x, y, η)(f (η xy ) -f (η)), f ∈ C(X V ), η ∈ X V . (3.7)
Here C(X V ) is the set of all real functions on X V . The corresponding Q-matrix a V ζ is given by

a V ζ = a V ζ (η, ξ) η,ξ∈X V = A V ζ 1 {ξ} (η) η,ξ∈X V . (3.8)

A c c e p t e d m a n u s c r i p t

For each probability measure ν 0 on X V , the operator A V ζ uniquely determines a continuous-time Markov chain with initial distribution ν 0 and state space X V which has cadlag paths. Definition 3.2. Fix V ⊂ S and choose a probability measure ν 0 on X

V . Let c V ζ = (c V ζ (x, y, •))
x,y∈V be the family of transition rates which is defined in (3.6). The corresponding continuous-time Markov chain on X V which has cadlag paths and initial distribution ν 0 is called cell sorting model in V with fixed boundary condition ζ and with start in ν 0 .

The dynamics of a cell sorting model in bounded volume V ⊂ S with fixed boundary condition ζ ∈ X is easy to illustrate. Given a present configuration η ∈ X V , imagine a clock at each (unordered) pair of neighboring sites {x, y}, x, y ∈ V, |x -y| = 1, which is set to an exponential time with parameter c V ζ (x, y, η). The system waits until the first clock rings, say it is the clock at the pair {x 0 , y 0 }. Then the transition η → η x0y0 is performed, that is the occupancy at x 0 and y 0 is exchanged, and all clocks are set anew to exponential times, now with parameters c

V ζ (•, •, η x0y0 ).
Actually, any cell sorting model in bounded volume develops in this way. One only has to substitute the respective transition rates for c V ζ .

In addition to the cell sorting model in bounded volumes it is necessary to have a well defined model for S = Z d as well. This is because important qualitative characteristics of the model cannot be revealed until one applies suitable spatiotemporal limit procedures. Therefore it is essential to have a well-defined stochastic process even for S = Z d . However, for unbounded S, the state space X is innumerable. Therefore standard results of Markov chain theory no longer apply and the question of existence and uniqueness of a Markov process corresponding to the transition rates in (3.2) becomes non-trivial. In the following we show that the cell sorting model fits into the concept of interacting particle systems (IPS) as defined in [START_REF] Liggett | Interacting Particle Systems[END_REF]. We apply general results on the existence and uniqueness of interacting particle systems which ensure that the cell sorting model is well-defined. To this end, we collect essential properties of the transition rates. Proposition 3.3. Let the family of transition rates (c(x, y, •)) x,y∈S be given by (3.2). Then the following properties are fulfilled.

(A) translation invariance:

c(x, y, η) = c(x + z, y + z, θ z η), x,y,z ∈ S, η ∈ X, where θ z : X → X : (θ z η)(u) := η(u -z), u ∈ S, for any z ∈ S. (B) finite range jumps: c(x, y, •) = 0 if |x -y| > 1. (C) local rates: c(x, y, η) = c(x, y, ζ) for any x, y ∈ S and η, ζ ∈ X with η(z) = ζ(z) for z ∈ {z ∈ S : |z -x| ≤ 1 or |z -y| ≤ 1}.

A c c e p t e d m a n u s c r i p t

(D) positive rates:

inf x,y∈S |x-y|≤1 inf η∈X, η(x) =η(y) c(x, y, η) > 0.
(E) irreducible rates: For any two different sites x, y ∈ S there exists a finite sequence x = x 0 , x 1 , ..., x n = y ∈ S such that

x i = x i-1 and min{c(x i-1 , x i , η) : η ∈ X, η(x i-1 ) = η(x i )} > 0 for all 1 ≤ i ≤ n. (F) uniformly bounded rates: sup x∈S y =x sup η∈X c(x, y, η) < ∞.
Proof. Properties (A)-(C) are found directly from the definition of the respective rates. We show that properties (A)-(C) imply (F): Let o = (0, ..., 0) ∈ Z d . By (B) and (C) we have

sup η∈X c(o, y, η) = ⎧ ⎪ ⎨ ⎪ ⎩ 0, if |y| = 1 sup η∈X η(z)=0,|z|>2 c(x, y, η) < ∞, if |y| = 1.
β ij > 0 > 0 for any x, y ∈ S with |x -y| ≤ 1 and any η ∈ X with η(x) = η(y). Finally, (D) gives that min{c(z 1 , z 2 , η) : η ∈ X, η(z 1 ) = η(z 2 )} > 0, if |z 1 -z 2 | = 1.
Since there is always a path along the edges of the lattice Z d which connects to given points x, y ∈ S, this implies (E).

By Liggett [1985, Ch. 1], each family of transition rates c = c(x, y, •) x,y∈S which satisfies (A)-(C) above uniquely defines a Markov process on X, which is a called a Feller process. We introduce some formalism now in order to make this statement sound.

Suppose that the set W is equipped with the discrete metric d. Then (W, d) is a compact metric space and the σ-field W = B(W ) of all Borel sets coincides with the power set of W . The product space X = W S is equipped with the product topology which is metrizable. Note that X is compact w.r.t. the product topology. The product σ-field F = W S coincides with the Borel σ-field B(X).

The space of all continuous real functions on X equipped with the supremum norm is denoted by C(X). This set is complete and separable. A real function

A c c e p t e d m a n u s c r i p t

f on X is called a tame function if there is a finite set T ⊂ S, |T | < ∞ such that f (η) = f (ζ) for any η, ζ ∈ X with η(x) = ζ(x) for x ∈ S\T .
Thus tame functions depend only on finitely many coordinates. It is easily verified that the (countable) set T (X) of all tame functions is a dense subset of C(X) [Georgii, 1988, (2.21)].

We are given an operator A which acts on T (X) by

(Af )(η) = x,y∈S c(x, y, η) (f (η xy ) -f (η)) , f ∈ T (X), η ∈ X. (3.9)
This operator is well-defined, that means, for f ∈ T (X), the above infinite sum converges absolutely and defines a continuous function on X. This follows from Liggett [1985, Prop.I.3.2], since the family of transition rates (c(x, y, •)) x,y∈S satisfies the conditions of Proposition 3.3 (A)-(C). The closure of A is the generator of a Markov semigroup (T t ) t≥0 on C(X) [Liggett, 1985, Th.I.3.9]. Given a probability measure ν 0 on (X, F), the operator A uniquely determines a Markov process (η t ) t≥0 with initial distribution ν 0 and state space (X, F), which has cadlag paths [Liggett, 1985, Thm.I.1.5]. Note that ν t , the distribution of η t , is given by the relation

ν t (f ) = ν 0 (T t f ), f ∈ C(X), t ≥ 0. Definition 3.4. Let c = c(x, y, •)
x,y∈S be the family of rate functions which is given by (3.2) and let A be the corresponding Markov pre-generator that is defined in (3.9). The associated Markov process which has cadlag paths and initial distribution ν 0 ∈ P(X, F) is called cell sorting model in Z d with start in ν 0 .

We end up with a stochastic, continuous-time, individual-based model of interacting cells on a regular square-lattice, which can be used to study the emergent dynamic properties on tissue scale. Based on the commonly accepted view that cell mobility and differential intercellular adhesion are essential for tissue segregation, we present a model where the two factors interact dynamically. We argue that cell adhesion dampens cell mobility which means that the cells are the less mobile the more they stick to each other. This assumption is translated into a precise mathematical model which belongs to the class of stochastic interacting particle systems. Any further details at the cellular level which might have an effect on cell mobility or intercellular adhesion such as the control of cell shape or cell size are considered to be of minor importance and are integrated into the stochastic component. We have two model parameters, the mobility parameter α, that regulates the rate of cell displacement, and the adhesion parameter β = (β ij ) i,j∈W , which describes the strength of intercellular adhesion. This allows to study the effect of each factor separately.

Our IPS-model for cell sorting is defined via the specification of (local) transition rates. This reflects our understanding that it is the concept of the dynamics at the cellular scale which will be reproduced by the model. Mathematically,

A c c e p t e d m a n u s c r i p t

our model for cell sorting is a Markov process that describes the temporal development of the spatial cell arrangement within a discretized fixed volume. In the case where the considered volume is bounded, the model simplifies to a finite Markov chain and can be described by standard mathematical methods.

In order to reveal essential qualitative properties when analyzing such systems, it is a fundamental necessity to go beyond finite systems and to study the model in its infinite-volume limit. However in the case of unbounded volumes, the state space becomes innumerable and the corresponding stochastic model would be no longer a Markov chain. Therefore, we have chosen the technically more involved approach to develop an IPS-model, which can be considered as an infinite-volume analog of the finite system and includes the finite model as a special case.

Mathematical Analysis

In this section, the key arguments of the analysis are presented. Predominantly technical and intermediate findings as well as detailed proofs are shifted to the Appendix. An elucidating summary of the analytic findings and a discussion of the consequences for the model behavior can be found in Section 5.

We concentrate on an analysis of the longtime behavior of our model together with a characterization of the steady states. Since we have a stochastic model, the steady states are probability measures which are invariant under the dynamics. We show that the invariant measures are reversible with respect to the cell sorting dynamics, that is we prove that the probability flux between any two states is balanced. The reversible measures are characterized further to be Gibbs measures. This property is exploited to answer the question of which typical spatial arrangements the model system exhibits in the longtime limit.

Gibbs measures are equilibrium models which originate from statistical mechanics (see, for instance [START_REF] Georgii | Gibbs measures and phase transitions[END_REF] for their mathematical foundation). They are probability measures on the configuration space X which are defined via potentials and their corresponding Hamiltonians. The potential is a postulate about the system's interdependence structure once the temporal evolution has stabilized. It is the starting point for equilibrium modeling. Hamiltonians and Gibbs measures are constructed from the potential in a standard way. By analyzing them the consequences of the postulated interdependence structure in equilibrium are explored.

Our approach is contrary. Instead of conjecturing an interdependence structure for the system's steady states, we prove that the model dynamics implies a special structure of the invariant measures. We show that the invariant measures are Gibbs measures with respect to a certain potential which is specified explicitly. Thus the interdependence structure in equilibrium emerges as a consequence of the assumptions on the system dynamics.

We start with an investigation of the infinite model in paragraphs 4.1 and 4.2, that is we study the cell sorting model in S = Z d . The model behavior in large but bounded volumes is analyzed in paragraph 4.3.

A c c e p t e d m a n u s c r i p t

Reversible measures are Gibbs measures

Suppose that S = Z d and let be given a family c of transition rates by (3.2). The set of all probability measures on (X, F) is denoted by P(X, F). We recall that we have a Markov generator A and a Markov semigroup (T t ) t≥0 on C(X) which are uniquely associated to the family c.

Definition 4.1. A measure μ ∈ P(X, F) is invariant w.r.t. (T t ) t≥0 if and only if μ(f ) = μ(T t f ), f ∈ C(X), t ≥ 0.
The set of all invariant measures w.r.t. (T t ) t≥0 is denoted by

I = I (c).
There is a subset of I that is often easier to characterize. This is the set of reversible measures on (X, F). Reversible measures are characterized by the property that the probability flux between any two states is balanced.

Definition 4.2. A measure μ ∈ P(X, F) is reversible w.r.t. (T t ) t≥0 if and only if μ(fT t g) = μ(gT t f ), f,g ∈ C(X), t ≥ 0.
The set of all reversible measures w.r.t. (T t ) t≥0 is denoted by R = R(c).

Note that R ⊂ I and I = ∅, but, in general, R = ∅ is possible [Liggett, 1985, I.1.8(f) and I.5.2].

The concept of canonical Gibbs measures that was developed within the context of statistical physics is well-suited to characterize the reversible measures in our model. Gibbs measures are defined with respect to some potential which determines the behavior inside bounded volumes conditioned on the outside configuration (see [START_REF] Georgii | Gibbs measures and phase transitions[END_REF] as an excellent general reference).

In detail, define, for T ⊂ S, the local configuration space X T := W T , the corresponding projection π T : X → X T : π T (η) = {η(x)} x∈T =: η T and the σ-algebra

F T = π -1 T (W T ). Let us agree to write π -1 T (v) for the set π -1 T ({v}) = {η ∈ X : π T (η) = v}, v ∈ X T . We recall the definition T = {T ⊂ S : |T | < ∞}.
Let us recollect the formula (3.5), where transformations τ T : X × X T → X are defined for T ∈ T . The configuration τ T (η, v), which is given by

τ T (η, v)(z) = η(z), z ∈ S\T v(z), z ∈ T, arises from η ∈ X by substituting v ∈ X T for the sub-configuration η T = (η(x)) x∈T . Definition 4.3. A family Φ = (Φ T ) T ∈T of F T -measurable real functions Φ T :
X → R is a finite range potential, if there is a real number r > 0 such that Φ T = 0 for any T ∈ T with max

x,y∈T |x -y| > r.

A c c e p t e d m a n u s c r i p t

Each finite range potential Φ can be assigned to a family (H Φ T ) T ∈T of Hamiltonians by

H Φ T (η) = V ∈T V ∩T =0 Φ V (η), T ∈ T , η ∈ X, (4.1)
and a family (h Φ T ) T ∈T of Boltzmann factors by the relation

h Φ T (η) = exp -H Φ T (η) , T ∈ T , η ∈ X. (4.2)
For w ∈ W, T ∈ T , denote

N T : X × W → N 0 : N T (η, w) = x∈T 1 {ζ:ζ(x)=w} (η). (4.3)
That is, N T (η, w) gives the number of type-w cells within the volume T . Put N T = (N T (w, •)) w∈W . The σ-algebra of events which are invariant under permutation of the sites in T is denoted by

E T = σ(π S\T , N T ) ⊂ F.
The term

Z T (ζ) = η∈O T (ζ) h T (η), ζ ∈ X,
μ(π -1 T (v)|E T )(ζ) = Z -1 T (ζ) h Φ T (τ T (ζ, v))1 N -1 T (N T (ζ)) (τ T (ζ, v)), μ -a.s., where T ∈ T , v ∈ X T , ζ ∈ X.
The set of all canonical Gibbs measures w.r.t. Φ is denoted by G c (Φ).

After having defined canonical Gibbs measures, we are able to state that the reversible measures of the cell sorting model in infinite volume S = Z d agree with the canonical Gibbs measures w.r.t the potential Φ.

Theorem 4.5. Let the family c = (c(x, y, •)) x,y∈S of transition rates be given by (3.2). Then there exists a finite range potential Φ = (Φ T ) T ∈T such that

R(c) = G c (Φ) = ∅.

A c c e p t e d m a n u s c r i p t

We prove this theorem by showing that the transition rates c satisfy a socalled detailed balance condition w.r.t. to the finite range potential Φ = (Φ T ) T ∈T given by

Φ T : X → R : η → -β η(x)η(y) , if T = {x, y} with |x -y| = 1 0, otherwise. (4.4)
See the Appendix, §AppendixA.1, for details. Note that the corresponding Hamiltonian (H Φ T ) T ∈T as defined in (4.1) satisfies

H Φ T (η) = - {x,y}:|x-y|=1 {x,y}∩T =∅ β η(x)η(y) , η ∈ X, T ∈ T . (4.5)
The potential Φ is not uniquely determined by the transition rates of the model. In case W = {0, 1}, we can characterize the set of reversible measures with the aid of an especially simple potential. Let be

β * := β 00 -2β 01 + β 11 , ( 4.6)
and define Φ = ( ΦT ) T ∈T by

ΦT (η) = -β * η(x)η(y), if T = {x, y} with |x -y| = 1 0, otherwise. (4.7)
It is clear that Φ is a finite range potential in the sense of Definition 4.3. The potential Φ is known in the mathematical literature under the name Ising potential.

As is shown in Lemma AppendixA.4 in the Appendix, §AppendixA.2, the family c of transition rates satisfies the detailed balance condition w.r.t. the finite range potential Φ = ( ΦT ) T ∈T given by (4.7). As an immediate consequence we find: Proposition 4.6. Let c = (c(x, y, •)) x,y∈S be given by (3.2) and Φ = ( ΦT ) T ∈T given by (4.7). Then

R(c) = G c ( Φ) = ∅.
Thus the findings of this paragraph are the following. The reversible measures of the cell sorting model on S = Z d are canonical Gibbs measures w.r.t. the potential Φ. If the fully occupied system is considered, we find that the potential Φ is equivalent to the simpler and well-studied Ising potential. This implies that the reversible measures of the fully occupied cell sorting model on S = Z d coincide with the canonical Gibbs measures with respect to the Ising potential.

I (c) = G c ( Φ) = ∅.
Proof. We apply Sakagawa's Theorem 2.1 [START_REF] Sakagawa | Equivalence between canonical Gibbs measures and stationary measures for stochastic lattice-gas models[END_REF]. This theorem is valid if the exchange rates satisfy 3.3 (B)-(F) as well as the detailed balance condition (A.1) w.r.t. a finite range potential Φ = ( ΦT ) T ∈T that fulfills

ΦT (η) = φ(T ) x∈T η(x) (4.8)
for some φ :

T → R. If we choose φ(T ) = β * , T = {x, y} with |x -y| = 1, 0, otherwise, (4.9)
then we find that the assertion of the theorem is implied by Proposition 3.3 and Lemma AppendixA.4.

There is no such result in dimension d ≥ 3 or for n ≥ 2.

Large finite systems

The cell sorting model on Z d is an idealization which is developed as a framework for spatio-temporal limit procedures. Actually, one is interested in the behavior of systems in large but finite volumes. The object of this paragraph is the longtime behavior of the cell sorting model in bounded volumes. We study the ergodic properties of those systems, that is we describe the set of invariant measures and their domains of attraction. The question of which spatial arrangement the model system shows in the long-time limit is addressed in Paragraph 4.4.

For reasons of clarity, we state the results on the ergodic behavior for the cell sorting model with fixed boundary conditions as defined in Definition 3.2. Note that all arguments apply to the cell sorting model with periodic boundary conditions as well.

We fix a bounded volume V ∈ T and an external condition ζ ∈ X. Suppose that the family c V ζ of transition rates is given by (3.6). We recall that c V ζ determines a Markov chain on X V = W V with Q-matrix a V ζ and generator A V ζ . Let us denote the set of all probability measures on X V by P(X V ).

Definition 4.8. (i) A measure μ ∈ P(X V ) is invariant w.r.t. A V ζ if and only if μ A V ζ f = 0, f ∈ C(X V ).
The set of all invariant measures w.r.t.

A V ζ is denoted by I c V ζ .
A c c e p t e d m a n u s c r i p t

(ii) A measure μ ∈ P(X V ) is reversible w.r.t. A V ζ if and only if μ(fA V ζ g) = μ(gA V ζ f ), f,g ∈ C(X V ).
The set of all reversible measures w.r.t.

A V ζ is denoted by R c V ζ .
We observe that the cell numbers within V are conserved under the dynamics induced by A V ζ . Indeed, we recall that N V (.) is the vector of cell numbers defined in (4.3) and we choose κ ∈ N V , where

N V := κ = (k 0 , ..., k n ) ∈ N W 0 : n i=0 k i = |V | .
Then the set

X κ V := {ζ ∈ X V : N V (ζ) = κ} is invariant (closed ) under the action of A V ζ .
It follows from proposition 3.3(D), that any two configurations with the same cell numbers per type communicate. The latter means that, for any η, ζ ∈ X κ V , there exist an n ∈ N such that (a V ζ ) n (η, ζ) > 0, that is ζ can be obtained from η by performing a finite number of nearest-neighbor exchange steps which each have positive rate. That way, the state space X V is partitioned into closed communicating classes X κ V , κ ∈ N V . This is the basis for the following assertions, which are implied by standard Markov-chain theory (see, for example, [START_REF] Brémaud | Markov chains. Gibbs fields, Monte Carlo simulation and queues[END_REF]). Proposition 4.9. (i) Assume that κ ∈ N V . There exists a unique probability measure μ κ ζ on X V which is invariant with respect to a V ζ and satisfies

μ κ ζ (X κ V ) = 1. (ii)
The set of invariant probability measures with respect to A V ζ is given by

I V := μ ∈ P(X V ) : μ = κ∈N V γ κ μ κ ζ , γ κ ≥ 0, κ∈N V γ κ = 1 .
After having described the invariant measures of the cell sorting model in bounded volume, we turn to the question of which probability distribution on X V is observed in the longtime limit limit.

Theorem 4.10. Let ν 0 denote the initial distribution of some Markov process generated by A V ζ . Then the distribution at time t, referred to as ν t , t ≥ 0, satisfies lim

t→∞ ν t (η) = κ∈N V ν 0 (X κ V )μ κ ζ (η), η ∈ X V .
Now we characterize further the invariant measures of the finite adhesion model. We will show that each invariant measure is actually reversible and coincides with a bounded-volume canonical Gibbs state with respect to the potential Φ defined in (4.4). To this end, we first define canonical Gibbs measures with As reversible measures are invariant, we can employ uniqueness arguments to show that the invariant measures coincide with the canonical Gibbs measures.

Definition 4.11. Suppose that Φ is a finite range potential and let H Φ V and h Φ V be given by (4.1) and (4.2). For V ∈ T and ζ ∈ X, the probability measure μ

V,ζ ∈ P(X V ) is a Gibbs measure w.r.t. Φ in V with boundary condition ζ, if μ V,ζ (η) = Z -1 V,ζ exp -H Φ V (τ V (ζ, η)) = Z -1 V,ζ h Φ V (τ V (ζ, η)), η ∈ X V , where Z V,ζ := η∈X V h Φ V (τ V (ζ, η))
is the normalizing factor. The family (μ V,ζ ) V ∈T is a grand-canonical Gibbs ensemble w.r.t. the potential Φ and the boundary condition ζ.

If the cell numbers per type are prescribed, we arrive at the concept of canonical Gibbs measures in bounded volume. Definition 4.12. Suppose that κ ∈ N V and that Φ is a finite range potential.

Let H Φ

V and h Φ V be given by (4.1) and (4.2). The probability measure μ V,ζ,κ ∈ P(X V ) is a canonical Gibbs measure with potential Φ, boundary condition ζ and cell numbers κ, if

μ V,ζ,κ (η) = Z -1 V,ζ,κ exp -H Φ V (τ V (ζ, η)) 1 N -1 V (κ) (η) = Z -1 V,ζ,κ h Φ V (τ V (ζ, η))1 N -1 V (κ) (η), η ∈ X V , where Z V,ζ,κ := η∈X V h Φ V (τ V (ζ, η))1 N -1 V (κ) (η)
is the normalizing factor. Given a net (κ(

V ) = (k 1 , ..., k n )(V )) V ∈T with κ(V ) ∈ N V and lim V ∈T |V | -1 k j (V ) = p j , j ∈ W, the family (μ V,ζ,κ(V ) )
V ∈T is a canonical Gibbs ensemble w.r.t. Φ, boundary condition ζ and cell densities (p j ) j∈W .

After having defined the canonical Gibbs states with fixed boundary condition, we state that they are reversible measures of the cell sorting model with fixed boundary condition. The proof of this proposition can be found in the Appendix, §AppendixA.3. Proposition 4.13. Suppose that a potential Φ is given by (4.4) and let be κ ∈ N V . Then the canonical Gibbs measure μ V,ζ,κ corresponding to Φ is reversible with respect to A V ζ . Since the canonical Gibbs measure μ V,ζ,κ is reversible, it is invariant w.r.t the dynamics A V ζ . In addition, it holds that μ V,ζ,κ (X κ V ) = 1. By Proposition 4.9(i), there is exactly one invariant measure

μ κ ζ that satisfies μ κ ζ (X κ V ) = 1. Thus μ κ ζ = μ V,

A c c e p t e d m a n u s c r i p t

Corollary 4.14. Let ν 0 ∈ P(X V ) denote the initial distribution of some Markov process generated by A V ζ . Suppose that ν 0 (X κ V ) = 1 for some κ ∈ N V . Then the distribution at time t, referred to as ν t , t ≥ 0, satisfies

lim t→∞ ν t (η) = μ V,ζ,κ (η), η ∈ X V .

Phase segregation

We consider the question of which spatial arrangement the model system shows in the long-time limit. It will be shown that the equilibrium measures of large finite systems exhibit a typical geometry, which is disclosed by a renormalization procedure.

Suppose that W = {0, 1} and d = 2. There is a critical parameter value β * c > 0, such that for β * > β * c the Gibbs ensemble has more than one limit point as the volume is enlarged to Z d . This non-uniqueness is reflected in a phase segregation phenomenon for the canonical Gibbs ensemble. That is, for large finite volume V and fixed β * > β * c the canonical Gibbs measure in V concentrates on those configurations that are characterized by high type-1-cell density p * + (β * ) inside some connected component and small type-1-cell density p * -(β * ) outside this component. The values β * c and p * ± can be explicitly calculated, see (4.10) and (4.11) below. The shape of the type-1 cluster can be detailed further by a variational problem, the Wulff construction.

We follow Dobrushin, Kotecky and [START_REF] Dobrushin | Wulff construction. A global shape from local interaction[END_REF] for a sketch of the main ideas. For reasons of clarity, we state the results for canonical Gibbs ensembles with periodic boundary conditions, as defined in detail in the appendix. This corresponds to the situation where boundary effects are neglected. Actually, phase segregation occurs for models with fixed boundary condition as well. The boundary conditions have a major impact on the shape of the segregated regions. This effect is discussed in Paragraph 4.5.

We recall that the cell sorting model with periodic boundary conditions is described in Definition 3.1. By considerations analogous to those in Paragraph 4.3, we find that the longtime behavior of the cell sorting model in Λ N = [0, N -1] ∩ Z d with periodic boundary condition and type-1-cell number R is described by the measure μ β * ,N,R , which is a canonical Gibbs measure on Λ N with respect to the potential Φ and with periodic boundary conditions. See §AppendixA.4 for the exact definition of μ β * ,N,R . We will derive an asymptotic result for the measure μ β * ,N,R on Λ N as the length-size N of the cube Λ N tends to infinity. To accomplish this, we embed Λ N into the continuous torus T = R d /Z d of diameter 1. The details of this embedding can be found in §AppendixA.5. The main idea is that a configuration on Λ N is identified with the empirical measure of the suitably scaled type-1-cell positions within this configuration. Carrying forward this idea to the measures on X N , a measure measure μ N on X N = W Λ N is transformed into a measure μN on the space M( T ) of bounded Borel measures on T . Further, we define The interpretation is as follows. The Borel measure δ γ+x on the torus T has the density function p * + 1 U1+x (y) + p * -1 U0+x (y), y ∈ T . This function attains the values p * + on the set U 1 + x and p * -on the set U 0 + x. The two regions are separated by the curve γ + x, which is a shift of γ. Thus, if the ensemble (μ β * ,N,R N ) N ∈N exhibits asymptotic phase segregation along some curve γ, then the canonical Gibbs measures μ β * ,N,R N , when rescaled to unit volume, are asymptotically concentrated on such configurations where a high-density phase with characteristic density p * + is separated from a low-density phase with characteristic density p * -by the curve γ or one of its shifts γ + x. Note that the absolute position of the boundary γ cannot be fixed since the canonical Gibbs measures with periodic boundary condition are translation invariant. Therefore it is evident that one finds in (4.14) an equiprobable distribution of all shifts γ + x of the curve γ.

β * c := 2 ln(1 + √ 2) = 2 sinh -1 (
We can now state the result on phase segregation in the cell sorting model with periodic boundary condition.

A c c e p t e d m a n u s c r i p t

Theorem 4.16. Suppose that d = 2. For any

β * > β * c , p ∈ (p * -(β * ), p * + (β * ))
, there exists some curve γ β * ,p such that for any sequence (R N ) N ∈N satisfying (4.13) the ensemble (μ β * ,N,R N ) N ∈N exhibits asymptotic phase segregation along the curve γ β * ,p . Remark 4.17. (1) The curve γ β * ,p is yielded by the Wulff construction. It is characterized by the property that a functional called total surface tension is minimized among all closed self-avoiding rectifiable curves which enclose a given volume [START_REF] Dobrushin | Wulff construction. A global shape from local interaction[END_REF]. In detail, there is a functional W, which assigns to each subset T ⊂ T with rectifiable boundary ∂T the total surface tension W(∂T ). The total surface tension W can be calculated directly from the Hamiltonian (4.5) though not in a simple way. Compare Dobrushin et al. [1992, §1.5] for the explicit formula. Among all subsets T ⊂ T of a given volume v and with rectifiable boundary, the Wulff droplet of size v is that set W β * ,v ⊂ T which minimizes the total surface tension, that is

W(∂W β * ,v ) = min T ⊂ T ,volume(T )=v, ∂T is rectifiable W(∂T ).
The density of type-1 cells in the Wulff droplet W β * ,v is p * + , while the density of type-1 cells outside of W β * ,v is p * -. Since the overall density of type-1 cells is p, we find that the volume of the Wulff droplet is

v = v(p) = (p -p * -)(p * + -p * -) -1 .
(2) The shape of γ β * ,p sensitively depends on the shape of the considered bounded volume and the boundary conditions. See the discussion in the next paragraph for details. (3) If β * > β * c but p < p * -or p > p * + , the overall cell density of the dominating cell type is already above the characteristic cell density p * + that would be expected in the segregated region. This situation is called super saturation. There are no clusters of the rare cell type observed macroscopically, because there are not enough cells of this type in the aggregate [START_REF] Dobrushin | Wulff construction. A global shape from local interaction[END_REF] (4) If 0 < β * < β * c , there is no phase segregation. The canonical Gibbs measures are almost product measures. This means that we typically observed disorder in the longtime limit.

Proof. Dobrushin et al. [1992, Thm.1.7, Thm.1.8] combined with [START_REF] Ioffe | Exact large deviation bounds up to T c for the Ising model in two dimensions[END_REF]; [START_REF] Ioffe | Dobrushin-Kotecký-Shlosman theorem up to the critical temperature[END_REF]. The main tool for the cited results is a surface order large deviation theorem for the defect of magnetization.

Impact of the boundary conditions

In the preceding paragraph, we have stated the results on phase segregation for periodic boundary conditions, thus neglecting the effects of potential cellboundary interactions. Now we explain what happens when we take such effects into account. For fixed boundary conditions, we also observe phase segregation. However the shape of the clusters as well as the decision of which cell-type sorts into the center of the volume depend on the boundary conditions. It will turn

A c c e p t e d m a n u s c r i p t

out that the geometry of phase segregation is very sensitive to the shape of the considered volume as well as to differences in the strength of cell-boundary interaction. We even observe a boundary phase transition, that is an abrupt change of the system's longtime spatial arrangement resulting from a smooth change at the boundary. Suppose that d = 2 and W = {0, 1, b}. We fix a volume V ⊂ S = Z d and consider the cell sorting model in V with fixed boundary condition as given in Definition 3.2. We allow only the two alternate states 0 or 1 at the lattice sites in V , thus we require that the cell number κ By Corollary 4.14, the longtime behavior of the cell sorting model with fixed boundary conditions is described by the canonical Gibbs measures μ V,ζ b ,κ as given in Definition 4.12. Consequently, we will consider phase segregation in the canonical Gibbs ensemble with fixed boundary condition. The relevant results in the mathematical literature that shall be applied are formulated in the context of thermodynamics. In particular, the theory on equilibrium crystal shapes as a result of the interatomic forces is relevant. These problems are studied with the help of Gibbs measures on the space E V := {-1, 1} V or E := {-1, 1} S . Therefore it is necessary to embed the measures on X V , in particular the canonical Gibbs measures μ V,ζ,κ , into the set P(E V ) of measures on E V . For this, we define a bijection F : X → E : η → ω η by

= (k 0 , k 1 , k b ) ∈ N 3 0 satisfies k 0 + k 1 = |V |
F (η)(x) = ω η (x) = 2η(x) -1, x ∈ V.
(4.17)

Obviously, the function

F -1 : E → X : ω → η ω satisfies F -1 (ω)(x) = η ω (x) = ω(x) + 1 2 , x ∈ V. (4.18)
The function F induces a bijection between P(X V ) and P(E V ) by

F (μ) := μ • F -1 , μ ∈ P(X V ) and F -1 (μ) := μ • F, μ ∈ P(E V ).
It turns out that the invariant measures μ V,ζ,κ of the cell sorting model are mapped by F onto the so-called canonical Ising models.

[-|V |, |V |], if μβ,γ,m (ω) = Z -1 β,γ,m exp -Hβ,γ (ω) 1 M -1 V (m) (ω), ω ∈ E V ,
where MV :

E V → R : M V (ω) = x∈V ω(x)
gives the excess magnetization,

Z β,γ,m := ω∈E V exp -Hβ,γ (ω) 1 M -1 V (m) (ω),
is the normalizing factor and

Hβ,γ (ω) := - β 2 x,y∈V, |x-y|=1 ω(x)ω(y) -γ x∈V,y∈V c |x-y|=1 ω(x) (4.19)
is the Ising Hamiltonian.

The basis for our discussion of the impact of the boundary conditions is the following embedding result. The proof of this proposition can be found in the Appendix, §AppendixA.6. This proposition allows us to apply the results of [START_REF] Shlosman | The droplet in the tube: A case of phase transition in the canonical ensemble[END_REF] and Bodineau, [START_REF] Bodineau | Rigorous probabilistic analysis of equilibrium crystal shapes[END_REF] on asymptotic phase segregation for canonical Ising models. All of these findings are proven rigorously in the mentioned literature. Because of the elaborate mathematical machinery which is necessary to state the results in full precision, we only sketch the main results.

Proposition 4.19. Suppose that κ = (k 0 , k 1 , 0) with k 0 + k 1 = |V | and define ζ b ∈ X by ζ b (x) = b, x ∈ S. Let μ V,
Fix

V = Λ N = [-N, N ] 2 ∩ Z 2 with N large and choose κ = (k 0 , k 1 , 0) ∈ N V . Define p = k 1 |V | -1 . Further, suppose that β * > β c and p ∈ (p * -, p * + )
, where p * ± are defined in (4.11). This is the parameter region where phase segregation is observed. We consider first the case of 1-boundary conditions, that is the case γ * = β * . In this regime the configurations with many type-1 cells close to the boundary are favored. Hence the typical configurations under μ V,ζ b ,κ show a large cluster of type-0 cells surrounded by type-1 cells. The 0-cluster does not touch the boundary of V , if V is large enough to contain the whole Wulff droplet of size p|V |. When the number of the 0-cells becomes so large that the Wulff droplet of this size cannot be placed inside the volume V , the shape of the 0-cluster is given by a modified Wulff construction. See Shlosman [1989, §1.4] for details. A qualitative picture illustrating both cases is given in Figure 2.

The case of 0-boundary conditions, where γ * = -β * , can be understood by symmetry arguments. When the roles of 0 and 1 are exchanged, then the cell In case of 1-boundary conditions, that is for γ * = β * > β * c , type-0 cells (dark grey) sort to the center of the aggregate. They are surrounded by type-1 cells (light grey). The left picture shows the sorting geometry if the type-1 cell density p satisfies p * -≤ p ≤ p * + but p is small enough that the whole Wulff droplet of size v(p) = (pp * -)(p * +p * -) -1 can be placed inside the unit square. In this case, the 0-cluster has a Wulff shape and does not touch the boundary. If p * -≤ p ≤ p * + , but the Wulff droplet of size v(p) does not fit into the unit square, the shape of the 0-cluster is given by a modified Wulff construction (right).

sorting model with type-1-cell density p and 0-boundary conditions is transformed into a cell sorting model with type-1-cell density 1p and 1-boundary conditions. The parameter β * does not change since it is symmetric in 0 and 1. Note that 1p ∈ (p * -, p * + ) if and only if p ∈ (p * -, p * + ). Consequently, for 0-boundary conditions, we see a type-1 cluster which is surrounded by type-0 cells. The shape of this cluster is again determined by the Wulff or modified Wulff construction.

The typical configurations that are observed in the long-time limit for empty boundary conditions, that is for γ * = 0, are shown in Figure 3. See Shlosman [1989, §1.4] for details.

If one relaxes the restrictions on the shape of V , that is if one allows that V is some bounded convex subset of Z 2 , then essentially the above results remain valid. However the shape of the Wulff droplet can change. See Bodineau et al. [2000, §IV] for details.

In principle, studies on phase segregation for spatially non-homogeneous boundary conditions are possible as well. However the Wulff construction requires the solution of a variational problem that might be very complicated. The case of rectangle-shaped V where the boundary condition at one wall is different from the type-1 boundary conditions at the other walls is considered in Bodineau et al. [2000, §IV]. It is shown that the shape of the type-0 cluster changes abruptly from the Wulff shape to the Winterbottom shape as the boundary conditions are varied. This phenomenon is an example for a boundary phase transition. For free boundary conditions, that is for γ * = 0, type-1 cells (light grey) separate from type-0 cells (dark grey) in the longtime limit, if the effective adhesion parameter is sufficiently high (β * > β * c ). The 1-cluster has the form of one fourth of the Wulff shape and is located in one of the corners of the unit square, if the type-1 cell density p is slightly bigger than p * -(left). If p is slightly below p * + , a type-0 cluster in one of the corners of the unit square is observed, which has the form of one fourth of the Wulff shape (right). For intermediate values of p, a serial arrangement is observed (middle).

Model Behavior

The mathematical analysis of our cell sorting model suggests that our model describes basic mechanisms in the process of sorting. In particular, phase segregation, which is interpreted as cell sorting, is observed in the longtime limit.

We have a fairly complete picture for two-dimensional, fully occupied twotype systems, that is for W = {0, 1} and d = 2. In this case, our cell system consists of cells of two different types where all sites are occupied by cells of either type. If the differential adhesion parameter β * = β 00 + β 11 -2β 01 is sufficiently high then the typical pattern arising in the longtime limit is characterized by two clear-cut separated regions consisting of virtually only one cell type.

The volume proportion of these two regions is determined by the ratio of type-0 and type-1 cell numbers while the asymptotic shape of the segregated domains is given by the solution of an isoperimetric problem, the Wulff shape. The latter is defined by that closed rectifiable curve which minimizes a functional, called surface tension or interfacial free energy, with respect to all closed rectifiable curves enclosing a fixed volume. The Wulff shape depends sensitively on the boundary conditions and the overall form of the given volume where the sorting takes place. The surface tension functional can be derived directly from the potential (4.7) of the cell sorting model and is therefore related to the individual cell-cell interaction based on intercellular adhesion. In this respect our cell sorting model provides a link how differential intercellular adhesion causes tissue surface tension. It shows that tissue surface tension can indeed guide the phenomenology of the cell sorting process as was postulated by the DAH

A c c e p t e d m a n u s c r i p t

1/2 if β * is close to β * c . If the density p of type-1 cells is not in the interval p * -, p * + , there is no phase segregation because the volume is statured by the dominating cell typ. Only for p ∈ p * -, p * + the system exhibits sorting with characteristic concentration p * + of the dominating cell type. Note that neither the critical differential adhesion parameter β * c nor the typical concentration p * + of the dominating cell type in the case of segregated volumes depend on the given concentrations of the two cell types. The latter factor solely influences the volume proportion of the segregated regions.

In summary, the model predicts a phase transition, that is an abrupt change from intermingled (unsorted) to segregated (sorted) longtime behavior, as the differential adhesion parameter β * increases. So the parameter β * regulates whether sorting occurs or disorder persists and it determines the characteristic cell density of the dominating cell type within the segregated regions. The shape of the segregated patterns for β * > β * c is controlled by the surface tension, which in turn is functionally related to the intercellular adhesion parameters. In contrast, the cell type that sorts to the center is determined by the order between the cell-boundary adhesion parameters β 0b and β 1b . These results on the qualitative longtime behavior can be summarized in a phase diagram (Fig. 5).

The migration parameter α neither plays any role for the classification of the longtime behavior nor for the geometry of the segregated patterns. However it influences the pace of relaxation to equilibrium. The higher the cell mobility is the faster the cell system converges to its stationary state provided that the other system parameters are kept constant.

Mathematically, the cell sorting model is a (generalized) stochastic lattice gas model. For a characterization of the qualitative longtime behavior, we could show that the system converges in the long run to a stationary measure which is even reversible and belongs to the class of Gibbs ensemble measures with respect to the (generalized) Ising potential. The study of essential properties of Gibbs measures is an active mathematical field of its own. For the case that W = {0, 1}, that is for the fully occupied two-type cell sorting model, we could exploit mathematical or statistical physics results on phase segregation for Gibbs measures w.r.t the Ising potential in the canonical ensemble.

This link between reversible measures and Gibbs measures could be made precise for the cell sorting model in an infinite volume as well. We showed that the reversible measures of the infinite cell sorting model coincide with the (infinite) Gibbs measures with respect to the generalized Ising potential. This presents the opportunity for ambitious further research on emergent dynamic properties which relies on the application of proper spatio-temporal limit procedures. Again the cell sorting model with W = {0, 1} is better understood, since in this case we know additionally that there are no non-reversible stationary states.

Further questions on stationary measures for the infinite model as well as phase transition and phase segregation in the cell sorting model with more than two cell types are still open. It turns out that hitherto existing methods for the analysis of such models heavily rest on the fact that there are only two possible states per spatial site.

Our cell sorting model is a stochastic particle based model which reflects our perception about individual adhesion-biased cell migration. Its emergent longtime behavior at the tissue level is characterized by two qualitatively different regimes: phase segregation and disorder. The mobility parameter does not influence the characteristics of the resulting pattern but guides the temporal scale of the pattern formation. Our cell sorting model is an IPS-model and therefore allows further rigorous mathematical analysis.

Discussion

Biological Foundation of the Modeling Assumptions

The proposed model is an interacting particle system which is constructed from our idea of adhesion-biased cell migration at the cellular scale. It is defined via the specification of local transition rates. Thereby our concept of the cell dynamics is (qualitatively) translated into the dynamical properties of the model system. We have concentrated on two influencing factors at the cellular scale, cell migration and intercellular adhesion, bundling additional influences into a stochastic component. The resulting Markov process model is the first genuine dynamical model of cell sorting as a result of intercellular interaction which can be analyzed rigorously. Although a caricature of the complex interactions within real biological cell aggregates, it is complex enough to reproduce cell sorting at the tissue scale. In this respect it is a 'minimal model', showing that our simplified picture of the processes at the cellular scale might contain the essential variables that explain cell sorting. On the other hand, the limited number of influencing variables allows a thorough mathematical analysis which reveals farreaching implications. And it opens up the road for further research on the dynamic properties of cell sorting. Our main hypothesis, that cell adhesion dampens cell mobility, that is the cells are less mobile the more they stick to each other, is formulated and analyzed on purely mathematical grounds without resort to arguable analogies with physical behavior.

Despite the undoubted importance of cell migration for cell sorting, there are only a few experimental or theoretical results, which describe the properties of individual cell mobility within cellular aggregates or estimate the impact of several types of intercellular interaction on cell locomotion. The intuition that individual cell mobility is lower the more a cell adheres to its neighbors seems to be plausible but shall be discussed explicitly. Our assumption is supported by the findings in [START_REF] Upadhyaya | Anomalous diffusion and non-Gaussian velocity distribution of Hydra cells in cellular aggregates[END_REF]. There Upadhyaya and coworkers monitored and analyzed the motion of single endodermal cells within homotypic endodermal and ectodermal cell aggregates. They show that the effective diffusion constant for individual endodermal cells in an endodermal environment is smaller than that in an ectodermal environment. Presuming that the strength of intercellular adhesion between ectodermal cells is less than that between endodermal cells [START_REF] Krieg | Tensile forces govern germ-layer organization in zebrafish[END_REF], we can conclude that cells are more motile in a less cohesive environment. Nevertheless it is not clear from this study that differences in the adhesive intercellular interaction cause the observed variation in single cell mobility. There might be other interactions between cells in contact which affect cell locomotion. Note that our model parameters β ij , i, j ∈ W , which we have called adhesion parameters out of habit, describe the effect of adhesive interaction as well as that of other short range intercellular processes with an impact on cell mobility. Here adhesion is certainly essential but other processes such as the modulation of cortical tension due to actin filament reorganization in intercellular zones might turn out to be relevant as well. So, on one hand, our model shows that quantitative differences in any shortrange intercellular mechanism which has an effect on cell locomotion results in cell sorting provided the differences are large enough. This confirms the critique by Harris [START_REF] Harris | Is cell sorting caused by differences in the work of intercellular adhesion? A critique of the Steinberg hypothesis[END_REF] of the DAH. On the other hand, we see that the determining which intercellular process is essential and which one is of minor importance is a problem that needs to be studied at the protein to cell level.

Here the challenge is to estimate the effect of molecular interaction between cells in contact on the ability of the cells to detach and move away. Once there are sufficiently high contact-specific differences in cell mobility, sorting occurs regardless of the microscopic sources of these variations.

Considering the concept of intercellular adhesion in cell sorting studies, we find that the prevalent picture is dominated by a static view which is rooted in the analogy to physical adhesion and cohesion, for instance between liquid molecules. The latter is the result of van der Waals forces, hydrogen bonds or electrostatic interaction and is characterized by the following three important properties [START_REF] Harris | Is cell sorting caused by differences in the work of intercellular adhesion? A critique of the Steinberg hypothesis[END_REF]. Firstly, physical adhesion is equivalent to close range attraction, that is the forces that pull the molecules of a liquid together are the same that hold them together. Secondly, physical adhesions are reversible bonds, which means that the breaking of an adhesive binding is just the reverse of its formation. Thirdly, the existence of adhesive bonds does not alter the internal structure of the involved molecules. Intercellular adhesion differs from adhesion between liquid molecules since these properties do not hold in general for intercellular adhesion in each mentioned respect. First, several adhesive bonds between cells develop after cells come into contact, so that the cells are held together by additional forces compared to those which pulled them together. The overall strength of intercellular bindings increases with time (on a scale that is relevant for cell migration) [START_REF] Krieg | Tensile forces govern germ-layer organization in zebrafish[END_REF]. Secondly, the ability of a cell to detach from a neighboring cell might follow different rules and is, in general, not just the reverse of the bond-forming process. And thirdly, intercellular adhesion is mediated by cell surface molecules which function in various other cellular processes [START_REF] Halbleib | Cadherins in development: cell adhesion, sorting and tissue morphogenesis[END_REF][START_REF] Marsden | Integrin-ECM interactions regulate cadherin-dependent cell adhesion and are required for convergent extension in Xenopus[END_REF]. This means that the formation of adhesive bonds can change the intracellular structure and thereby alter several cellular properties such as cell locomotion characteristics. Therefore we suggest abandoning the physical analogies and adopting a dynamical viewpoint by directly assessing the effect of intercellular bindings on individual cell mobility.

Our consequent dynamical viewpoint allows us to overcome another chal-

A c c e p t e d m a n u s c r i p t

lenge connected with the DAH and several existing cell sorting models such as [START_REF] Glazier | Simulation of the differential adhesion driven rearrangement of biological cells[END_REF]. Employing the liquid-tissue analogy, it is assumed there that cell aggregates are closed thermodynamic systems, since the applied results on minimization of surface tension in liquids rely on the principle of energy conservation. However cells can produce energy by themselves, so cell aggregates are thermodynamically open. Therefore, even when one postulates the existence of some 'adhesion energy' or 'surface tension energy' for cell aggregates, it cannot be assumed a priori that cell migration is such that the 'free (adhesive or interfacial) energy' is minimized. Instead, this property might result at the tissue scale from the characteristics of the individual cell dynamics. Besides, 'adhesion energy' or 'surface tension energy' are properties at the tissue scale. So the postulate that cell migration is directed towards the minimization of these quantities implies that each individual cell must gain information about the spatial arrangement within the whole cell aggregate in order to decide on the next migration step. Are there cell-intrinsic mechanisms that can accomplish this task? If one excludes global diffusive signalling, then cells can only collect information about their local environment and behave according to this information. Therefore, the real question is whether individual cell behavior based on local information actually leads to global minimization of the proposed quantities. This cannot be answered theoretically until one specifies a dynamical model where the individual cell behavior within cell aggregates is based on biologically realistic assumptions about the local intercellular interaction.

By proposing an IPS-model we have chosen a well-established model class where powerful tools exist for the analysis. IPS models require that the space is discretized, which is a disadvantage with regard to both the biological situation and some technicalities in the model specification and analysis. However, currently the mathematical theory for analogous modeling problems in continuous space is not developed enough to set up and analyze a cell sorting model within this framework [START_REF] Kondratiev | On non-equilibrium stochastic dynamics for interacting particle systems in continuum[END_REF]. It remains open -and is actually a subject of current mathematical research -to what extent the spatial discretization affects the model predictions.

The model can easily be amplified to incorporate more detail at the cellular scale. For instance, concerning cell shape and orientation, we implicitly assume that the cells have a spatially symmetric form (spheres or cubes) and that they do not show a preferred orientation of movement. If one suspects that the cells polarizations affect cell sorting, than these details can be included into the cell sorting model by changing the local state space to e d , -e d }}, where e i denotes the i-th unit vector of R d , i = 1, ..., d, with the interpretation that the pair (k, v) represents a cell of type k oriented in direction v. Adjusting the transition rules such that they reflect the assumed effect of orientation on cell migration one can start to analyze the model. Another generalization of the model would be to allow larger interaction neighborhoods, thus reproducing the effect of filopodia in intercellular adhesion. Of course, the more detail is included, the harder the analysis becomes. So far, only models with comparatively simple local state space W and small interaction neighborhood can be analyzed rigorously.

W = {0} ∪ {(k, v) : k = 1, ..., n, v ∈ {e 1 , -e 1 , ...,

A c c e p t e d m a n u s c r i p t

Since our cell sorting model is still well-defined for unbounded volumes, it becomes possible to explore its dynamical behavior at characteristic scales. There exist already several mathematical results in this direction, for instance [START_REF] Bertini | Stochastic interacting particle systems out of equilibrium[END_REF], but for a detailed application to the cell sorting model these results must be developed further.

Model Predictions for Cell Sorting Behavior

Concerning the analysis of the longtime behavior in the cell sorting model we would like to discuss two main results. Firstly, there is a phase transition from unordered to ordered patterns in the longtime limit as the effective parameter β * increases. Since the results of our analysis are of qualitative nature, the critical value β * c cannot be quantified precisely but it is clear that the condition 2β 01 < β 00 +β 11 does not automatically imply that cell sorting occurs. This fact is often ignored in studies concerning cell sorting [START_REF] Shi | Biophysical properties of cadherin bonds do not predict cell sorting[END_REF][START_REF] Steinberg | Reconstruction of tissues by dissociated cells[END_REF]. In particular, whenever the impact of specific binding molecules on cell sorting is to be estimated one should take into account the existence of this threshold. Note further that disorder is found in the longtime limit for each β * = β 00 + β 11 -2β 01 that is non-negative but below the critical value β * c , but not in the parameter regime 2β 01 < β 00 + β 11 . In the latter case, the effective differential adhesion parameter β * is negative, which means that we have essentially repulsion between unlike cell types. The patterns that are expected in the long run are of checkerboard type [START_REF] Blöte | Accurate determination of the critical line of the square Ising antiferromagnet in a field[END_REF][START_REF] Van Den Berg | A uniqueness condition for Gibbs measures, with application to the 2-dimensional Ising antiferromagnet[END_REF][START_REF] Dobrushin | Phase diagram of the two-dimensional Ising antiferromagnet (computer-assisted proof)[END_REF][START_REF] Georgii | Percolation and number of phases in the twodimensional Ising model[END_REF], and can therefore be easily mistaken with unordered patterns. However, for checkerboard-type patterns, the degree of order, measured for instance by the absolute values of the spatial correlations, is high. When one requires that β * > 0, then just the relative distances between the adhesion parameters β ij count, but not the position of zero within the ordered chain. This means that our results on cell sorting apply as well to a wide parameter regime with intercellular repulsion.

Secondly, once we see phase segregation for β * > β * c , an important question is which cell type sorts into the center of the cell aggregate. The DAH predicts that it is the cell type with the higher homotypic intercellular adhesion. However recent studies disagree with this hypothesis [START_REF] Krieg | Tensile forces govern germ-layer organization in zebrafish[END_REF]. Our model suggests that the type-specific differences in the strength of cellular interaction with the surrounding material determine the sorting hierarchy in two-type sorting experiments. This agrees with the experimental findings in [START_REF] Ninomiya | Epithelial coating controls mesenchymal shape change through tissue-positioning effects and reduction of surfaceminimizing tension[END_REF], where it was reported that otherwise identical cell aggregates show different sorting hierarchy when one aggregate is wrapped by an epithelial layer. In our model, the parameters β 0b , β 1b which control the sorting hierarchy quantify the overall effect of those processes between the cells at the boundary of the cell aggregate and the surrounding medium which have an impact on cell locomotion. Here adhesive interaction is an important candidate but other mechanisms such as the interface-specific intracellular modulation of actomyosin filaments might have an effect as well [START_REF] Krieg | Tensile forces govern germ-layer organization in zebrafish[END_REF]. Depending on the properties of the surrounding medium, the processes that dominate the

A c c e p t e d m a n u s c r i p t

cell-medium interactions can be different from those that are important for the effect of intercellular contact on cell locomotion. Note further that the shape of the underlying volume where the sorting takes place influences the geometry of the sorted pattern. Therefore the study of sorting within a growing cell aggregate might produce unexpected results which cannot be inferred directly from the static volume geometry. Another result of our analysis is the discovery that the temporal evolution of the system is directed towards the minimization of a functional that can be interpreted as surface tension. That means, starting with an idea about the impact of the local environment on the individual cell dynamics, we can describe the emergent cooperative longtime behavior. So the postulates of the DAH on the existence of tissue surface tensions and on the ability of cell aggregates to arrange themselves in configurations of minimal surface tension are now results of the analysis instead of assumptions. The surface tension functional is defined in relation to the adhesion parameters without to resorting to physical analogies.

The mobility parameter α does not influence the characteristics of the phase segregation patterns observed in the longtime limit. However it does affect the pace of the sorting process. The process of phase segregation is slower for smaller α. Note that sorting is also slowed down by an increased β * . Therefore, for small α and very large β * , the process of cellular rearrangements might not reach the fully sorted regime within biologically relevant time scales. In this case, we expect only partially sorted cell aggregates, where small to medium-sized cell clusters of one type are scattered within a sea of alternate cell types.

Applied Methods for the Model Analysis

We have focused on the rigorous analysis of the two-dimensional fully occupied two-type cell sorting model (W={0,1}, d=2). For this situation, we show that our model is equivalent to a stochastic lattice gas with Ising potential and exploit existing mathematical tools and results in this field. In the case where more than two cell types are involved in a cell sorting experiment or when the system is three-dimensional, the proposed model is still applicable. A mathematically pure analysis is much more complicated and might be accompanied by an empirical study of the model implications. In this respect it is important to notice that the temporal development of an IPS-model can be easily simulated on a computer by suitable time discretization. The trajectories of the resulting time-discrete Markov chain show the same longtime behavior as the continuous-time cell sorting model. See [START_REF] Klauss | Modelling and simulation by stochastic interacting particle systems[END_REF] for details.

We do not expect qualitatively different results for the analysis of our cell sorting model in the three-dimensional case. Cerf and Pistora showed in [START_REF] Cerf | On the Wulff crystal in the Ising model[END_REF] that also for d ≥ 3 the typical configurations that appear in the long term are those with minimal surface tension, where again the surface tension is a functional at tissue scale that can be derived from the intercellular adhesion coefficients.

A c c e p t e d m

a n u s c r i p t

Summary

Our model bridges the gap between ideas concerning cell migration and intercellular interaction at the cellular scale that are currently debated and tissue scale predictions on the longtime behavior of the whole cell aggregate which were put forward by the DAH. By concentrating on simplified dynamic properties, we observe the existence of a phase transition in the onset of sorting, highlight the role of cell-boundary interaction in the cell sorting hierarchy and derive a concept of surface tension that is related to the biological properties at the cellular scale. We conclude that the interaction between the cells and the boundary of a confining vessel has a major impact on the sorting geometry. model allows further studies of the emergent dynamic characteristics of the cell sorting process.

  , which causes one cell to increase its volume by one lattice site and another cell to decrease its volume, accordingly. This volume fluctuation induces shifts in each cell center of mass and was therefore claimed to represent random cell motion.

Figure 1 :

 1 Figure 1: Type-specific cell migration of on a square lattice is biased by nearest-neighbor intercellular adhesion.

  Note that this means essentially that we identify the points of Λ N with the points of the discrete d-dimensional torus of diameter N . Replacing in the above equations (3.1) and (3.2) the Manhattan norm |.| by the norm |.| Λ N , we define a family c N p = (c N p (x, y, •))x,y∈S of transition rates on X N := W Λ N . For each initial state, the family c N p of time-homogeneous transition rates uniquely determines a Markov chain on X. In detail, let be C(X N ) the set of all real functions on X N and define a Markov p : C(X N ) → C(X N ) by the equation

  Hence y =o c(o, y, η) < ∞. Property (F) follows now from the translation invariance (A). Proposition (D) follows from the observation that c(x, y, η) ≥ min j∈W \{0} α j exp -4d max i,j∈W

  is referred to as the canonical partition function. The sum is taken over the set O T (ζ) = {η ∈ X : η S\T = ζ S\T , N T (η) = N T (ζ)} of those configurations that satisfy the outside condition ζ S\T and have got the same particle numbers as ζ inside of T . Definition 4.4. Suppose that Φ is a finite range potential. A probability measure μ ∈ P(X, F) is a canonical Gibbs measure w.r.t. Φ, if

  states Now we ask whether non-reversible invariant measures exist for the cell sorting model on S = Z d . This question can be answered for at most twodimensional, fully occupied systems, that is for W = {0, 1} and d = 1 or d = 2. Theorem 4.7. Let W = {0, 1} and suppose d = 1 or d = 2. If the family of finite range transition rates c = (c(x, y, •)) T ∈T is given by (3.2), then

  . Next we show that the canonical Gibbs measures are reversible measures of the cell sorting model with fixed boundary condition.

  ζ,κ and one obtains the following corollary of Theorem 4.10.

Definition 4. 15 .

 15 The canonical ensemble (μ β * ,N,R N ) N ∈N exhibits asymptotic phase segregation along some curve γ, if the measures μβ * ,N,R N weakly converge to m(β * )μ γ for any sequence R N satisfying (4.13).

  and k b = 0. The boundary of V shall be composed solely of typeb cells, that is we concentrate on spatially homogeneous boundary conditions ζ ∈ X, where ζ b (x) = b, x ∈ S. The parameters β 0b and β 1b describe the strength of the cell-boundary interaction. Let us define γ * := 2β 1b -2β 0bβ 11 + β 00 . (4.16) We distinguish several special situations. If γ * = β * , where β * is defined in (4.6), we have a cell sorting model with 1-boundary. Note that we have γ * = β * if, for instance, β 0b = β 01 and β 1b = β 11 . If γ * = -β * , for instance for β 0b = β 00 , β 1b = β 10 , we have a cell sorting model with 0-boundary. In case that γ * = 0, we have a cell sorting model with free boundary.

  A measure μβ,γ,m on E V is a canonical Ising model with parameters β, γ ∈ R and excess magnetization m ∈

  ζ b ,κ denote the canonical Gibbs measure w.r.t. the potential Φ given by (4.4). Then F (μ V,ζ b ,κ ) = μβ,γ,m , where m = k 1k 0 , β = β * /4, γ = γ * /4 and μβ,γ,m is a canonical Ising model.

Figure 2 :

 2 Figure 2: Phase segregation for 1-boundary conditions and varying type-1 cell density p.

Figure 3 :

 3 Figure 3: Phase segregation for free boundary conditions and varying type-1 cell density p.

Figure 5 :

 5 Figure 5: Phase diagram for the longtime behavior of the cell sorting model. The effective adhesion parameter β * and the density p of the dominating cell type determine the longtime behavior of the cell sorting model. For 0 < β * < β * c , the typical configurations that are observed in the longtime limit are disordered (intermingled states). Above the critical value β * c , phase segregation is observed for 0.5 ≤ p < p * + (β * ). In the parameter region β * > β * c , p > p * + (β * ), the dominating cell type outnumbers the alternate cell type so much that phase segregation is impossible (super saturated states). The green line is the graph of p * + (β * ) for β * > β * c .

  Now assume that x, y ∈ S, |x -y| = 1, y, η)h Φ {x,y} (η) = c 0 (x, y, η) expβ η(x)η(z) + z:|y-z|=1 β η(y)η(z)β η(x)η(y) ⎫ ⎬ ⎭ =c 0 (x, y, η) exp -β η(x)η(y) = (α η(x) + α η(y) )p(x, y) exp -β η(x)η(y) , if η(x) = η(y), 0, otherwise.The last term depends on η only via N {x,y} (η). Therefore the function c(x, y, .)h Φ {x,y} is E {x,y} -measurable and the detailed balance condition is satisfied.Proof of Theorem 4.5. The assertion of Theorem 4.5 follows directly from Lemma AppendixA.2, Lemma AppendixA.3 and Proposition 3.3 (E).AppendixA.2. Detailed balance condition w.r.t ΦLemma AppendixA.4. Let the family of transition rates c = (c(x, y, •)) x,y∈S be given by (3.2). Then the family c of transition rates satisfies the detailed balance condition w.r.t. the finite range potential Φ = ( ΦT ) T ∈T given by (4.7).Proof.If |x -y| = 1 then c(x, y, η) = c(x, y, η xy ) = 0, hence c(x, y, η)h Φ {x,y} (η) = c(x, y, η x,y )h Φ {x,y} (η x,y ). Now assume that x, y ∈ S, |x -y| = 1 and η ∈ X. For any u, v ∈ S with |u -v| = 1, we have Φ {u,v} (η) = -β η(u)η(v) = -β 00 (1η(u))(1η(v))β 01 (1η(u))η(v) β 10 η(u)(1η(v))β 11 η(u)η(v) = -(β 00 -2β 01 + β 11 )η(u)η(v) + (β 00β 01 )(η(u) + η(v))β 00 = Φ{u,v} (η) + (β 00β 01 )(η(u) + η(v))β 00 . because Q is E V -measurable. Therefore μβ,γ,m (F (η)) = μ V,ζ,κ (η), η ∈ X V . Since F (μ V,ζ b ,κ ) = μβ,γ,m if and only if μ V,ζ b ,κ = F -1 (μ β,γ,m ) = μβ,γ,m • F ,this implies the assertion. We show (A.3). For η ∈ X V , we find-H Φ V (τ V (ζ b , η)) 00 + 2β 01 + β 11 ) + (ω η (x) + ω η (y))(β 11β 00 ) + ω η (x)ω η (y)β * ] + 1 2 x∈V,y∈V c |x-y|=1 [ω η (x)(β 1bβ 0b ) + β 0b + β 1b ] , since 4β η(x)η(y) =4β 00 (1η(x))(1η(y)) + 4β 01 [(1η(x))η(y) + η(x)(1η(y))] + 4β 11 η(x)η(y) =β 00 (1ω η (x))(1ω η (y)) + β 01 [(1ω η (x))(1 + ω η (y)) + (1 + ω η (x))(1ω η (y))] + β 11 (1 + ω η (x))(1 + ω η (y)) =β 00 (1ω η (x)ω η (y) + ω η (x)ω η (y)) + β 01 (1ω η (x) + ω η (y)ω η (x)ω η (y)) + β 01 [(1 + ω η (x)ω η (y))ω η (x)ω η (y)] + β 11 (1 + ω η (x) + ω η (y) + ω η (x)ω η (y)) =(β 00 + 2β 01 + β 11 ) + ω η (x)(-β 00 + β 11 ) + ω η (y)(-β 00 + β 11 ) + ω η (x)ω η (y)(β 00 -2β 01 + β 11 ) =(β 00 + 2β 01 + β 11 ) + (ω η (x) + ω η (y))(-β 00 + β 11 ) + ω η (x)ω η (y)β * and 2β η(x)b =2β 0b (1η(x)) + 2β 1b η(x) =β 0b (1ω η (x)) + β 1b (1 + ω η (x)) =ω η (x)(β 1bβ 0b ) + β 0b + β 1b . (x)(β 1bβ 0b ) + β 0b + β 1b ] (x)(β 1bβ 0b ) + β 0b + β 1b ] ,

The parameters β 0b and β 1b specify the strength of the interaction between a cell of the respective type and the boundary. The value γ * = 2β 1b -2β 0bβ 11 + β 00 determines the sorting hierarchy. Type-0 cells (dark grey) wrap type-1 cells (light grey) for γ * = β * (left) but sort to the center of the aggregate if γ * = -β * (middle). For γ * = 0 a serial structure is observed (right).

[ [START_REF] Steinberg | On the mechanism of tissue reconstruction by dissociated cells. I. Population kinetics, differential adhesiveness and the absence of directed migration[END_REF].

Considering the question of which cell type sorts into the center of the given volume, we find that the boundary conditions are fundamental. The situation at the boundary is characterized by the adhesion parameters β 0b and β 1b which specify the strength of adhesive interaction between a cell of the respective type and the boundary. The model predicts that the question of which cell type sorts into the center of the volume can be determined from the value of the parameter γ * = 2β 1b -2β 0bβ 11 + β 00 . For γ * = β * , we observed that a type-0 cluster is surrounded by type-1 cells, while type-1 cells sort to the center of the aggregate for γ * = -β * . In the case that γ * = 0 only partial engulfment is observed (see Fig. 4).

In contrast to the phase segregation observed for a sufficiently high differential adhesion parameter, there is a parameter region where disorder persists. Below a critical value β * c , there are just unorganized states possible in the longtime limit, that is our cell sorting model predicts that there is no cell sorting for these parameter values. This means in particular, that the difference (β 00 + β 11 ) -2β 01 has to be larger than the threshold value β * c in order to achieve phase segregation and that it is not sufficient to require that β * is just positive.

We find an intermediate behavior for a differential adhesion parameter β * only slightly above the critical value. When we observe asymptotic phase segregation in this parameter regime, then the decomposition is not complete but characterized by a typical concentration p * + = p * + (β * ) of the dominating cell type. This coefficient p * + depends exclusively on the differential adhesion parameter β * and not on the given total cell numbers. Its value is slightly above

Proof. It follows from Georgii [1979, Thm. 1.21 and Remark 1.28] that G c (Φ) = ∅ for any finite-range potential Φ. The proposition R(c) ⊃ G c (Φ) as well as the equality R(c) = G c (Φ) for c irreducible are implied by Georgii [1979, Thm. 2.14].

Lemma AppendixA.3. Let the family of transition rates c = (c(x, y, •)) x,y∈S be given by (3.2). Then the family c of transition rates satisfies the detailed balance condition w.r.t. the finite range potential Φ = (Φ T ) T ∈T given by (4.4).

Proof.

A c c e p t e d m a n u s c r i p t

Hence

The last term depends on η only via (η(z)) z∈S\{x,y} and N {x,y} (η). Therefore H Φ {x,y} -H Φ {x,y} and h Φ {x,y} /h Φ {x,y} are E {x,y} -measurable. Thus we conclude from Lemma AppendixA.3 that c satisfies the detailed balance condition with respect to Φ.

AppendixA.3. Proof of proposition 4.13

Since X V is a finite discrete space, any continuous function on X V is a finite linear combination of indicator functions 1 {ξ} , ξ ∈ X V . Thus, by the linearity of A V ζ , it is enough to show that

Using (3.8), this is equivalent to

Hence the above equation is trivially satisfied for any pair η, ξ ∈ X V with ξ ∈ {η xy : x, y ∈ V, |x -y| = 1}. It remains to show that

A c c e p t e d m a n u s c r i p t

for any η ∈ X V and x, y ∈ V . The validity of this equation is inferred from the following consideration. Fix x, y ∈ V . From the definition (4.1) of the Hamiltonian H Φ V we find that

where the summand

Applying Lemma AppendixA.3 and using (4.2), we obtain

Hence

. By Definition 4.12, this is exactly the equation (A.2).

AppendixA.4. Canonical Gibbs measures with periodic boundary conditions

We recall that the cell sorting model with periodic boundary conditions is described in Definition 3.1. We consider cubes Λ N = [0, N -1] ∩ Z d and identify the points of Λ N with the points of the discrete torus of diameter N . For N ∈ N, we define the map τ * N :

for any y ∈ N Z d + x. Since there are only two alternate states per lattice site, we have for the particle numbers defined in (4.3) that N T (., 0) = |T | -N T (., 1), T ∈ T . Therefore let us agree to abbreviate N T (η) := N T (η, 1), η ∈ X N .

A c c e p t e d m a n u s c r i p t

Definition AppendixA.5. The probability measure μ β * ,N,R ∈ P(X N ) is a canonical Gibbs measure w.r.t. Φ in Λ N with periodic boundary conditions and type-1-cell number R, if

AppendixA.5. Embedding of the discrete torus into the continuous torus

We recall that the points of the cube Λ N = [0, N 1 ] ∩ Z d are identified with the points of a discrete torus of diameter N . This torus can be embedded into the continuous torus T = R d /Z d . The norm in T is given by |x| T := min{|y| : y ∈ x + Z d }. We use M( T ) to denote the space of bounded Borel measures on T equipped with the topology of weak convergence. We transform measures on the torus Λ N into measures on (M( T ), B(M( T ))), where B(M( T )) denotes the Borel σ-algebra of M( T ). For this, it is useful to represent configurations as measures on T . Namely, denoting by δ x the Dirac measure supported by the point x ∈ T , we associate with each configuration η ∈ X N the measure

Thus a measure measure μ N on X N = W Λ N is transformed into a measure μN on M( T ) by μN (A) = μ N ({η : μ η ∈ A}), A ∈ M( T ).

AppendixA.6. Proof of Proposition 4.19

For η ∈ X V , we observe that F (η) ∈ M -1 V (m) if and only if η ∈ N -1 V (k 1 ), since MV (F (η)) = 2N V (η) -|V |, η ∈ X V and 2k 1 = m + |V |. Therefore it is sufficient to show that

where H Φ V is defined in (4.5) and Q : X V → R is some E V -measurable function. Indeed, we obtain from (A.3) that μβ,γ,m (F (η)) =Z -1 β,γ,m exp -Hβ,γ (F (η))

Comparing this formula with the Definition 4.12, we find that