Poissonian cellular Potts models reveal nonequilibrium kinetics of cell sorting

Cellular Potts models are broadly applied across developmental biology and cancer research. We overcome limitations of the traditional approach, which reinterprets a modified Metropolis sampling as ad hoc dynamics, by introducing a physical timescale through Poissonian kinetics and by applying principles of stochastic thermodynamics to separate thermal and relaxation effects from athermal noise and nonconservative forces. Our method accurately describes cell-sorting dynamics in mouse-embryo development and identifies the distinct contributions of nonequilibrium processes, e.g. cell growth and active fluctuations.

The dynamics of many nonequilibrium systems can be described by a time-dependent phenomenological Hamiltonian which actively controls transitions through a sequence of target states.Widely adopted frameworks of vertex, cellular Potts, and other methods rely on this effective energy-based principle to explain spatial organization in living systems [1][2][3][4][5][6][7][8][9][10][11].
Whereas the optimum of the system's energy specifies a target state of such an active transformation, the unfolding of the modeled process in time is determined by its kinetic parameters.In vertex or subcellular-element models with a continuous phase space, these parameters correspond to the transport properties-the damping coefficients.However the traditional cellular Potts models (CPMs), whose discrete-state dynamics are implemented by a modified Metropolis sampling, lack an explicit control over such kinetic parameters.
Transport properties and the timescales they control are especially important when multiple processes evolve interdependently.In the course of embryonic development, numerous cellular and tissue-level processes require precise mutual coordination [12].For example the sorting of cell types in the early mouse embryo must be completed before the subsequent morphogenetic events commence [13,14].
To introduce kinetic parameters into CPMs we invoke the theory of stochastic thermodynamics [15,16], which comprehensively describes discrete-state physical processes driven by changes of free energy.As shown
further, the transport properties control the system's frenetic activity [17], which constitutes the time-symmetric component of a stochastic action-the complement of entropic changes of a system's trajectory.While being less demanding than the subcellularelement method, CPMs can treat composite materials and describe more intricate shapes than vertex models [4,[18][19][20][21].Each cell in three dimensions corresponds to a contiguous collection of voxels with the same "spin" value-labels distinguishing individual objects in the system [Fig.1(a)].
CPMs were first introduced by Graner and Glazier [4] to study how differences of surface energy between homotypic and heterotypic contacts cause cell sorting in ).Tau-leap simulations in discrete time (DT, dt = 10 −4 arb.u.) match the exact distribution (Theory); p-value of the multinomial test 0.997 [22,23].The continuous-time simulations (CT) and the Metropolis algorithm (MA) produce comparable results (p-values 0.966 and 0.986 respectively).In the DT and CT models the chain is inhomogeneous with action rates (α2i+1 = 0.1 for the odd indices, α2i = 0.3 arb.u. for the even).Error bars are given by three standard deviations of 10 4 realizations.(b) Relaxation of the chain energy between two equilibrium states with average values ⟨H⟩0 and ⟨H⟩∆ at temperatures T (t < t0) = T0 = 1.8|J| and T (t ≫ t0) = T0 + ∆T = 2.0|J|/kB respectively (J = −1 arb.u.).The results of MA sampling are reported alongside the CT simulations of a slow (α2i = 0.05 and α2i+1 = 0.08 arb.u.) and fast (α2i = 0.5 and α2i+1 = 0.7 arb.u.) kinetics.Each curve traces an average over 10 4 trajectories with a standard-error band.
development.Using the modified Metropolis algorithm they found that clusters of cells emerge in a typical configuration favored by the system's energy function in which the first sum runs over spin pairs σ i and σ j with symmetric coefficients J ij (σ i , σ j ) = J ji (σ i , σ j ) encoding the surface interactions, whereas the second term penalizes deviations of the volume V k of the k th cell from its preferred value Vk .Usually J ij (σ i , σ j ) are identically zero unless the spins σ i and σ j are in direct contact.As the method of Graner and Glazier [4] evolved beyond a mere proof of concept, it was further generalized to include nonequilibrium aspects, such as cell division and active motility [4,11,[18][19][20][21][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40].
The modified Metropolis algorithm of modern CPMs is not the most general kinetic model for discrete systems [41][42][43][44][45] and has limitations [31,46].In fact, the Metropolis scheme was originally designed to bypass slow simulations of systems' dynamics when sampling equilibrium ensembles [47,48,Chapter 7].Since the early days of CPMs, questions have therefore been raised regarding the interpretation of time, temperature, dissipation effects, and nonequilibrium aspects of this approach [46].
As discrete-state Markovian systems, which describe continuous-time physical phenomena, CPMs can be most naturally regarded as Poissonian processes (see Sec. 1 in Supplemental Materials).The framework proposed here leverages this result of queuing theory [49] to address the above questions.Through Poissonian kinetics and stochastic thermodynamics we introduce interpretable time and energy scales, account for response co-efficients and forces not incorporated in the Hamiltonian, and finally separate thermal and athermal fluctuations.
Framework.-As an illustration of our approach we first consider a paradigmatic example of discrete systems and a special case of the Potts model [50]-an Ising chain σ = (σ 1 , σ 2 , ..., σ N ) with nearest-neighbor interactions.Arranged on a one-dimensional lattice, N spins σ i=1,2...N ∈ {−1, 1} with a periodic boundary condition σ N +1 = σ 1 are described by the Hamiltonian with an interaction constant J.
To define dynamics of the Ising chain we assume that each spin flips its sign with a Poissonian transition rate k i (σ), which in general depends on the current state σ [Fig.1(b)].Within a sufficiently small time dt at most one spin can change its value.In equilibrium, the detailed-balance condition for such a spin in which H(σ i ) and k i (σ i ) are respectively the spin's energy and transition rate, given the values of σ j̸ =i .From Eq. (2) it follows then The chain's stochastic kinetics is then given by a master equation, once a common factor between k i (σ i ) and  3) is specified [41,[51][52][53].In a general context each spin may be characterized by a statedependent action rate α i (σ i ), which determines the probability for the i th spin to attempt a sign change.When such an attempt occurs, the transition probability p(σ i → σ ′ i ) is determined by a directing function L(σ ′ i ) with a normalization constant Z [54] For the two possible outcomes-no change, σ ′ i = σ i , and a transition, σ ′ i = −σ i -the normalization constant expands to Z = e L(σi) + e L(−σi) and Eq.(4) yields With the transition rate given by the product of the attempt rate and the transition probability we find from Eq. ( 3) A complete specification of the Ising chain now requires both the Hamiltonian and the spins' action rates.Such a system can be simulated exactly in continuous time by the standard techniques for master equations, or approximately by using a tau-leap algorithm with a step dt [55].
The action rates do not compromise the canonical distribution of the Ising chain in equilibrium [Fig.2(a)].These parameters control the unfolding of dynamical processes, such as relaxation of transients.For example, chains with larger action rates, initially prepared in equilibrium at temperature T 0 and subject to a sudden temperature change ∆T , relax to the new steady state faster [Fig.2(b)].By design, the Metropolis scheme renders samples of the target equilibrium ensemble after a very short transient trajectory, which can not be controlled by algorithmic or system parameters.
Model analysis.-Toanalyze the Poissonian dynamics of the Ising chain we apply the theory of stochastic thermodynamics.Any given trajectory of the system θ from an initial state σ 0 to a final state σ M can be decomposed into a sequence of M elementary paths Each path T consists of n quiescent intervals of arbitrarily small time dt in the same configuration σ, followed by a sign change of the i th spin which produces the next state σ ′ .The probability of this change is p i ≃ k i (σ i )dt, whereas the probability of the quiescent period lasting n steps is With these definitions, the probability of the elementary path from a given initial condition is Now we can decompose the stochastic action A of the elementary path into the entropic and frenetic components [17], ∆S and D respectively: Indeed the probability of a time-reverse trajectory T -a change of the i th spin conditioned on the initial configuration σ ′ and followed by n quiescent steps-is in which p ′ i = k i (−σ i )dt.Due to Eqs. ( 7) and ( 8) the time-asymmetric part of the action yields The time-symmetric part renders a more involved expression approximated by for small dt.The total action of the whole trajectory is Without compromising the entropic activity, action rates control the system's frenesy through transition rates k j , cf.Eqs ( 7), ( 14), and (16).The entropy change of a relaxation process is entirely determined by the energy difference between the initial and final configurations of the system [Eq.( 13), Fig. 2(b)].In contrast, the frenesy depends on the kinetics of each state transition in a system's trajectory.
Poissonian cellular Potts models.-Wenow construct a three-dimensional Poissonian CPM [56].Voxels on a cubic lattice describe the state of K distinct cells and a medium, taking values σ i ∈ {0, 1, 2, ..., K} (Fig. 1).The coefficients J ij (σ i , σ j ) in Eq. ( 1) vanish when voxels i and j are both occupied by the same object, or when the voxel i is not within the Moore neighborhood of the voxel j [31].Otherwise J ij assume constant positive values encoding the surface interactions between objects.
For each of the total ν object types, our framework introduces a Poissonian state-dependent action rate α(σ i ) ∈ {α 0 , α 1 , ..., α ν }, with which an i th voxel attempts to change its current value σ i .Its possible target values σ (j) i are chosen from the Von Neumann neighborhood like in the standard CPMs [31], with the transition probabilities given by a general version of Eqs. ( 4) and ( 8) e ∆L(σ (j) i ) = α(σ Temperature in traditional CPMs is a fictitious parameter manipulated to adjust the level of fluctuations [57].In contrast, our approach regards it as a physical variable that is set at an experimentally controlled value.To prevent cell fragmentation, usually suppressed by a periodically applied annealing, we adopt the local-connectivity test of Durand and Guesnet [31] in a modified form (see Sec. 3 in Supplemental Materials).
An extension of the directing-function formalism [54, Appendix B] can also incorporate more general nonequilibrium forces into Eq.(17) as in which the active exponents ϕ(σ i , σ i ) are functions associated with a specific transition σ i → σ (j) i , and in general ϕ(σ i , σ i , σ i ).When this perturbation breaks the detailed-balance condition, such a transition incurs an irreversible thermodynamic work.
Active exponents can introduce non-Hamiltonian forces and, by violating the fluctuation-dissipation theorem, athermal noise (see Sec. 2 in Supplemental Materials).Here we focus on noise amplification by nonequilibrium processes, which are usually present inside cells, such as the chemically-driven polymerization of cytoskeletal filaments or molecular motor activity [58][59][60].If we set all transitions σ i → σ (j) i , except for the trivial ones σ (j) i ≡ σ i are promoted.Nonconservative forces are not generated by active exponents, whose asymmetric components vanish ϕ(σ i , σ Specific active processes can be modeled more explicitly as well.Cell growth is typically implemented by a time-dependent preferred volume Vk in Eq. ( 1).Persistent cell motility can be incorporated by additional terms of the Hamiltonian [32,[34][35][36]61], or by asymmetric active exponents.Cell sorting during embryonic development.-Asa biophysical example we consider the sorting of epiblast (EPI) and primitive endoderm (PrE) cells in the early mouse embryo.These cells form the inner cell mass (ICM) aggregate and sort into an outer single layer of PrE cells separating the epiblast from the medium [62] [Fig.3(a)].
Recent advances provide unprecedented experimental access to the dynamics of isolated ICMs [63][64][65][66][67].We quantify the segregation of the two cell types by a sorting score computed from distances of EPI and PrE cells (r EPI i and r PrE j respectively) from their common geometric center: in which N PrE and N EPI are the numbers of PrE and EPI cells.By definition the score s ∈ [−1, 1] is close to zero for unsorted cells, and −1 or 1 when all PrE cells are inside or outside the aggregate respectively.To model the sorting process, we chose the five interaction constants J medium:EPI , J medium:PrE , J EPI:EPI , J EPI:PrE , J PrE:PrE from a physiologically relevant range of the EPI and PrE surface tensions, set the temperature to the experimental value at 310.15 K, and calibrated the growth parameters to match the observed proliferation dynamics (see Sec. 3 in Supplemental Materials).
Almost perfect sorting is achieved within 480 min of CPM simulations for a wide range of parameters [Fig.3(b)].The action rates control the relaxation dynamics of the sorting process [Fig.3(c Furthermore, our simulations predict distinct shape dynamics of the two cell types, which quantitatively agree with the experimental data in a parameter-free comparison of cellular aspect ratios [Fig.3(d)]: EPI cells tend to more rounded shapes, whereas PrE cells stretch normally to the radial direction at the outermost shell of ICM.For more details on ICM sorting in vivo see Ref. [67].
Faster kinetics of PrE cells promote sorting [Fig.4(a)].This effect of inhomogeneous action rates can neither be modeled in the traditional CPMs (see Sec. 1 in Supplemental Materials), nor can it be compensated by time rescaling, as the curves parameterized by α 0,EPI,PrE in general belong to different families [Fig.4(a)].
Without cell growth and division [Fig.4(b)], or active fluctuations [Fig.4(c)] sorting is hindered.Both mechanisms do a thermodynamic work on the system: the growth of cells generates stresses, and new cell boundaries increase the total surface energy, whereas the active fluctuations inject energy by breaking detailed balance.Responding to these nonequilibrium processes, the system rapidly acquires the energetically favored sorted state.Modeling active fluctuations by an effective temperature is also viable [68,69], but may misrepresent the system's response to thermodynamic forces (see Sec. 2 in Supplemental Materials).
In fact, the response coefficients of cells' surfaces are directly related to action rates (see Sec. 2 in Supplemental Materials).Because energy is a derived unit of time, "independent" rescaling of time or temperature, afforded by the traditional CPMs at the account of fictitious energy scales, is forbidden once the experimental data fix k B T and J ij .
Conclusions.-PoissonianCPMs provide a physically consistent framework to study complex materials with active properties, which is generally applicable to other discrete-state systems [51][52][53]70].Its kinetic parameters control transport coefficients and permit an unambiguous interpretation of time.Active fluctuations and nonequilibrium processes are clearly separated from thermal effects and passive relaxation.We applied this framework to examine the roles of distinct nonequilibrium processes in embryonic cell sorting, and show that either growth and division, or active shape fluctuations are required for successful segregation of cell types.Equation (S23) represents an implicit constraint imposed by the traditional CPMs through the choice of g 0m and A 0m , which is not dictated by any physical laws or empirical observations.Historically this choice was convenient to sample equilibrium ensembles when the dynamical details are ignored.However, there are infinitely many Poissonian processes, whose dynamics generates the same ensemble, and among which the traditional CPMs select implicitly only one.
Our framework affords a more general form of a Poissonian process given by Eqs. ( 17) and (18).In equilibrium conditions, these processes generate a Boltzmann distribution with the designated Hamiltonian, whereas the action rates α i (σ) can be arbitrary nonnegative functions of the state vector σ.In the following, we assume that the Hamiltonian H(σ) is given.With all voxels' states except σ i fixed, we denote for brevity α 0 = α i (σ i = s 0 =i ).Then i th voxel changes its current state with an escape probability [cf.Eqs.(17)] in which we introduced weights w m = exp(∆L m0 ).The individual rates can then be expressed as which allows us to rewrite Eq. ( 17) in the form Furthermore, because the transition rate k 0m and its reverse k m0 satisfy the detailed-balance condition, we can derive a formula for the weights, which is analogous to Eq. ( 18): with κ m being the escape rate of the state σ i = s m .As follows from Eq. (S25) there exists a time constant such that τ 0 k 0m = w m .Hence from Eq. (S24) we can express the action rate in terms of τ 0 as In summary the above procedure requires as input the system's Hamiltonian and all the transition rates and yields the escape rates κ j=0,1...M and the time constants τ j , which then determine all action rates α j .Whereas the dynamics of traditional CPMs implies a particular choice of the action rates α(σ), which can be reconstructed by following the above procedure, our approach offers their general functional form α i (σ), which can be adjusted to relevant physical assumptions or learned from empirical data.In the main text we have proposed to set α i (σ) = α i (σ i ), which relies on local information about the current state of the i th voxel.As discussed in the next section, this choice implies a constant response coefficient of cell surfaces to the thermodynamics forces.
Notably, the general theory of the Metropolis algorithm allows one to use any rule of the acceptance probability A jm , as long as it ensures the detailed balance condition p jm /p mj = (g jm A jm )/(g mj A mj ) = exp(β∆H jm ).Using this property, one may replace the traditional rule (S22) by the softmax function .
Then g jm = g mj can be interpreted as the attempt probability with homogeneous action rates Under these conditions the transitions pjm = g jm Ãjm satisfy the detailed balance condition [cf.Eqs. ( 17) and ( 18)].Hence, the Metropolis scheme can be modified to recapitulate special cases of the Poissonian dynamics, e.g. with homogeneous action rates α.S5.Schematic of a surface element perpendicular to the x axis at the interface between two states (s0 and s1), which undergoes a random walk.

Response coefficients and nonconservative forces
The action rates α i (σ) and active exponents ϕ(σ i , σ (j) i ), introduced in the Poissonian framework of CPMs, can also be interpreted in the context of cell-surface mechanics [38][39][40].In particular, action rates together with the spatial discretization give rise to the response coefficients of cell-surface elements, whereas the active exponents lead to nonconservative forces-not derived from the Hamiltonian-and amplify fluctuations, as shown below.
Consider two voxels, σ 0 = s 0 and σ 1 = s 1 , which share one face perpendicular to the spatial coordinate axis x, in different states s 0 ̸ = s 1 (Fig. S5).Such a face represents an interface, e.g. between two cells, which is positioned with a probability p(t, x 0 ) = 1 at time t and coordinate x 0 .If at time t + dt the voxel σ 0 (σ 1 ) copies the state from σ 1 (σ 0 ), the interface is moving left (right) on the discrete lattice to a new position x −1 = x 0 − ∆x (x 1 = x 0 + ∆x).For simplicity we assume that within some larger 3D neighborhood of σ 0 and σ 1 there are no voxels in other states s / ∈ {s 0 , s 1 }.Then we can express the probability of finding the interface at a given coordinate x i at time t + dt as .

(S28)
Given that the spatiotemporal scales ∆x and dt are sufficiently small, the following expansions may be applied: Equation (S26) defines a random walk, which has been thoroughly analyzed in the context of diffusion problems [71][72][73][74][75].As we are following closely the notation of Ref. [54], we may immediately borrow its results to express into which we can substitute Here we decomposed ϕ(s 1 , s 0 ) = 2Φ + βF ∆x and ϕ(s 0 , s 1 ) = 2Φ − βF ∆x into symmetric (frenetic) and asymmetric (entropic) parts, Φ and F ∆x respectivly, in which the latter represents the work done by the force F when the interface moves right.In the continuous limit ∆x → 0, dt → 0, and Eqs.(S31)-(S33) yield a Fokker-Planck equation [54] with the diffusion constant emerging as α(s 0 )α(s 1 )∆x 2 /2 → D, and the gradient ∆H/∆x → ∂ x H.This diffusion equation stands in correspondence to a Langevin dynamics of the form [76] ẋ with the mobility coefficient µ = βD and the standard Brownian motion dB [⟨dB⟩ = 0, ⟨dB(t)dB(s)⟩ = δ(t − s)].
A few remarks about Eqs.(S34) and (S35) are worth emphasizing: • The action rates together with the spatial discretization ∆x give rise to the response coefficient µ between thermodynamic forces and their conjugate currents.
• These equations contain a nonconservative force F , which does not stem from the Hamiltonian H and emerges due to the antisymmetric part of the active exponents.This is a nontrivial extension of the CPMs afforded by the Poissonian framework.Our derivation also gives a statistical-physics justification of the conservative forces −∂ x H considered in Refs.[38][39][40].
• The frenetic part of the active exponents amplifies noise and, thus, breaks the fluctuation-dissipation theorem.This effect can be summarised by the effective diffusion coefficient • Equations (S34) and (S35) generalize the Smoluchowski equation and the Brownian dynamics in a potential H(x) [54], which emerge in the limit of vanishing nonequilibrium forces.Note, that in 1D the active force F can always be effectively reincorporated into the Hamiltonian [54, Appendix B], but in more dimensions not all forces can be expressed through gradients of an effective Hamiltonian.
The relation between the diffusion coefficient D ∝ α∆x 2 and the spatial discretization ∆x has an important consequence.If one changes the spatial resolution ∆x = ξ∆x, the action rates should also be rescaled to α = α/ξ 2 to preserve the transport properties of cell surfaces.
Lastly, we discuss the concept of effective temperature [68,69], which one may invoke to account for the fluctuationallowance parameter in the traditional CPMs.If we were to simulate the CPM at the effective temperature T eff = k B T e Φ , with the detailed balance implied, the fluctuation-dissipation theorem would be restored at the expense of altering the response coefficient µ eff = D/(k B T eff ) = µe −Φ , and thus either speed up or slow down the descent of the system along the energy gradients.This approach therefore misrepresents the kinetic details of the modeled systems.
However, by using an effective temperature, Poissonian CPMs are capable of extending the equilibrium fluctuationdissipation theorem for the cell-surface dynamics consistently through manipulation of the transport coefficients.In particular, we set T eff = k B T e Φ , and rescale all action rates α eff = αe Φ/2 .Thereby the effective level of noise is amplified through D eff = D(α eff ) = De Φ , whereas the mobility of cell-surface elements remains unaffected Whereas the rescaling outlined above preserves the mesoscopic dynamics of cell surfaces by design, other mesoscopic and macroscopic variables may still be affected by these modifications [69].

Computational details
In all simulations of CPMs we used a 3D voxel grid with a resolution ∆x 3 = 1 µm 3 and a discrete tau-leap algorithm with a time step 0.1 min [55].To convert surface tension constants of the EPI and PrE cells with the medium, γ medium:EPI and γ medium:PrE respectively, to the interaction parameters J medium:EPI and J medium:PrE we used a formula J c = 6∆x 2 γ c /26 for the contact type c, which can be interpreted as follows.A single voxel in our model has 6 faces of area a = ∆x 2 and interacts with 26 sites in the surrounding Moore neighborhood.The surface energy of such a voxel in a medium is then E surface = 6aγ c = 26J c , from which the formula for J c follows.By using the same formula the tension constants of cell contacts γ EPI:EPI , γ EPI:PrE , γ PrE:PrE can also be converted to the spin interaction constants J EPI:EPI , J EPI:PrE , and J PrE:PrE respectively, and vice versa.The total surface area of a cell was energetically constrained [77].
Because experimental observations at our disposal are not sufficient to fit parameters of our CPM, we had to estimate and adjust manually their values to achieve a similar dynamics of the sorting score in simulations.Thus we chose J medium:EPI = 30.0×10 3 k B T and J medium:PrE = 9.2×10 3 k B T , which are close to the maximum experimentally observed tension γ medium:EPI = 669 pN/µm and the minimum observed tension γ medium:PrE = 178 pN/µm (see Experimental details).Such a large difference of tension is necessary to compensate additional sorting mechanisms, which are not accounted for by our model and are discussed in the companion paper [67].The cell-cell tension constants were chosen in the interval (γ medium:PrE , γ medium:EPI ) to favor the interface energy of PrE cells with the medium over EPI and homotypic contacts over heterotypic ones: J EPI:EPI = 18.5 × 10 3 k B T , J EPI:PrE = 25.4 × 10 3 k B T , and J PrE:PrE = 13.8 × 10 3 k B T .The experimental temperature is 37 • C.
The volume elasticity constant κ in Eq. ( 1)is perhaps the most difficult parameter to choose, because we do not have experimental observations of cell-volume dynamics.A very large elasticity constant leads to artificial cuboidal cell shapes, whereas cells with a small elasticity and preferred volume may collapse to a single voxel due to the surface tension.Furthermore cells with a large value of κ grow faster (see growth and division below).Combined with a large action rate such cells may acquire unrealistic shapes with long protrusions and no bulk volume.Empirically we found that κ = 1154 k B T /∆x 6 allows us to observe realistic cell shape with action rates below 3.57 min −1 .
At the present quite little is known about regulation of the cell mitosis and death in the inner cell mass of the mouse blastocyst.Experimentally we could however measure volumes of a few mitotic cells (see Experimental details).Their mitotic cycle lasts approximately 12 h [67].In the view of this information we decided to neglect the cell death and use the following model of growth and division.
For each cell we sample a division volume Ṽ from the empirical distribution found experimentally [Fig.S6(a)].The cell's preferred volume Ṽ is growing with a constant rate g until the cell's actual volume V reaches or exceeds Ṽ .The mother cell is then divided approximately into two halves perpendicular to the major principal axis of the cell's gyration tensor [78][79][80].From experiments we found that average of ICM cells is about 1500 µm 3 , whereas the average volume of mitotic cells is approximately 3000 µm 3 -the double of the average cell volume.Hence the cells should roughly double in size within one cell cycle.Given the chosen elasticity constant κ, we adjusted the growth rate of the preferred volume g = 3.36 µm 3 /min to match the 12-hours duration of the mitotic cycle.The initial conditions for our simulations were generated from a single cell to match the observed median value of 25-26 cells per ICM at day E3.5, by setting the cell-medium and cell-cell tension constants and the cell action rate to the average values: γ medium:cell = (γ medium:EPI + γ medium:PrE )/2, γ cell:cell = (γ EPI:EPI + γ PrE:PrE + γ EPI:Pre )/2, α cell = (α EPI + α PrE )/2.The cell fates were assigned randomly with 60% PrE cells [67].

Sorting score
Action rates were adjusted to match the sorting-score dynamics in experiments.Examples of dynamics modulated by the spanned range of action rates is illustrated in Fig. S6(b).
To prevent cell fragmentation we adopted the method of Durand and Guesnet [31] in a modified form.As we discussed with one of the two authors [81], their approach does not guarantee that the cells remain simply connected in 3D models.In order to check that the cell is locally connected at a given site, we find all voxels occupied by the cell in the adjacency neighborhood of this site: all possible configurations of the adjacency neighborhood can be reduced to a total of 10 symmetrically nonequivalent situations to consider (Fig. S7).One special case illustrated by Fig. S7(h) passes the connectivity test described in Ref. [31], but changing the value of the central site in this configuration would create a hole.Therefore we prohibit such changes in our simulations.In the other situations we apply the method of Ref. [31] in the original form.
The above modification of the connectivity test prevents a trivial scenario, in which one cell would "drill" a handle in another cell, whereas the connectivity of cells is always ensured by the test of Durand and Guesnet [31].For the The central site has the same value as four of the six voxels in its adjacency neighborhood.Whereas the case (g) requires the continuity test to ensure that changing the value of the central site is allowed, in the case (h) changes are not allowed, even though the red voxels in the adjacency neighborhood are connected.Panels (i) and (j): The central site has the same value as five and all six voxels of its adjacency neighborhood respectively.In the latter case no changes are allowed.
parameter values of interest we did not observe in our simulations appearance of holes or handles.Nonetheless our modification may not guarantee that such geometric defects would not occur for some extreme values of parameters.In 3D systems it appears that a hole at a given site could be created through a change of adjacent voxels.

Cell tension measurements
Micropipette aspiration was performed as described in [82] to measure the surface tension of ICM cells.In brief, a micro-forged micropipette coupled to a microfluidic pump (Fluigent, Microfluidic Flow Control System) was used to measure the surface tension of ICM cells.Micropipettes were prepared from glass capillaries (Warner Instruments, GC100T-15) using a micropipette puller (Sutter Instrument, P-1000) and a microforge (Narishige, MF-900).A firepolished micropipette with a diameter of 7-8µm was mounted on an inverted Zeiss Observer Z1 microscope with a CSU-X1M 5000 spinning disc unit, and its movement was controlled by micromanipulators (Narishige, MON202-D).Samples were maintained at 37 • C with 5% CO 2 .A stepwise increasing pressure was applied on ICM surface cells using the microfluidic pump and Dikeria software (LabVIEW), until a deformation with the same radius as that of the micropipette (R p ) was reached.The equilibrium pressure (P c ) was measured, images were acquired in this configuration and then the pressure was released.At steady state, the surface tension γ at the cell-fluid interface is calculated based on Young-Laplace's law: γ = P c /2(1/R p − 1/R c ), in which P c is the net pressure used to deform the cell of radius R c .Image analysis and measurement of the pipette radius R p and R c was done in FIJI, and calculation of surface tension was done using custom scripts in Python v3.9.

Live-imaging of externalized inner cell mass
For externalization of the inner cell mass from mouse blastocysts, the zona pellucida (ZP) was first removed from blastocysts with pronase (0.5% w/v Proteinase K, Sigma P8811, in global medium containing HEPES (LifeGlobal,

FIG. 3 .
FIG. 3. Sorting of epiblast (EPI) and primitive endoderm (PrE) cells in ICMs isolated from mouse blastocysts.(a) Immunostaining images on embryonic day 3.5 and 8 h later (see Sec. 4 in Supplemental Materials).(b) Cellular-Potts simulation in the initial configuration and 8 h later.(c) Sorting score of live-imaging data averaged over 8 experiments compared to a mean trend of 500 simulations [α0,EPI,PRE = (0.95, 0.69, 2.12) min −1 ]; standard-error band shown only for the experiments.(d) Average cellular aspect ratios characterized in the maximum-projection planes of experimental data and by the eigenvalues of the 3D gyration tensors in the simulations.

kFIG. 4 .
FIG. 4. Kinetic and nonequilibrium aspects of Poissonian CPMs.(a) Convergence of sorting log(1 − s) for different values of heterogeneity δα = (αPrE − αEpi)/α0, α0 = 1 min −1 ].Inset: Scores at t = 480 min show that faster PrE kinetics promotes sorting.Heterogeneous kinetic properties produce non-rescalable differences in the sorting dynamics-note the two-point intersection between the black dashed curve with homogeneous action rates (α0 = αEPI = αPrE = 0.2 min −1 ) and the green curve (δα = −0.5 min −1 ).(b-c) Either cell growth and division, or active fluctuations ( φ), or extreme values of an effective temperature kBTe are required for complete sorting.Curves in panels (a) and (b-c) represent averages over 1000 and 500 trajectories respectively [parameters set to values as in Fig. 3(c) unless noted otherwise].
), Fig. S1(b) in Supplemental Materials].We sampled 100 combinations of the values {α 0 , α EPI , α PrE }, with each entry chosen from the interval (0.10, 3.57) min −1 .The best match is closest to the experimental curve in the least-squares sense.
FIG. S7.Scheme of the 3D adjacency neighborhood in the connectivity test of Durand and Guesnet [31] for 10 symmetrically nonequivalent configurations.Panel (a): the central site (red) contains no voxels of the same value in the six sites of the adjacency neighborhood (gray).If the red color represents a cell, the value of the central site cannot be changed, because the corresponding cell would disappear from the simulation.If the red color represents medium, in our simulations such a change may be admitted, although it would violate the detailed balance conditions and must be therefore interpreted as an active process, e.g.absorption of a fluid by cells.Panel (b): The central site has the same value as one of the six voxels in its adjacency neighborhood.Panels (c) and (d): The central site has the same value as two of the six voxels in its adjacency neighborhood.Panels (e) and (f): The central site has the same value as three of the six voxels in its adjacency neighborhood.Panels (g) and (h):The central site has the same value as four of the six voxels in its adjacency neighborhood.Whereas the case (g) requires the continuity test to ensure that changing the value of the central site is allowed, in the case (h) changes are not allowed, even though the red voxels in the adjacency neighborhood are connected.Panels (i) and (j): The central site has the same value as five and all six voxels of its adjacency neighborhood respectively.In the latter case no changes are allowed.