Adjusting the range of cell–cell communication enables fine-tuning of cell fate patterns from checkerboard to engulfing

During development, spatio-temporal patterns ranging from checkerboard to engulfing occur with precise proportions of the respective cell fates. Key developmental regulators are intracellular transcriptional interactions and intercellular signaling. We present an analytically tractable mathematical model based on signaling that reliably generates different cell type patterns with specified proportions. Employing statistical mechanics, We derived a cell fate decision model for two cell types. A detailed steady state analysis on the resulting dynamical system yielded necessary conditions to generate spatially heterogeneous patterns. This allows the cell type proportions to be controlled by a single model parameter. Cell–cell communication is realized by local and global signaling mechanisms. These result in different cell type patterns. A nearest neighbor signal yields checkerboard patterns. Increasing the signal dispersion, cell fate clusters and an engulfing pattern can be generated. Altogether, the presented model allows us to reliably generate heterogeneous cell type patterns of different kinds as well as desired proportions. Supplementary Information The online version contains supplementary material available at 10.1007/s00285-023-01959-9.


Introduction
Cell fate decisions play an essential role in establishing cellular function during development.In this process, previously indeterminate cells specify themselves into one of several different cell types.In many cases, there is a strong correlation between gene expression patterns and subsequent cell fate.Therefore, it is necessary to understand the dynamics of different genes to unravel the secrets of differentiation.
One prime example of this differentiation process is the differentiation towards neural and epidermal cells in Drosophila.Characteristically, epidermal cell progenitors express high levels of transmembrane protein Notch, whereas neural progenitors express low levels of the same [1,2].A similar example is found in the inner cell mass (ICM) of the preimplantation mouse embryo.There, the transcription factors NANOG and GATA6 have been identified as the earliest markers for the segregation of the ICM into epiblast and primitive endoderm cells, respectively [3,4].Apart from the spatial cell fate distribution, the correct cell fate ratio is also of particular interest [5,6,7].
In mathematical models, cell fate decisions are often described by systems of ordinary differential equations (ODE) based on a gene regulatory network (GRN).At the single cell level, toggle switches as models of interactions of two genes have been investigated in great detail [8,9].These represent mutual inhibition of two proteins combined with auto-activation.As a result, three stable steady states arise with regard to gene expressions that represent the different cell fates.It depends on the initial conditions which state a cell will be attracted to.At the tissue level, experimental studies hint towards the importance of paracrine signals with regards to differentiation [10,11].
Lateral interaction models have already found their way into the current research landscape.For the Delta-Notch signaling pathway, patterns of alternating cell types have been reconstructed [12].For the mouse embryo, models including cell-cell communication due to fibroblast growth factor signaling have been employed to create similar salt-and-pepper/checkerboard patterns [13,14].So far, these studies are concerned with an averaged nearest neighbor signal, i.e. cells do not communicate beyond their nearest neighbor.Further studies suggest that in fact cell fate patterning in the mouse embryo is the result of a complex interplay of cell signaling, cell division, cell sorting and apoptosis [15,16,17].
Mathematical modeling allows untangling the individual components and investigating their pattern formation potential.It was previously shown that cell division alone yields cell fate clusters [18].Simulations of cells sorting due to differential adhesion have been shown to generate engulfing patterns [19].This resembles the result of the minimization of the total contact energy [20].Here, we focus on the potential of intercellular signaling.In addition to nearest neighbor signaling, we consider signaling that can reach further across a tissue.This builds upon previous ideas for Drosophila [21,22,23] as well as the mouse embryo [24,25].In addition to that, they activate themselves.The signal s is an external factor activating v.
Based on methods from statistical mechanics [26,27,28], we derived a model describing the temporal development of the expressions of two genes.A generalized signal incorporates external influences on cell fate decisions.Performing a detailed stability analysis of the ODE system, we obtained necessary conditions in the form of a parameter interval to always generate a mixture of two different cell types in a tissue.Numerical simulations for an averaged nearest neighbor signal as well as a distance-based signal demonstrate the potential of our model to establish different spatial cell fate patterns ranging from checkerboard via clustering to engulfing patterns.To quantify the different resulting patterns, we employed individualized pair correlation functions (PCFs).A cell type proportion analysis revealed which proportions our model can create, but also which restrictions there are.Our work introduces an easy to control mathematical model for gene expression and our analysis results provide insight into signaling driven pattern formation and cell type proportioning.

