A computational model of stem cells’ internal mechanism to recapitulate spatial patterning and maintain the self-organized pattern in the homeostasis state

The complex functioning of multi-cellular tissue development relies on proper cell production rates to replace dead or differentiated specialized cells. Stem cells are critical for tissue development and maintenance, as they produce specialized cells to meet the tissues’ demands. In this study, we propose a computational model to investigate the stem cell’s mechanism, which generates the appropriate proportion of specialized cells, and distributes them to their correct position to form and maintain the organized structure in the population through intercellular reactions. Our computational model focuses on early development, where the populations overall behavior is determined by stem cells and signaling molecules. The model does not include complicated factors such as movement of specialized cells or outside signaling sources. The results indicate that in our model, the stem cells can organize the population into a desired spatial pattern, which demonstrates their ability to self-organize as long as the corresponding leading signal is present. We also investigate the impact of stochasticity, which provides desired non-genetic diversity; however, it can also break the proper boundaries of the desired spatial pattern. We further examine the role of the death rate in maintaining the system’s steady state. Overall, our study sheds light on the strategies employed by stem cells to organize specialized cells and maintain proper functionality. Our findings provide insight into the complex mechanisms involved in tissue development and maintenance, which could lead to new approaches in regenerative medicine and tissue engineering.

role in preserving the tissue's functionality.To address the first goal of regulating non-genetic diversity, we acknowledge the crucial role of controlled stochasticity in generating population heterogeneity [3][4][5][6][7][8][9] .Therefore, our model systematically integrates stochasticity at every level, aiming to faithfully capture the intricate dynamics observed in living tissues.Besides, regulatory networks have been studied as the main decision-makers in a wide range of biological systems [10][11][12][13][14][15][16][17][18][19][20][21][22] and in the presence of stochasticity 4,6,11,17,[23][24][25][26][27] .Therefore, a multi-stable regulatory network in the form of a set of ordinary differential equations (ODEs) is considered as the decision-making unit for stem cells in our model 2,3,20,28 .Traditionally, the Gillespie algorithm is considered to be a "golden standard" for describing the behavior of systems with a small number of determinants driven by inherent fluctuations 3,6,29 without having to deal with complex mathematical equations.By considering probabilities as an integral part of any living system, here we used the Gillespie algorithm to simulate the time evolution of our designed stochastic system.
The regulatory networks orchestrate proliferation/differentiation balance to provide population heterogeneity, and maintain a homeostasis state 3 .To address the second and third objectives, we acknowledge that intercellular communication stands out as the primary factor in establishing and maintaining the structure of the population 9,30 .Consequently, our model is enriched with an intercellular communication mechanism.This means that the decision-making process of stem cells is influenced not only by intrinsic factors but also by extracellular signals capable of diffusing among population cells 3,31 .To achieve this, a reaction-diffusion process is incorporated into the model, facilitating the formation of a spatial pattern within the population and its subsequent maintenance in a steady state.
Focusing on the impact of stochasticity, the proposed model is defined based on six basic principles discussed in previous studies 2,3,9 as follows: (i) stochasticity is an inevitable part of any living cell 2,8,9,26,[32][33][34][35][36] .(ii) two major sources of stochasticity, the non-deterministic position of the cell division plane and nonuniform distribution of determinants in the cell lead to the random distribution of the cytoplasmic molecules among daughter cells during cell division 2,9,[37][38][39][40][41][42][43] .(iii) cell fate is determined based on the number of determinants in the offspring upon the completion of cell division and assumed to be fixed during cell life cycle 2,9 .(iv) cell determinants interact with each other via an internal switch 44,45 .(v) the decision bias in the internal switch is determined by model parameters representing interactions between the switch elements 2,9 .(vi) the switch parameters could also be affected by the cell location in its environment, and it is the key to the spatial pattern in the population 9 .
In this project, we propose a computational model to describe the stem cells' internal mechanism leading to self-organization by considering signal diffusion in cell-cell communication, and stochasticity at all levels of simulations.This project demonstrated the ability of our designed model to function as a stem cell mechanism underlying the formation and maintenance of desired spatial patterns in the population (fitted to population functioning).Besides, in the presence of controlled noise our model could easily reach the required population heterogeneity and homeostasis state.Our previous model was modified since it was able to manipulate cellular decision-making biases in order to allow daughter cells to be born in their destined territory, then defend it.Furthermore, we show that the desired spatial order of cells in the population could be determined by a "leading" initial signal corresponding to that desired pattern.In addition, we discuss the impact of the death rate as a key factor in our system.Finally, to further illustrate the strength of the model, we explore the behavior of the system in the absence of one of its constraints, the dish wall.