Protein interaction model
In this section, we derive a model to describe cell fate decisions.As a basis for this we choose methods from [26,27,28] which allow us to describe transcriptional regulation on the level of the DNA.We consider a simple system of two different transcription factors u and v together with an external signal s describing the cell-cell communication.To this end, we consider a gene regulatory network (GRN) characterized by the mutual inhibition of u and v, as well as their auto-activation and the signal s activating v and inhibiting u (Fig. 1).

Gene regulation
To describe the dynamical system underlying transcriptional regulation, we consider two basic assumptions: 1. Transcription determines the production of new protein.

Decay describes the lifetime of the protein.
These assumptions are translated into a generic ordinary differential equation (ODE) describing the concentration of a protein u over time: The second term is the exponential decay with decay rate γ u .The first term describes the rate of transcription of the corresponding gene.Here, p u denotes the probability that RNA polymerase (RNAP) is bound to the promoter of u.The production rate r u describes how much protein can be produced while RNAP is bound.

Binding probability
Following [28,26,27], we consider the different binding events of a gene regulatory network (GRN).However, we assume that the auto-activatory part of u is dominant, such that the base activity of the RNA polymerase will be neglected.This means that the production of u mainly depends on its binding close to its own promoter.Now the system can be in two different states.Either u is bound or it is not.First we count the number of possibilities how these states might arise.We divide our space into Ω different lattice sites and describe the total number of protein via U = uΩ.The binomial coefficients yield the number of possible states Number of unbound states: Assuming different energies whether a protein is unbound ε unbound u or bound ε bound u , the two states have total energies Using Boltzmann statistics, the energy of the two states enables us to describe the probability that the system is in either of these states via e −βε unbound and e −βε bound .The partition function is given by the sum of all possible Boltzmann weights over every microstate, i.e.
Using the partition function, we are able to calculate the binding probability p u by the ratio of bound states Z bound and all states combined as Assuming Ω U , we use the approximation Ω! (Ω−U )! ≈ Ω U .We divide the numerator and denominator of (9) by Z unbound and define the energy difference ∆ε For simplicity, we introduce the energy coefficient η u := e −∆εu and use u = U/Ω to get again the volume fractions.This leads to

Interactions
The crucial parts in transcriptional regulation are the interactions between constituents.In the following, we consider that an additional species v interacts with the promoter of u.This results in a system, with the following microstates: The binding energy differences remain as before with an additional factor for the interaction η uv = e −∆εuv .The binding probabilities for U and V are then given by The advantage or disadvantage given by the interaction energy difference now determines the nature of the interaction.For η uv = 1, (12) can be simplified using factorization The binding probability reduces to the case without interaction.Consequently, the cases where η uv = 1 describe binding probabilities that are either lower or higher than the case with no interaction: • η uv = 0 ⇔ ∆ε uv = ∞: complete inhibition / blocking Case η uv = 0 was listed separately, because it represents a special case of inhibition in which u and v cannot be bound at the same time.

Describing the cell fate decision between two fates
We imagine a system, where two antagonistic proteins u and v are the deciding factors for the decision of a cell's fate.Therefore, u and v mutually inhibit each other.We use a blocking type of inhibition, i.e. the promoter of u is not active as soon as v is bound in the vicinity of u's promoter and vice-versa.
We can also interpret this as u not being able to bind if v is already bound.Both u and v are assumed to be dominantly auto-activating, such that the base activity of the RNAP can be neglected.Finally, an external signal s influences u and v in different ways.It activates v by cooperatively binding with v.At the same time, u is inhibited by s and the cooperative binding of u and s.We assume that for both promoters, the respective energy coefficients are equal.The above considerations lead to interaction coefficients Any single bound state results in the terms η α α with α ∈ {u, v, s}.The remaining state has v and s bound simultaneously, yielding the term η v η s η vs vs.For the binding probability of u, we collect all the terms including u and divide them by the combination of all other terms, resulting in Likewise, the probability of v is given by Using (1) together with (15) and ( 16) for a total of N different cells, we end up with a system of ordinary differential equations (ODEs) We note that so far, the cell-cell interactions are not further specified.This means that the absorbed signals of each cell s i are provisionally considered as a generalized function of the expression values of all cells s : 3 Steady State Analysis

Existence of steady states
In order to get a better understanding of our ODE system, we want to delve further into the resulting steady states of the system.This means, we consider Consequently, we get When rearranging (19) and (20), we find two possible solutions for u i and v i , respectively.These solutions are Taking every combination of u i and v i from ( 21) into account, we end up with four different steady states.For three of the steady states, we can get either no expression of u and v or high expression of one transcription factor and none for the other: These steady states share the lower bound 0. Additionally, a rough estimate for an upper bound is given by the ratios of reproduction and decay r u /γ u and r v /γ v .For parameter combinations such that the left hand sides of the inequalities provide a reliable estimate for the steady state values.The fourth steady state is an oddity that arises by combining the non-zero solutions for u i and v i from (21).When combined, the corresponding variables u i and v i cancel out and we find the relation This also leaves our system to be over-determined and the values of u i and v i cannot further be identified.However, by using (26) in the steady state solution v i = 0 in (21), we obtain the following state: Isolating s i in equation ( 26) leads to a critical signal value s * for which this steady state will always occur This critical signal value is also responsible for a switching behavior in our system (Fig. 2).For values below or above s * , a cell ends up in states ( 23) (u + v − ) and ( 24) (u − v + ), respectively.At exactly s * , u and v move towards the straight line defined by (27) with no unique steady state.Altogether, we have successfully identified the relevant steady states ( 22)-( 24) of our ODE system (17) as well as the condition to force a switch in the cell's fate.

Linear stability analysis
In the following sections, we investigate the steady states in further detail.We employ linear stability analysis to determine the parameter regime that allows us to find a desired steady state for the overall system.At the single cell level, we rule out (22), since it is not relevant to cell fate specification.At the tissue level, we distinguish between homogeneous and heterogeneous steady states.A homogeneous equilibrium state consists of cells of a single type only.This means that either all of the cells in the tissue are in state ( 23) (u + v − ) or all of them are in state (24) (u − v + ).Best case scenario, is a mixture of the two cell types.Therefore, we aim at excluding the homogeneous steady states as well.We follow the definition of linear stability for an ODE system   23), ( 27) and (24).
We say, an ODE system is linearly stable in x * , if its linearization matrix L ODE = f (x * ) has only eigenvalues with negative real part.Using the N -dimensional identity matrix I N , we can write the linearization matrix of (17) as Using the chain rule, the block matrices A xy , x, y ∈ {u, v} can be written in terms of the partial derivatives where we define ∂pu ∂u := ∂pu ∂uj , (u i , v i , s i ) i,j=1,...,N . The other block matrices are defined analogously.For our purposes, we only need to focus on the following derivatives As usual, the eigenvalues of matrix L ODE are defined as the roots of the characteristic polynomial where I 2N denotes the identity matrix in 2N dimensions.At first glance, this determinant seems impossible to calculate.However, when inserting the respective steady states, we are able to reduce the matrix tremendously.