Mathematical model of the system
We consider building on our previously designed model 3 to study the capacity of that to imitate the behavior of human embryonic stem cells (hESCs) in early development.To be specific, we aim to model the production and rearrangement process through which a mass of homogeneous cells, as the initial state of the system, gives rise to a desired spatially ordered sequence 46 as the final state of the system.A spatially ordered population of cells determines the final state of a system, which consists of stem cells (S), two specialized cells (A, and B) and progenitor cells (P) as intermediate cells.It is worth mentioning that the presence of progenitor cells prevents the accumulation of mutations in the stem cells, by taking the responsibility of producing specialized cells through a large number of cell division cycles 3 .The stem cells are capable of self-renewing and/or differentiating into progenitor cells, and they appear to be a bi-stable system.Progenitor cells, on the other hand, can self-renew and/or give rise to two non-dividing terminally-differentiated cells termed A, and B 3,46 .Like the previous version of the model, the progenitor cells are considered as a tri-stable dynamical system with a limited capacity of proliferation and restricted potential of differentiation 47 .To be able to study the potential reactions in the model we should emphasize the existence of three phases in the system: the first phase which is demonstrated by a homogeneous mass of stem cells, the second phase which starts with the first division of stem cells and ends with a population of stem cells, progenitor cells, and differentiated cells in a desired spatially ordered pattern, and the third phase which begins with pattern formation, as the final state of the second phase, and terminates to the tissue homeostasis state.
Pattern formation could only happen in response to the signaling pathways 30 .As BMP4 breaks symmetry in inner cell mass during development, our model also requires a signal, referred to as the "leading" signal, which performs the same function.As we prefer to challenge the model to be able to control the overall population behavior without any outside controller, we assume that in the first, and second phases of the simulation, population cells secrete the primary signal needed to form the desired ordered array of cells in the final state.We are aware of the fact that in early development, pattern formation is much more complex and occurs along a hierarchical path.A first signal establishes the foundation for initial symmetry breaking.This leads to the birth of distinct differentiated cells, and then the presence of distinct signaling pathways.These pathways also could lead to producing even more differentiated cells by breaking the symmetry in their corresponding territory.Finally in three weeks the complex structure of Grastula emerges as a result of a finite loop of secreting a signal and breaking the existing symmetry 30,46 .However, in this model, we skip the hierarchical aspect, and the final pattern is formed in one step.A model that can establish and maintain a desired population pattern in two phases could be utilized by the new-born cells in the population, in order to facilitate the above-mentioned hierarchical process.
In the first phase, stem cells do not divide.However, they are involved in the process of producing the required "leading" signal.When the division of the population is triggered in the second phase of the simulation, progenitor cells (P) emerge from the population.Subsequently, differentiated cells (A and B) are born and finally formed in a spatial pattern corresponding to the "leading" signal.In the third phase, the system reaches its stable steady state and maintains it.
Because of the remarkable degree of flexibility that reversible transfer of cells could provide in the case of injury, and even under normal conditions, here, we believe that stem cells can transit to progenitor cells at a fixed rate and vice versa 3,[48][49][50][51][52][53][54][55][56][57][58] .Within the described system, the dynamics of the model are described as follows: where η , η s , and η p represent the rates of three different division types of stem cells, ω p/s denote the transition rates of S/P and p , A , B , µ d , µ A , and µ B denote the rates of six different division types that progenitor cells can go through.In addition, γ A and γ B indicate at which rate A and B cells diminish from the population.Consider- ing n S , n P , n A , and n B , as the average densities of cell types S, P, A, and B respectively ( n S,P,B,A are cell numbers normalized by volume), the time evolution of n S,P,B,A is given by It was previously proven that this system can reach its stable steady state on (n * A , n * B , n * P , n * ) , as long as the fol- lowing conditions are satisfied 3 : where n is the average total density of cells, computed as n = n S + n P + n A + n B , and η S = η S (n) .We set the parameters in the model in a way that these two essential conditions are satisfied.

Stem cells' internal mechanism
The following set of ODEs used in several projects 2,20,28 to describe a two-element bi-stable regulatory network, here is employed as the stem cells' internal mechanism: It actually provides the regulatory mechanism which takes care of the part of the model dynamics in Eq. ( 1).The form of the above-mentioned ODEs is stand on some the fundamental assumptions of the model: first, X s and Y s are the relative amount of two cell determinants, and as it is deducible from the phrase that their values determine the cell's final fate, second, the mutual repression effect of the determinants (in the form of Hill functions), and their degradation rate describe the dynamical behavior of SCs.Here, n, β s , ι X s , ι Y s , and γ are set as the Hill coefficient, the effective rate of determinants synthesis, inhibition rates of X s and Y s , and the degradation rate, respectively.
In general, the parameters are tuned (by try and error method) to build the SCs' regulatory switch with two stable steady states corresponding to two main cell fates in the model, stem cell type (S) and progenitor cell type (P).In addition, to have the first condition in Eq. ( 6) be satisfied, as ι x s is the parameters which directly controls the rate of symmetric division of S −→ η s S + S , we need to set ι x s as a function of n, where ι x s (n) ′ < 0 3 .

Progenitor cells' internal mechanism
The regulatory mechanism of the progenitor cells is described by the following set of ordinary differential equations (ODEs) as a two-element tristable system 2,20,28 :  (3) Likewise, it is assumed that X p and Y p are two cytoplasmic determinants whose values determine the final fate of progenitor cells' offspring 3 .Here in this system, n, α x p /y p , β p , ι X s /Y s , and γ are studied as the Hill coefficient, activation rates, the effective rate of determinants synthesis, inhibition, and degradation rate, respectively.There are two more parameters, ε s 1 , and ε s 2 as "signaling effect coefficients" that will be discussed in this section.The parameters are set to build a tristable steady-state system with three fixed points corresponding to three cell fates, one progenitor (P), and two differentiated cell types (A and B).The phase plane of this tristable steady-state system represents the domains of the three cell fates.These domains indicate the decision boundaries governing the daughter cells' final fate right after their birth 3 .As discussed before, in the first and second phases of the simulation, a "leading" signal, coined as S l , is pro- duced to organize differentiated cells in order to form the desired pattern in the population.Besides, to maintain the population pattern in steady states, it is essential for each differentiated cell, A/B, to secrete their distinct signaling molecules, S 1 , S 2 to conquer their territory 3 .As shown in Eq. ( 8), the impact of the signals S l , S 1 , and S 2 is performed on progenitor cells' internal switch via α x p , and α y p through ε s 1 andε s 2 parameters.The former's impact arranges the population's cells, while the latter maintains/fixes them in their territory.
The "signalling effect coefficient", ) , where i ∈ 1, 2 .The model primarily relies on S l to estab- lish a spatial pattern in the population.This implies that it is unnecessary to retain it in the model once the pattern has been established.Therefore, here, f (S l ) = a S 12 l 1+S 12 l , in the absence of S i , and f (S l ) = 0 , otherwise.Furthermore, it is biologically reliable, since that is observed during development, and it minimizes the energy costs of the system.g(S i ) = aS i b , if S i ≤ b , and g(S i ) = a , otherwise 3 .Finally, to have the second condition in Eq. ( 6) be satis- fied, the parameter β p is tuned (by try and error method) in such a way that the rate of symmetric division, P P −→ P + P is much less than divisions in the forms of P