Excluding steady state (22)
In the following, we elaborate on how to exclude the first steady state (22) as solution for our ODE system (17).Without loss of generality, we assume u 1 = 0 = v 1 .This way, in row N + 1 all entries but one of the matrix L ODE − λI 2N become 0. The remaining entry with index (N + 1, N + 1) is Laplace expansion then enables us to write the determinant of the whole matrix as a product of (39) and the determinant of the remaining submatrix.Thus, it suffices to focus on the first eigenvalue given by This translates to the eigenvalue λ being Although the signal thus far has not been further specified, we propose a realistic physical representation by assuming s i ≥ 0. Furthermore, we consider an activation of v by the signal s, i.e. η vs > 1 and therefore, inequality and consequently provide the necessary condition for instability.The exclusion of this steady state strengthens our focus on ( 23) and (24), which represent the two different cell types u + v − and u − v + , respectively.

Instability of tissue-wide homogeneous steady state (23)
With steady states (23) and (24), we aim to find a parameter region for which we achieve a heterogeneous steady state, i.e. we get a tissue with a mixture of cells in the two states.To this end, we derive conditions for instability of the homogeneous steady state.We start with state (23) and set u i = ru γu − 1+ηssi ηu and v i = 0 for all i.Inserting these expressions into the derivatives (32)-(37) results in a simplification of L ODE .Since (34) and (37) are zero for every i, j, the off-diagonal block matrix A vu = 0.This means the determinant is given by the product of the determinants of the block matrices on the diagonal.Again, since (37) is zero, A vv becomes a diagonal matrix with diagonal entries Inserting Using this, we determine N factors of the characteristic polynomial N eigenvalues are given by the second factor in (45).For instability, it is sufficient that only one of these is greater than zero.In other words, this results in the inequality After appropriate rearranging, we obtain a sufficient condition for our parameters At this point, the general case cannot be simplified further.Depending on the cell-cell interaction and therefore the incoming signal s i , one can find an even more accurate description of this relation.Alternatively, we can formulate this condition in terms of energy differences as which allows us to see the maximum allowed deviation of the difference between ∆ε u and ∆ε v .Keep in mind that for this condition, we only relied on the first N eigenvalues.In truth, this condition might be even more relaxed than what we derived.
3.5 Instability of tissue-wide homogeneous steady state (24) We set u i = 0 and v i = rv γv − 1+ηssi ηv(1+ηsηvssi) .Using the same approach as before, we find that (33) and (36) are zero for all i, j and thus A uv = 0.In addition to that, we get a diagonal matrix for A uu .For u i = 0, its diagonal entries are Inserting v i yields As before, this allows us to determine N factors of the characteristic polynomial We exploit again the instability condition that any eigenvalue must be positive and find the inequality This yields another condition for η u .As before, it is necessary to fulfill this inequality for a single value s i , i.e. the minimum of all possible signal values suffices in that regard Again, we write this in terms of energy differences

Steady state summary
The stability conditions (47) and (54) define an interval for −∆ε u , with The reproduction rates r u , r v and decay rates γ u , γ v shift this interval by ln ruγv rvγu .The length of the interval is determined by the minimum and maximum signal values combined with the associated energy differences −∆ε s and −∆ε vs .The results of our stability analysis are summarized in figure 3.At the single cell level, we are able to exclude u − v − cells using inequality (41).Therefore, at the tissue level, we can distinguish between three different states.The stability interval (55) yields the exact parameter regime for the transition of the homogeneous states to the heterogeneous ones.These elegant lower and upper bounds for −∆ε u incorporate every parameter in our ODE system (17).Finally, we know that the lower bound in (55) is associated with the homogeneous u − v + state, whereas the upper bound is associated with the homogeneous u + v − state.Therefore, we expect a monotonous increase in the number of u + v − cells as the energy difference −∆ε u increases.
−∆ε u The states we are aiming for are highlighted with higher opacity.Nodes and their corresponding number on the axes reference the relevant equation for the transition from one state to another.
4 Tissue organization

Cell graph
In our context, cells are represented by 2D/3D points in space with a fixed radius which is equal for all cells.The Delaunay cell graph provides a reliable indication of the neighborhood relationships of the cells [29].Therefore, we initialize our graph G using the Delaunay triangulation.If the Euclidean distance between two cells exceeds the sum of their two radii, then the edge is removed from G, i.e.only cells in direct contact with each other are connected via an edge in G (Fig. 4).Edge weights are collectively set to 1.We then define the cell distance d ij as the length of the shortest path between cells i and j.

Pair correlation function
Cell differentiation patterns in our case are the result of two different cell types arising in a tissue.Patterns with the same condition have already been quantified using pair correlation functions (PCFs) [30].We use a similar approach to quantify our patterns with a PCF depending on the cell distances d ij .This requires counting different types of cell pairings for certain distances.Therefore, we introduce the sets: The pairings of all u + v − cells with distance k, S u k , are related to all possible pairings of the same distance S k by forming their ratios.Analogously, we perform the routine for u − v + cell pairings S v k to get These ratios alone will not suffice to compare the patterns for varying cell type proportions.Therefore, we normalize these by the probabilities of randomly picking two equal types of cells using the total number of u + v − cells T u and u Combined, the PCFs measure the ratios of u + v − or u − v + cell pairs within every possible distance normalized by the probability of finding these cell pairs, i.e.
For a uniformly distributed amount of u + v − or u − v + cells, the correlation function returns a value close to 1 for every cell distance k.Consequently, deviations from 1 yield information about how much more or fewer equal cell pairs are found in certain ranges.

Numerical results
In this section, we present the numerical solutions of ( 17).The explicit Euler method is used to solve the ODE until a steady state is reached.We consider two different types of signaling.Paracrine signals that exhibit low diffusivity can be described by a nearest neighbor signal.For larger diffusivities, the signal disperses throughout the tissue such that its intensity decreases with the distance traveled.

Signal construction
A signal that is secreted by one cell and diffuses slowly throughout the tissue will likely end up only affecting neighboring cells.We investigate a signal that gets activated by u Here, we used the notation N G (i) from graph theory to denote the neighbors of vertex i in the graph G.
We can also write the whole signal in terms of an adjacency matrix.For (67) this matrix will be The signal can ultimately be written as s = Au.From the steady state ( 23) we know u i = 0 for some of the cells in a heterogeneous tissue.Therefore, the minimum of the signal will also be 0. The non zero steady state has a rough upper bound Therefore, the maximum signal also obeys Using parameter combinations, such that increase the binding probability for v, tipping its fate towards u − v + .Analogously, s < s * will lead to u + v − .By definition, s i is the mean of a cells neighboring u j values.Assuming the neighbors to be in steady state and using the steady state approximation u i ≈ r u /γ u , the signal can be written as a fraction where l i denotes the number of u + v − cells adjacent to cell i.From this, we can determine the maximum number of u + v − cells in a neighborhood for the cell to still adopt the fate u + v − .Therefore, we replace s i with s * and solving the equation for l i to find Since we are looking for a natural number, the final result is Here, x describes the floor function, i.e. the nearest lower integer of a number x. Small differences between energy coefficients η u = η v + δ with δ > 0 being small, will lead to l max = 0. Therefore, a single u + v − cell will have no neighbor of equal type.At the same time, cells without any received signal, i.e. s i = 0 will adopt u + v − fate (Fig. 2).In conclusion, a cell surrounded only by u − v + cells will adopt u + v − fate, whereas a cell with a single u + v − in its neighborhood has to adopt u − v + fate.On an ideal hexagonal grid, i.e. each cell has exactly six neighbors, an ideal arrangement would amount to 1/3 of the cells being u + v − .This estimate nearly fits the simulated proportion jumps of 27% at both ends.An exact number cannot be determined, as the number of neighbors varies from cell to cell with an average of 5.5 ± 1000000 neighbors.Further increases of η u only lead to discrete increases of l max , explaining the different jumps in cell type proportions.