Signalling dynamics
Here, the "leading" signal, S l , has been studied in three ways.A fixed deterministic pattern is first considered.Second, the dynamic of the S l signal is described as follows: where D, α s l , and k are defined as diffusion, production and degradation rates.Third, to generate S l signal, a Turing system is used as follows 59 : It is the general model of a Turing system that promises to generate spatially heterogeneous patterns with two diffusive chemicals ( s l and s a ) from uniformly distributed sources.Here, S l demonstrates the number of the "leading" signal, and S a is the number of its associated substance which differs in diffusivity.The parameters, d l/a , and γ s l are diffusion coefficients and a common factor multiplied by the reaction terms.In Eq. ( 10), it is assumed that f (s l , s a ) = A s l − s a + C , and g(s l , s a ) = B s l − s a − 1 as the linear reaction terms, where A, B, and C are constants.To observe the formation of stable patterns in a system with linear reaction terms, it is necessary to constrain the S l value within a finite range 59 .To be specific, it is assumed that s lower ≤ S l ≤ s upper , where s lower s upper are constants.Since the presence of stochasticity is inevitable in our model, the Gillespie algorithm is used to create Turing patterns.Therefore, the model's parameters had to be tuned (by try and error method), compared to the ones set in Shoji's model 59 .Figure 1 shows the effect of parameters, s upper , and d u/a on the Turing patterns gen- erated by Gillespie's algorithm.In Fig. 1, from right to left, (d l , d a ) = (5e3, 1e5), (1e4, 2e5), (2.5e5, 5e6), (1e6, 2e7) , and from top to bottom, s upper = 5, 10, 30 ( s lower = 0 ).Patterns in the first, second, and third rows of Fig. 1 coined as reversed-spot, stripe, and spot patterns respectively 59 .A set of reaction-diffusion equations describe the dynamics of the signaling molecules, S 1 , and S 2 as follows 3 : Here D, k, α s i ( i ∈ 1, 2 ) represent the diffusion coefficient, the decay, and the production rate of the signaling molecules, respectively.Besides, n is the Hill coefficient, and β is the effective rate of signaling molecules' synthe- sis.The parameters are set to describe the dynamics of signaling molecules, S 1 , and S 2 , in the form of a bi-stable system 3 .In this model, for generating reversed-spot, stripe and spot signalling patterns, (d l , d a ) is set equal to (2.5e5, 5e6), (1e4, 2e5), (1e4, 2e5) and s upper = 5, 10, 30 ( s lower = 0 ), respectively.( 8)  11)), and second, by increasing the birth rate of A (B) cells through the ε s 1 ( ε s 2 ) in Eq. ( 8).

Gillespie algorithm
The Gillespie algorithm is used to capture the time evolution of the system 2,3,9,29,60 .Table 1 illustrates 28 potential reactions that can occur in each time step, as well as their corresponding propensity functions and parameters' values.As the introduced propensity functions represent high-order reactions, they are used as an approximation to the Gillespie algorithm 61 .The endpoint of the simulation is when the stem cells have gone through 50 divisions on average.
The simulation starts with an initial population of stem cells distributed randomly in a hypothetical dish.Each individual mesh of the dish can be potentially occupied with one of the four cell types, S, P, A, or B with or without the presence of signals S l , S 1 or S 2 .The number of determinants is chosen randomly from the stem cells' territory in the phase plane describing its dymanics 3 .In the first phase of the simulation, a mesh occupied with a stem cell can be the scene for eight (11 in the third scenario) potential reactions: production/degradation of determinants X S , and Y S , production/degradation/diffusion of S l signal (plus those of S a signal in the third scenario), and stem cells' movement.As shown in rows 1 to 4 in Table 1, the propensity function of production/degradation reaction of determinant X S ( Y S ) is determined based on the positive/negative term in the first (second) ODE in Eq. (7).The movement is provided in the model to prevent the system from being blocked (see Khorasani, and Sadeghi 3 for details).As it is depicted in rows 7, and 14, the movement can occur at a constant rate ( m s , and m p for stem cells and progenitor cells, respectively), and this involves randomly selecting an empty mesh, and then placing the cell from the current mesh into that empty space.Rows 23 to 25 illustrate the propensity functions of production/degradation/diffusion of S l signal in the second and third scenario captured by Eqs.(9), and (10), respectively.In the same manner, rows 26 to 28 show the propensity functions of production/ degradation/diffusion of S a signal in the third scenario.At each time step of the Gillespie algorithm, one of the above-mentioned reactions is chosen and based on that in the current mesh, one of the determinants X S , and Y S , or signal S l ( S a in the third scenario) decreases/increases in value or in one of the neighboring meshes, the value of signal S l ( S a in the third scenario) decreases/increases.It is necessary to mention that all these reactions also occur in the second and third phases of the simulation.However, they are the only ones with a chance to happen in the first phase.The second phase of the simulation starts when the S l signal is organized in the desired pattern.In this model, it is triggered after a certain number of iterations from the start point.As we move into the second phase, two reactions with prominent roles are revealed: the division of SCs and the transformation of the SCs with constant rates of r s , and w p , respectively.With these two reactions, progenitor cells and as a result, specialized cells (A and B cells) and their corresponding signal molecules emerge in the population.In other words, all the reactions listed in Table 1 can play a role in the system's scene from now on.
When a stem cell divides into a progenitor cell, the number of determinants X P and Y P in a new-born cell is initiated randomly from the middle domain of it's corresponding phase plane 3 .The propensity functions of production/degradation reaction of these two determinants are captured based on Eq. 8 (see rows 8 to 11) 3 .A progenitor cell can be divided with a constant rate of r p (row 12) to P, A, and B cells.With A and B cells in the population, S 1 and S 2 signaling molecules can set foot in the system, and in a mesh occupied with S 1 , and S 2 , five reactions, production/degradation of S 1 , and S 2 , and diffusion got the chance to happen.Rows 17 to 22 demonstrate their corresponding propensity functions based on Eq. (11) (see [Khorasani, and Sadeghi 2022] 3 for details).
For the rest of the reactions, the propensity functions are chosen as constant rates 3 .In the model, death is defined as the omission of the chosen cell from the population.During the division, two offsprings are located in the current mesh and one of the empty neighboring meshes.The number of determinants X 1 , and X 2 ( Y 1 , and Y 2 ) in two new-born cells are set as , where #X ( #Y ) is the number of determinant in the mother cell 3 .The transition from a stem cell to a progenitor cell occurs with the constant rates of w p , where during that the S cell in the current mesh is replaced with a P cell.The number of determinants X P and Y P in transformed cell is initiated randomly from its specific domain of it's corresponding phase plane 3 .The transition from a P cell to a S cell occurs with the constant rates of w s with the same manner.