Signal construction
Depending on how a signal disperses in space, not only directly neighboring cells can have an impact on a cell's fate.It is possible, that the collective effect of cells that are further away might also influence its fate decision.Again, the secreted signal of a cell is activated by u i .We define the received signal s i as the weighted sum of secreted signals over all other cells Here, we use the distances d ij from our cell graph.The weights q dij −1 define the fraction of the signal that gets transported from cell to cell.Let e.g.q = 0.1, then second nearest neighbors of a cell receive only 10% of the signal of the direct neighbors (Fig. 7).The denominator in (79) is used for normalization.It describes the weights of the cell that gets the highest possible signaling weights.In a perfectly arranged circular tissue, this would be the cell right in its center due to the mean of cell distances d ij being lower.The dispersion parameter q enables us to describe the transition from a direct neighbor signal to an equally dispersed signal.For q = 0, the weights become Hence, the weights for all cells that are not directly in contact with the respective cell are 0 and we obtain a mechanism similar to the local signal (67).Alternatively, q = 1 yields This describes the case of every cell having the same impact on other cells independent of the distance between them.In summary, there is a continuous transition from a next neighbor signal at q = 0, through a distance-based global signal for q ∈ [0, 1], to an evenly distributed signal at q = 1.In matrix representation we get with the normalization factor For the estimation of the stability interval, we again use the upper bound u i < r u /γ u , such that At this point, we realize that the estimation follows the exact same procedure as before, leading to (73).
u v s q signal influence 0.1 90% 9% 1% 0% 0% 0% 0.5 51% 25% 13% 6% 3% 2% 0.9 21% 19% 17% 16% 15% 13% Additionally, v gets activated by an extracellular signal, whereas u is inhibited by the same.The signal received by the first cell on the left of the line is the sum of all cell-cell communication between one cell and any other cell in the system.The table highlights how much each cell contributes to the received signal for different dispersions q ∈ {0.1, 0.5, 0.9}.Percentages are rounded to the nearest integer.

Pattern formation
We want to investigate the effect of the distance-based signal on the formation of the patterns.Therefore, we showcase nine simulation results of organoids with different cell type proportions and different signal dispersions (Fig 8).The patterns generated for q = 0.1 can mostly be considered of the checkerboard type.In contrast to the averaged nearest neighbor signal, the signal in this case is not averaged over the number of neighbors.Cells at the boundary typically have three to four neighboring cells, whereas cells in the bulk area have a mean of six neighbors.Therefore, cells at the boundary will potentially not be able to get the same amount of signal as cells in the bulk area.The received signal however, is the deciding factor with regard to the cell fate decision in our model.The low amounts of signal received at the boundary make them more likely to adopt the u + v − fate.
As q increases, we see a higher accumulation of u + v − cells near the boundary with a slight clustering behavior in the bulk.For q = 0.9, the signal disperses strongly enough to generate an engulfing pattern, where u − v + cells are completely surrounded by u + v − cells.
The pattern formation with respect to q can be quantified using the PCFs for both u + v − and u − v + cells (Fig. 9).For comparison, we used a bisection on the stability interval to find values for −∆ε u that lead to a ratio of 89 : 88 u + v − and u − v + cells for every single q.We discover that an increase in q leads to a decrease of ρ v for large distances, i.e. less and less pairs of u − v + cells pairing in the boundary regions.Simultaneously, it increases for small distances due to the cells accumulating in the center.For ρ u , we see a slight increase for large q for small distances and a tremendous one for large distances for all q.
The slight increase at small distances comes from the fact that the u + v − cells arrange in layers at the boundary.The values for intermediate distances slightly decrease as the corresponding regions become more and more devoid of u − v + pairs.In conclusion, a distance-based signal according to (79) generates patterns ranging from checkerboard to engulfing by increasing the dispersion parameter q.Additionally, the PCFs capture the characteristics of these patterns, making it a powerful tool for pattern identification and comparison.

Cell type proportion
For different dispersion parameter values q, the proportions of u − v + show a monotonous decrease with increasing energy difference −∆ε u (Fig. 10).For low values of q, the proportions show some similarities to the local model due to individual larger jumps (Fig. 10 (a)).These jumps become less pronounced for medium (Fig. 10 (b)) and high dispersions (Fig. 10 (c)).Altogether, we have established full control over the cell type proportions.

Discussion
In this study, we have derived and analyzed a model that allows us to generate cell differentiation patterns based on a system of mutual inhibition of two transcription factors, auto-activation and cell-cell communication.The model was thoroughly analyzed and simulated patterns were characterized.

Derivation of the model from statistical mechanics
Statistical mechanics has already proven its usefulness in biological model systems like ion channel opening and closing as well as oxygen hemoglobin binding [28].These ideas have further been investigated for transcriptional regulation and were successfully applied for a wide variety of examples [26,27].To our knowledge, cell fate decision models have not been combined with statistical mechanics to date.
We derived a specific model based on two mutually inhibiting transcription factors u and v with autoactivation and an external signal inhibiting u and activating v. Assuming that auto-activation is the dominant factor in transcriptional regulation, we assume that RNA polymerase binding corresponds to the binding of u and v, respectively.Based on this, we were able to derive binding probabilities of RNA polymerase to the respective promoter.A system of ordinary differential equations was generated by combining these probabilities with constant production rates and exponential decay.As long as the auto-activation remains unchanged, minor changes in the GRN such as the removal of either the signal activation or the signal inhibition can still be managed by adjusting the equations accordingly.

Analysis of the model allows accurate determination of the stability of heterogeneous steady states
On the single cell level, we identified that the received signal determines the fate of a cell.There is a critical value of this signal, such that the cell will adopt u − v + fate if this value is undercut, or u + v − −∆ε u q 0.1 0.5 0.9 6.5 7 7.5     fate if the value is exceeded.This leads to the signal being the relevant factor of the switching behavior in this system.This describes a different point of view compared to systems that utilize differences in initial conditions to generate a cell fate switch [8,9].At the same time, models that incorporate a signal dependency, have not yet been analyzed in such great detail [13,14,24].Exact expressions for all possible steady states were derived.A stability analysis enabled us to identify parameter values, such that only the states corresponding to u + v − and u − v + fates are stable.Thus, we were able to limit the system to these two cell fates.On the level of multiple cells, we found analytical expressions of parameter bounds guaranteeing heterogeneous steady states.This means that within these bounds the pattern created by the system will always be a mixture of two different cell types.In conclusion, we have provided the necessary analytical tools to guarantee the generation of heterogeneous patterns of two different cell types.

Averaged nearest neighbor signaling leads to checkerboard patterns
In some biological systems, cell communication is hypothesized to be limited to direct neighbors.An example for this is the lateral inhibition of the Delta Notch signaling pathway in epithelial tissue of Drosophila, which has been studied in great detail [12].A different example is found in the preimplantation development of the mouse embryo, where transcription factors NANOG and GATA6 decide the fate of cells in the inner cell mass.Computational studies have investigated the effects of activation by an external signal in this biological system [31,13,14,24].A great common feature in all of these systems is the formation of checkerboard patterns, i.e. patterns in which cells of one type minimize the number of equal neighbors.Fittingly, we also found this type of pattern in our simulations using an averaged nearest neighbor signaling.We took this one step further and analyzed the possible cell type proportions one can create using this model.An analytical expression for the maximum number of equal cell types in a cell's neighborhood tells us that the cell type proportions are highly linked to the average number of neighboring cells in the system.In our 2D simulations cell type proportions below 30% and above 70% are not possible.