Scoring algorithm
To evaluate the maintenance of the spatial patterns in the population, we used the "scoring algorithm" 3 .To implement this algorithm, we need to construct two filter matrices corresponding to the "leading" signaling pattern, namely template and penalty.For each initial "leading" signaling matrix, if a mesh value is greater than 250/2 (250 is the maximum value of leading signaling molecules in each mesh) the corresponding mesh's values of the template matrix are equal to 1, and it is equal to −1 , otherwise.The values of the outer-dish meshes are set to 0. The values of the penalty matrix are equal to 1 in out-of-dish meshes and are equal to 0 otherwise.Each simulation starts with an initial state with a desired spatial pattern and continues until all progenitor cells have undergone an average of 50(100) divisions.Out of all the states encountered during the simulation, we select 500(1000) snapshots with dimensions of d × d .For each snapshot of the system state a corresponding 2d × 2d -dimensional matrix is created as follows: In each snapshot, meshes corresponding to yellow/red pixels (A/B  To measure the score values for each snapshot, we slide the template, and penalty filters over each of the 500(1000) 2d × 2d-dimensional matrices, from top-left to bottom-right (mesh by mesh), multiply their cor- responding values one by one and compute the summations.Each time the filters are moved over the snapshot matrix (one mesh at a time), two values are computed-one for the template filter and another for the penalty filter.The subtraction of these values is calculated, and the resulting score assigned to each snapshot is determined as the maximum among all subtracted values obtained by filters' sliding over the snapshot's matrix.Finally, we normalize the score values in the range of [0,1] by dividing to the maximum score value (where the system state perfectly matches the template and penalty matrices).

No "leading" signal
The computational model of SCs' mechanism designed by Khorasani et al. 3 can maintain the initial spatial pattern in the population.What if there is no initial pattern, and the initial state of the system is a population of randomly distributed stem cells?Two potential scenarios could happen: first, the four types of cells, S, P, A, and B, are randomly distributed in the population, second, in some parts of the dish, P cells divide into more A/B cells, certainly by chance, and then they maintain their territory by secreting S 1 /S 2 signaling molecules to the end of system time.To address this question, the first simulation starts with a population of stem cells in the dish as shown in Fig. 2A, in which the S cells are shown in cyan and the blue and gray pixels represents the empty and out of dish positions as indicated in the colorbar.By implementing Gilespie algorithm and the reactions described in Table 1, we let the system be updated to the point that all the P cells have been through at least 50 divisions.Figure 2B illustrates the final state of the system containing all four cell types S, P, A, and B in cyan, green, yellow and red, respectively.The results support the second scenario as we expected.The specialized cells (A, and B) are randomly born in an area and defend that as their territory.It is worth-mentioning here that the initial state of all simulations are represented the same as the population of S cells randomly distributed in the dish.

Fixed "leading" signal
At the beginning of the development process, in the presence of BMP4, hESCs are manipulated to produce an ordered array of specialized cells.Finally, they generate spatial patterns in the population 30 .It is evidently concluded that the presence of a "leading" signal, imitating BMP4's role, is essential for the formation of patterns in the population.Here, the second experiment is designed to demonstrate that starting with a population of stem cells, our designed system has the potential to generate an ordered array of differentiated cells in response to an initial "leading" signal.In the second experiment, fixed patterns of S l signals are assumed as the "leading" signal as shown in the first and third columns of Fig. 3 in which the "Red" and "Magenta" colors represent the maximum and zero values of signalling molecule in the corresponding mesh.The results in the second and fourth columns show that the model can effectively self-organize in response to the "leading" signal.Spot and reversedspot patterns could not be perfectly formed in the population.Although it is expected 3,31 , since the territory of one side is dominant 3,31 , which makes it easier for dominant cells to invade the apposite side.