Distance-based signaling enables a range of patterns from checkerboard to engulfing
In addition to the nearest neighbor signal, we investigated the effects of a signal that is capable of being dispersed throughout the tissue.This global cell-cell communication enables a range of patterns.From two cell types in a checkerboard like arrangement to one cell type engulfing the other depending on the signal dispersion.The introduced dispersion parameter q allows us to artificially vary between a signal that only reaches the neighboring cells and a signal that spreads evenly in the tissue.Simulations have shown that for low signal dispersion u + v − and u − v + cells tend to avoid being adjacent to the same cell type, hence we again recovered the checkerboard pattern.Furthermore, when increasing the signal dispersion, u + v − cells accumulate more at the boundary such that overall larger clusters of equal cell types are formed.High signal dispersion leads to an ideal segregation of cells with u + v − engulfing u − v + cells.Engulfing patterns are often believed to be the result of differential adhesion of two cell types.Indeed, it has already been demonstrated that the minimization of the energy as a function of differential adhesion leads to this type of engulfing [20].Not only have we found an alternative way to generate these patterns, but at the same time we were able to unify the formation of both checkerboard and engulfing patterns under the notion of differently dispersing signals.

Conclusion
We have provided a new model to describe transcriptional regulation for a system of mutually exclusive transcription factors.Furthermore, the model was analyzed in great detail with respect to parameters and stability.The model was extended by signaling mechanisms describing the cell-cell communication.
The local and global signaling obey a simple mathematical rule depending on the number of cells it has to travel across in order to reach its destination.A detailed description of the signaling transport mechanism, possibly including diffusion and advection mechanisms, provides room for further research.Additionally, signal production and uptake of cells play a crucial role in how effective different means of signal transport might be.Another perspective can be achieved by incorporating cell growth and cell division into the model and analyzing their effect on the resulting patterns.With this in mind, our study paves the way for numerous subsequent studies regarding signal-based pattern formation in developmental systems.

Figure 1 :
Figure 1: Illustration of the GRN considered in this study.Inside the cell u and v inhibit each other.In addition to that, they activate themselves.The signal s is an external factor activating v.

Figure 2 :
Figure 2: Streamline phase portraits of ODE system (17) for a single cell and three different values for s i .Arrows show the path from the initial condition towards the respective steady states (23), (27) and(24).

Figure 3 :
Figure 3: Illustration of the different steady states at the single cell level (left) and the tissue level (right).The states we are aiming for are highlighted with higher opacity.Nodes and their corresponding number on the axes reference the relevant equation for the transition from one state to another.

Figure 4 :
Figure 4: Visualization of an tissue with 177 cells (a) and its corresponding cell graph (b).Black lines represent the cell membranes.The cell centroids are shown as black dots in both pictures.Red lines represent the edges, which provide information about which cells are in contact with each other.

Figure 6 :
Figure 6: Simulated cell type proportions for 20 equidistant values of −∆ε u spanning over the stability interval (55).Cell proportions for u + v − are colored in cyan, u − v + in magenta.

Figure 7 :
Figure 7: Illustration of the GRN represented by our model as well as an exemplary representation of the signaling in a one-dimensional cell line.Inside the cell, u and v mutually inhibit each other.Additionally, v gets activated by an extracellular signal, whereas u is inhibited by the same.The signal received by the first cell on the left of the line is the sum of all cell-cell communication between one cell and any other cell in the system.The table highlights how much each cell contributes to the received signal for different dispersions q ∈ {0.1, 0.5, 0.9}.Percentages are rounded to the nearest integer.

Figure 8 :
Figure 8: Different patterns generated by the model on a tissue geometry with 177 cells.Colors depict the values of v i in steady state.High values of v i correspond to low values in u i and vice-versa, i.e. cyan and magenta represent u + v − and u − v + cells, respectively.From left to right, −∆ε u increases.From bottom to top, the dispersion q increases.

Figure 9 :
Figure 9: PCFs for u + v − cells (left) and u − v + cells (right) for different dispersion parameters q.Any PCF represents a tissue with a ratio of u + v − :u − v + = 88 : 89.The dashed black line at 1 resembles the PCF values of an ideal uniform distribution of two different cell types.If values lie above 1, this means there are more pairs found at that distance.Consequently, values below 1 resemble fewer pairs.

Figure 10 :
Figure 10: Simulated cell type proportions with respect to −∆ε u .Simulations were performed by dividing the stability interval for −∆ε u into 20 equidistant values.Dispersion parameter q increases from (a) to (c) resulting in different scenarios.