Dynamic "leading" signal
In the third experiment, we challenge our model with dynamic, biologically reliable "leading" signals.Besides, in this experiment, the "leading" signal does not last to the end point of the simulation but for a limited duration to initiate pattern formation in the population (at the end of phase #2 ).It is in agreement with biological experiments done previously to recapitulate spatial patterning in early development 30 .For this purpose, here we used four dynamic patterns, Gaussian, Spot, reversed-Spot, and Stripe, as the "leading" signal, shown in the first row of Fig. 4 59 in which the colorbar indicates the number of signalling molecule is each grid.
To study the capability of the model in initiating and maintaining the desired spatial pattern in the population and in the presence of dynamic "leading" signaling molecules, the medium is populated with randomly distributed SCs.The dish radius is initially set to 100 in the sense that any cell would be restricted to a maximum distance of 100 pixels the center point.To compute the time evolution of the cell populations, a stochastic simulation using the Gillespie algorithm is applied.Simulation is terminated when the progenitor cells have gone through 50 divisions on average.A total of four simulations were run corresponding to four different types of "S" signals, Gaussian, Spot, reversed-spot, and Stripe.In the first phase of each simulation, SCs do not divide.Though, they are involved in producing the required"leading" signal based on Eq. ( 9) for Gaussian pattern and Eq. ( 10) for Turing patterns.The first row in Fig. 4 demonstrates the state of the "leading" signal, S l , at the end of phase #1 .Until the end of phase #1 , the state of SCs in the dish remains unchanged and the same as the dish state in Fig. 2A.In the second phase, it is different.By triggering division in the second phase of the simulation, progenitor cells, and subsequently differentiated cells, A and B, emerge.A and B cells in our model are self-organized to generate the desired spatial pattern in the population.The internal mechanism of P cells is designed to achieve this (Eq.8).Based on Eq. ( 8), in the parts of the dish with maximum/minimum value of "leading" S l signal, the number of X P /Y P determinants exceeds the number of Y P /X P ones in P cells, and as a result with a great probability, A/B cells are born in these domains.Therefore, a pattern corresponding to the initial S l signal is generated in the population.The state of the system at which all progenitor cells have been through at least 50 divisions is shown in the second row of Fig. 4. In all four simulations, specialized cells, A, and B properly self-organize to create the desired pattern.The subplots in the third row of Fig. 4 show the diagrams of the abundance of four types of cells from the beginning to the end of the simulation.All diagrams have saturated so fast, and it indicates that the system can easily reach its stable, steady state on (n * A , n * B , n * P , n * ) as proved and promised in Section "Mathematical model of the system".

Death rate
Here, we study the death rate values of A and B cells ( γ A and γ B ).In this model, we have tuned more than 40 parameters (by try and error method), where all of them except one ( γ A/B ) orchestrate the production rates in the population.The balance between production and death rates determines the abundance of cells in the tissue.The parameter γ A/B is the only parameter to control cells' removal from the system.Therefore, the γ A/B value is critical in determining the n value in our system.On the other hand, it was previously proven that this system could reach its stable, steady state on (n * A , n * B , n * P , n * ) , as long as the conditions in Eq. ( 6) are satisfied 3 .The parameters in the model are set in a way that both of these essential conditions are satisfied.However, the satisfaction of these conditions is independent of the γ A/B value.In other words, different values of the γ A/B cannot affect the system's stability.Therefore, theoretically any positive value for the γ A/B is acceptable.However, the γ A/B value directly impacts the abundance of different cell types in steady state, i.e., n * A , n * B , n * P , and n * values.As to form a desired pattern, a specific num- ber of n * A , n * B , n * P , and n * is certainly needed, one could easily conclude that the formation of a desired pattern directly depends on the γ A/B value in our model (keeping in mind that the rest of the parameters, with a key role in the production rate, are fixed).On the other hand, based on the time evolution of n A , and n B in Eq. ( 5), It makes it clear that γ A/B , and n * P values are both mutually dependent.However, at the beginning of the simulation, we do not have access to the value of any of them to fix the other one.Even knowing one of these values, we cannot determine the exact value of the other one.To be specific, due to the stochastic nature of the system, it is impossible to obtain the exact values of A/B , µ d , and 2 µ A/B through simulation.
Looking on the bright side, we can approximate the n * A and n * B values based on the initial pattern of S l , the "leading" signal.Assume that in an initial pattern (the first row of Fig. 4) p percent ( (1 − p) percent) of the total positions are red/yellow (green/cyan/blue/magenta), the expected territory of A (B) cells.By utilizing trial and error method, to get a correct pattern in the populating corresponding to the initial S l pattern, the number of stem cells and progenitor cells with empty positions can not exceed 15 percent of the total dish's positions.It is instantly concluded that the desired pattern occurs in the presence of p * 0.85 * (N) (  where i ∈ A, B .It is worth mentioning here that the γ A , γ B values reported in Table 1 are the average values computed through the whole simulation.The results shown in Fig. 4 indicate that with this method of determining γ A/B value, the formation of the desired pattern as well as reaching and maintaining the stable, steady state on (n * A , n * B , n * P , n * ) are guaranteed in our model.As the stability of the system was proved previously, one can instantly conclude that if we could calculate the average value of γ A/B , accurately, based on Eq. ( 12), and substitute it for the death rate of A/B cells in the model, we would reach the same steady state of our system on (n * A , n * B , n * P , n * ) .However, the dynamics of the system would certainly differ.Aiming to find a suitable solution to determine γ A/B value, we studied two equations of According to these equations, γ A/B can be determined based on the ratio of n A/B (calculated by multiplication of n P , and production rate of A cells in the population), to n * A .To imitate the behavior of this dynamic, in the fourth simulation, we set the γ i = n i n * i * v g , where i ∈ A, B , and v g = 0.0125 (Running different simulations, we reached the v g value.Using this v g value, the system reaches almost 15 percent of total dish pixels occupied by stem cells, progenitor cells, and empty positions.The rest of the fourth experiment is designed and run similarly to the third one.Figure 5A, the first column from top to bottom, shows the states of the system after (almost) 30, and 100 divisions, and the abundance of four cell types (through time) respectively.Here, the "leading" signal is spot pattern and the γ A/B value is set based on the ratio of n A/B , and n * A/B .The result shows that although the desired pattern is formed, it is not maintained properly in the system.The narrow regions of each domain (A/B domains) are merged.The reason is that in this last method of determining γ A/B value, there is a meaningful delay between the detection of a change in the number of cells, and the corre- sponding regulation of γ A/B value.In other words, a great difference between n A/B , and n * A/B is needed to regulate the γ A/B value.Therefore, there is always this chance that the cells on one side grow in number and capture the other side's domain before the cells on the other side get the chance to defend their territory.For this reason, we modified the γ A/B equations to increase the sensitivity of the system: the experiment with this new formulation of γ A/B value for spot pattern as the "leading" signal in the population.
The results are shown in the second column of Fig. 5A where the first, and second subplots show the state of the system after (almost) 30, and 100 divisions and the third subplot demonstrates the abundance of four cell types through the simulation.The results illustrate that this modification affects the formation of more proper Spot patterns in the population compared to the former one (first column in Fig. 5A).It is also worth mentioning here that the traces of the cell types' abundance (third row in Fig. 5A) are saturated and it confirms the stability in the model which was promised in Section "Mathematical model of the system".

Maintenance of the Spatial Pattern in the population in addition to the cells' abundance
To evaluate our model's capacity for pattern maintenance by numbers, we leveraged the scoring method introduced by Khorasani et al. 3 .Figure 5B shows the template and penalty matrices for the spot pattern.For both experiments in Fig. 5A, during the simulation, we captured 500 snapshots of the model's state.Then for all the selected snapshots, we slide the template and penalty matrices on the snapshot matrices and the corresponding scores (discussed in section 1.6) are calculated.Fig. 5C indicates the scoring traces of the first and second simulations in Fig. 5A in blue and magenta.In other words, the scoring algorithm compares each snapshot with the filter chosen corresponding to the final state of the S l signal in phase #1 , and the scoring traces demonstrated in Fig. 5C indicate how much the structural pattern in the population differs from the "leading" signal's pattern through the simulation.Comparing two traces in this subplot also indicates that with the second formulation of gamma, the spatial pattern is maintained more properly in the population.
To evaluate the maintenance of the Spatial Pattern in the population with Gaussian, reversed spot and stripe signaling patterns, three more sets of simulations was run.In Fig. 6, from top to bottom, the template matrices (for the scoring algorithm), the final state of the system after almost 100 divisions, the abundance of four cell types, and the score traces through the simulation are shown.For reversed-spot and Gaussian patterns, the scoring traces are saturated (Fig. 6D,H).It can be interpreted as the proper maintenance of the established pattern in the population.The interpretation is also confirmed by the results shown in the Fig. 6B,F.The scoring trace for stripe patterns decreases rapidly after pattern formation (Fig. 6L).This behavior matches the system's state shown in the Fig. 6J.In other words, in the face of the stripe pattern, pattern formation could occur however, it cannot be maintained in the population.As shown in Fig. 5C, the system behavior in the case of spot pattern differs three other cases (after pattern formation).The spot pattern scoring diagram is not saturated like reversed-spot and Gaussian patterns, and it also does not decline rapidly like stripe patterns.It decreases so slightly.So slightly in the sense that even after 100 divisions, the loss of pattern in the population cannot be visually detected (Fig. 5A, second column).The scoring traces for the second and third experiments are also calculated and show in Supplementary Figures S1, and S2.

Self-organization in the absence of fixed border
Our goal in this section is to test the model's capabilities by studying system behaviour without a dish wall as a model's constraint.Dish walls act like a template in the simulations.Biological cells, however, do not require templates to create structures.More precisely, in this experiment, the cells are not restricted within the boundaries of the dish.In this simulation, it is assumed that the physical boundary of the medium is infinite.In a part of the medium, stem cells are planted in a circle, then the first phase of the simulation begins.
To perform this simulation, we equipped the model with a pseudo-adhesion factor.In biological experiments, adhesion is a factor to ensure both cell and cell, and cell, and tissue cohesion (through cell-cell adhesion, and cell-matrix adhesion).Cadherins are trans-membrane proteins that mediate cell-cell adhesion.Multi-protein structures mediate cell-matrix adhesion.These protein-based macromolecules are mainly produced locally by cells.With mediating attachment they can control proliferation, division, apoptosis, and travel, and through these processes facilitate tissue homeostasis.
For simplicity, we do not consider an additional factor as a multi-protein structure responsible for adhesion in the model.Rather, we assume that the signaling molecules produced in the model also plays the role of these proteins.The simulation is mainly similar to what was discussed in Section "Dynamic "leading" signal" with minor changes: • Cells can move and be born outside the hypothetical borders of the initial population, • It is assumed that initial cells in the population make an extracellular matrix, • The cell-matrix adhesion is mediated by signaling molecules, • With no cell-matrix adhesion division is inhibited, apoptosis is triggered, and migration is restricted (cells cannot easily move outside the matrix).
The result is shown in Fig. 7. Here, the "leading" signal, S l , is considered in a Gaussiann pattern.Although the outer border is crooked and a limited number of cells are observed outside the borders, the initial desired pattern in the population is clearly observed.These results show the capability of the introduced model in formation of the desired pattern, even in the absence of a template.

Discussion
Here, in this project the overall behaviour of the population is orchestrated by the internal decision-maker of stem cells, in the presence of stochasticity at the entire levels of the system, and without using any template or considering any movement process.An organism's life depends on the precise functionality of its organs.Besides, a proportionally organized structure is essential to the organs' proper functioning 3 .There is still much to be done to develop a comprehensive model that controls a stem cell's decisions in starting from one cell, and ultimately generating and maintaining a proper structure in the final developed tissue.In our last project, we took a step toward this tempting big picture and introduced a simple model that could properly maintain the pattern in the tissue.Here, we took another step to imitate the early development processes of production and rearrangement.Through these two main processes a mass of homogeneous cells gives rise to a desired spatially ordered sequence 46 as the final state of the system.In this model the initial mass of homogeneous cells is presented as a population of stem cells (S).
The model is described in Eqs.(1, 2, 3, 4), in details, and is simulated in four phases.The computation in the rest of the section 1.1 proves that in the second phase of the simulation, the described system could reach, and maintain the homeostatic state on (n * A , n * B , n * P , n * ) under some reachable conditions.Besides, Eq. ( 5) indicates that theses four values are determined by the set of parameters in the model.By beginning with a population containing source cells, and choosing the right set of parameters, it proves that any desired spatial pattern can be organized using the right ratio of stem cells, specialized cells, and progenitor cells ( n * A , n * B , n * S , n * P ).However, achieving this desired pattern in the population requires a corresponding signal pattern to break the symmetry in the population.The Turing patterns were used as "leading" signals for this purpose.We also equipped the model with a factor to receive this "leading" signal and translate it into the intended pattern.The results reflect the capacity of our model to describe the behaviour of hESCs in pattern formation during early development.
In this model, two sets of ordinary differential equations are defined to describe the internal regulatory mechanism of stem cells (Eq.7), and progenitor cells (Eq.8) 3 .Previous studies 3 reported that regulatory networks can provide a balance between proliferation and differentiation, as well as heterogeneity in populations.However, pattern formation could only happen in response to a "leading" signal 30 to break the symmetry in the initial cell mass, and direct the population toward the desired pattern.The results support this fact.The results indicate that first, with no "leading" signal a random pattern is formed in the population (Fig. 2), and second, the formation of any desired pattern in the population is possible in the presence of its corresponding "leading" signal (Fig. 3).Although the "leading" signals in Fig. 3 prove the capability of the model in pattern formation, the signalling patterns are not biologically reliable.Studying the system behaviour in Fig. 4, indicates the capability of the designed model in pattern formation, as well as reaching and maintaining the stable steady state of the system on (n * A , n * B , n * P , n * ) even in the face of stochastic, and dynamic "leading" signals.In the scene of multi-cellular organisms, the production of cells is not the only key role in the development, and maintenance of the organizes tissues 46 ..Therefore, in this project, we assigned a critical role to cell death rate as a control parameter.The first, and second rows of Fig. 5A show the corresponding states of the designed system (after 50, and 100 divisions, respectively) that approve this prominent role.Regarding the stability of the stripe pattern within our model, upon studying the results, it becomes evident that while our model is successful in pattern formation, it encounters challenges in maintaining the stripe pattern.This deficiency arises from the inherent vulnerability of narrow cell territories, where cells struggle to defend their domains, providing opportunities for rival cells to penetrate.Notably, due to the same diffusion rates of signaling molecules in the borders between the domains of cell types A and B, there is a facile interpenetration of territories, as visually depicted in the figures.In the context of the stripe pattern, wherein both cell types exhibit narrow territories, these territories essentially act as borders.Consequently, under certain conditions, A cells encroach upon the territory of B cells, and vice versa.This incursion disrupts the proper maintenance of patterns within the population.
In the third row of Figs.5A, and 6 the abundance of all four cell types from the beginning to the end of the simulation is shown (and also in the last row of Fig. 4).In all cases, the cell types diagrams are saturated, which confirms the stability in the model which was proved and promised in the section 1.1.
The last row of Figs.5A and 6 (as well as sumlementary Figures 1, and 2) shows scoring diagrams based on the borrowed algorithm introduced by Khorasani, et al. 3 .The diagrams of Gaussian, and reversed-spot patterns are saturated (after pattern formation) which indicates the capacity of our designed model in maintenance of the organized pattern in the population (in addition to pattern formation, and reaching and maintenance of demanded cell abundance).However, the diagrams of both spot and stripe patterns decline over time.As discussed earlier in this section, the territory corresponding to these patterns is not wide enough to be maintained.
It is worth mentioning here that the magenta scoring trace of the spot pattern declines slightly in a way that it cannot be visually detected in the Fig. 5C.From the stability point of view, we need to facilitate the model to a different reaction-diffusion equation to describe the S 1 , and S 2 behaviour in order to provide a great "defensive shield" for narrow territories (i.e. to provide saturated scoring diagrams).Biologically, this slight decline indicates that the pattern is not perfectly maintained but still satisfactory (in a spot pattern).Satisfactory in the sense that it maintains the pattern to the degree that the tissue can faithfully maintain its functionality.This interpretation is biologically reliable because the form of organs is constantly changing in living organisms.There is no organ in your body that remains completely similar.However, until some point it remains satifactorily the same and fulfills its functioning and not afterwards.
Lastly, a noteworthy point pertains to Fig. 7. Our model's distinctiveness lies in its ability to coordinate the entire system's behavior without relying on external control parameters.In our latest experiment, we examined the system's behavior in the absence of a dish wall.In multicellular organisms, direct interactions and the extracellular matrix are pivotal in holding cells together, fostering cohesion.This cohesion, mediated by proteins secreted in the extracellular matrix, governs the formation of organized multicellular structures.In the absence of a dish wall, our model considers the simplest form of cohesion, as detailed in Section 2.6.The outcome illustrated in Fig. 7 highlights that even without the constraint of a dish border, we can effectively regulate the overall behavior of the system.
In order to ensure the proper functionality of a tissue, it's fascinating to design a model that can form and maintain a well-organized structural pattern in the tissue.This will guarantee life.We are aware of the fact that here, we designed a simplified replica of two main processes of production, and rearrangement during early development.Pattern formation is much more complex and it occurs along a hierarchical path starting from one totipotent cell.Thus, in future work, it is essential to investigate stem cells' strategies to reach a developed, organized structure along a hierarchical path.This strategy should provide a finite number of cycles with the following instructions: first, producing the demanded abundance of cells, second, secreting a signal third, breaking the existing symmetry in the population and fourth, producing the differentiated cells as the source of the next level of development in a spatially ordered sequence.It could be a promising step toward the big picture of this project: generating organoids starting from one stem cell.

1 − ks 2 Vol
https://doi.org/10.1038/s41598-024-51386-zwww.nature.com/scientificreports/Here it is critical to address one fundamental question: How can A (B) cells maintain their territory?Promptly after the birth of A (B) cells in the population, the secreting of S1 (S2) signaling molecules is triggered.They defend the A (B) territory by first, preventing the production of S2 (S1) (based on the second term in Eq. (

cells) are to 1 /
− 1 , and the rest of the elements are set equal to 0. This d × d matrix makes the middle part of the corresponding 2d × 2d-dimensional matrix.The rest of the values are set equal to 0.

Figure 2 .
Figure 2. The initial and final state of the system in the absence of "leading" signal.Each pixel indicates one individual cell, S in cyan, P in green, (A) in yellow, and (B) in red, an empty position (EP) in blue, or an out-ofdish space (OD) in gray.

Figure 3 .
Figure 3.The presentation of static "leading" signal ( S l ) (first, and third columns) and their corresponding patterns formation (second, and fourth columns) in the population.In the leading signal patterns, "Red" represents the value of 250, while "Magenta" represents the 0 value.In the second and fourth columns, as it is shown in the colorbar, each pixel indicates one individual cell, one individual cell, S in cyan, P in green, A in yellow, and B in red, an empty position (EP) in blue, or an out-of-dish space (OD) in gray.

Figure 4 .
Figure 4.The system behavior in the face of dynamic "leading" signals.The first row represents the Gaussian, spot, stripe, and reversed-spot patterns as the "leading" signals, S l .The values of the signaling molecules are scaled in the range of [0, 250] to be in agreement with the values of static leading signals.The second row demonstrates the system state after 50 divisions (of progenitor cells).Each pixel in the second row can potentially demonstrates S, P, A, B cell types, an empty position or an out-of-dish space as shown by the colorbar in the second row.The third row shows the abundance of four cell types in the population through the simulation.
(1 − p) * .0.85 * (N) ) number of A (B) cells in homeostasis state.In other words, n * A = p * 0.85 * (N) ( n * B = (1 − p) * .0.85 * (N) ), where N is the number of all pixels in the dish.With this approximated value of n * A , and n * B we set the γ A/B value as follows:

Figure 5 .
Figure 5.The system behaviour in the face of dynamic "leading" signals, S l with spot pattern.A. The first and second rows represent the system state after the 50 and 100 divisions (of progenitor cells), respectively.The first and second columns represent two different formulations of death rate: γ i = n i n * i * v g , and γ i = 3 n i n * i * v g , where i ∈ A, B .As it is shown in the colorbar, each pixel can indicate one individual cell, S in cyan, P in green, A in yellow, and B in red, an empty position (EP) in blue, or an out-of-dish space (OD) in gray.The third row shows the abundance of four cell types through the simulation.B. The template and penalty filters of spot signal pattern.C. The scoring traces corresponding to the first and second columns in panel A. It indicates how the pattern in population in different states of the system through time is aligned with the desired initial signalling pattern.

Figure 6 .
Figure 6.The system behaviour in the face of dynamic "leading" signals, S l with Gaussian, reversed-spot and stripe patterns.The first and second rows represent the template matrices and the system state after 100 divisions (of progenitor cells), respectively.As it is shown in the colorbar, each pixel can indicate one individual cell, S in cyan, P in green, A in yellow, and B in red, an empty position (EP) in blue, or an out-of-dish space (OD) in gray.Here, γ i = 3 n i n * i * v g , where i ∈ A, B .The third row shows the abundance of four cell types through the simulation.The fourth row demonstrates the scoring traces to evaluate the pattern maintenance in the population.It indicates how the pattern in population in different states of the system through time is aligned with the desired initial signalling pattern.

Figure 7 .
Figure 7.The system behabior in the face of Gaussian "leading" signal pattern, and in the absence of physical dish borders.As it is shown in the colorbar, each pixel can indicate one individual cell, S in cyan, P in green, A in yellow, and B in red, an empty position (EP) in blue, or an out-of-dish space (OD) in gray.

Table 1 .
The potential reactions of the system together with their corresponding propensity functions.