{'fraction_non_alphanumeric': 0.04517583656410679, 'fraction_numerical': 0.01132832907146196, 'mean_word_length': 4.271778518625887, 'pattern_counts': {'":': 1, '<': 5, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Reaction-diffusion models for morphological patterning of hESCs

In this paper we consider mathematical modeling of the dynamics of self-organized patterning of spatially confined human embryonic stem cells (hESCs) treated with BMP4 (gastruloids) described in recent experimental works. In the first part of the paper we use the activator-inhibitor equations of Gierer and Meinhardt to identify 3 reaction-diffusion regimes for each of the three morphogenic proteins, BMP4, Wnt and Nodal, based on the characteristic features of the dynamic patterning. We identify appropriate boundary conditions which correspond to the experimental setup and perform numerical simulations of the reaction-diffusion (RD) systems, using the finite element approximation, to confirm that the RD systems in these regimes produce realistic dynamics of the protein concentrations. In the second part of the paper we use analytic tools to address the questions of the existence and stability of non-homogeneous steady states for the reaction-diffusion systems of the type considered in the first part of the paper. We find sufficient conditions on the data of the problem under which the system has an universal attractor.


INTRODUCTION
The use of reaction-diffusion (RD) system to model morphological patterning of developing organism started with a seminal work of Turing [19], who observed that in some cases the homogeneous distributions of chemical with coupled reaction rates are not stable, and upon small perturbations give rise to sizable, non-homogeneous patterns of chemicals. Cells at the peaks of Other mechanisms of morphogenesis such as of mechanical and mechanochemical types had been developed, see for example Murray at al. [13], Murray [12]. These models take into account not only chemical but also mechanical properties of treated cell samples, such as, for example, the motion of chemical in the sample induced by the growth of cells.
In this paper we will discuss recent experimental findings by Warmflash et al. [20], Chhabra et al. [4] and Heemskerk et al. [7] on cell fate differentiation of human embryonic stem cells (hESCs) that do not fit into the Turing paradigm of morphogenesis, or the other mechanisms mentioned above. In a typical experiment, a spacial confined colony of cells is treated with BMP4 (bone morphogenetic proteins) which leads to differentiation of cells into the three embryonic germ layers: endoderm, mesoderm, and ectoderm, surrounded by an outer ring of extraembryonic cells.
It has been established in earlier works, see for example Arnold et al. [1], that BMP4 results in expression of Wnt and Nodal proteins that are essential for formations of germ layers.
The evidence accumulated in [4] indicates that fate differentiation occurs not after the formation of stable patterns of BMP4/Wnt/Nodal, but concurrently, during the propagation of signaling waves of Wnt and Nodal, unlike the pre-patterning in Turing-type process. Additionally, the terminal, stable distributions of these proteins do not correspond (decisively) to location of the germ layers. For example, Wnt and Nodal signaling, which synergize to generate mesendodderm, both spread into the middle region where ectoderm forms. We mention in passing, that the distributions of proteins do tend to non-homogeneous, almost radial, steady states. Moreover, Chhabra et al. [4] performed a series of experiments ruling out the cell motion and cell growth as effective mechanisms of morphogenesis. The later facts points strongly in favor of a RD system as correct mathematical model.
We summarize below some of the characteristic features for the dynamics of distribution of proteins during the patterning, obtained in Chhabra et al. [4] and Heemskerk et al. [7], that we address in this paper.
(1) Initially high and uniform over the entire domain, BMP4 signaling activity evolves into a region of high activity near the boundary and low activity in the middle and central parts of the domain; (2) boundary of BMP4 initiates waves of Nodal and Wnt that move into the interior towards the center of the domain; (3) propagation of the Nodal wave proceeds independently of BMP4 (and Wnt), after a certain activation period; (4) distributions of BMP4, Wnt and Nodal activities tend to steady states, by the end of the experiment, with the peaks at the boundary, middle of the domain and the center of the domain, respectively. The difference between the peak value of Wnt and its value at the center smaller than the difference between the peak value and the values near the boundary.
The final distributions of all proteins appear to be radial, non-homogeneous steady states.
While the understanding of all mechanisms is far from being complete, an effort was made to generate dynamics consistent with (1)- (4) in the framework of activator-inhibitor RD systems. Tewary et al. [18] developed a RD model model that produces realistic patterns of BMP4, with realistic dependence on parameters, such as the size of the cell colony. However, due the fact that the reaction part of the model is linear, the model requires selection of matching boundary and initial conditions, that do not reflect the state of the problem at the beginning of the experiment. For example, it assumed that BMP4 inhibitor does not diffuse from the sample boundary, and, initially, peaks at the center of the sample. Thus, the model is partially "exogenous".
Chhabra et al. [4] addressed the mathematical modeling of dynamics described in part (3) of the above list. In this model, BMP4 acts as an activator for Wnt, which in its turn activates Nodal. The key assumption is a reaction term the equation for Nodal, that incorporates a threshold parameter, depending on the concentration of Nodal, that switches production from being Wnt dependent to auto-catalytic. Other assumptions include structurally different mechanisms of activation/inhibition for Wnt and Nodal.
The analysis of RD systems in the above mentioned papers relies on the numerical solutions of the corresponding systems of PDEs.
The purpose of the present paper is to develop closed form ("endogenous") models based on the classical Gierer-Meinhardt activator-inhibitor system, see equations (4), (5) below, that reproduce behaviors (1)-(4). By the closed form we mean, a solution to an initial-boundary value problem where the initial and boundary conditions reflect the actual experimental setup.
In particular, we assume that all substances can diffuse off the colony edge, and use "the Newton's law of cooling" with appropriate background values of the substances.
In the first part of the paper, using numerical simulations, we show that realistic dynamics of the proteins can be obtained by choosing suitable reaction coefficients in the activator-inhibitor system. The dynamics of BMP4 is best described by a system with a single, stable node, to which we will refer as type 1 system, see section 2. The dynamic of Wnt and Nodal systems fit to the patterns produced by type 2 systems, which have a stable node, a stable focus and a saddle point in the phase plane. Moreover, we show that different behaviors of Wnt and Nodal, as described in (5), can be attributed to the size of the reaction coefficients alone. This is due to a general fact that scaling reaction coefficients in a RD systems, which leaves the phase portrait unchanged, while retaining the same diffusion coefficient, results in a different dynamics and, in particular, in different steady states. Thus, we provide another explanation of phenomenological properties (1)-(4) based on the dynamical differences structurally similar system of PDEs. The analysis can be useful in providing an estimates on the ranges of the reaction coefficients that distinguish BMP4, Wnt and Nodal at the level of activator-inhibitor systems.
In the second part of the paper we address the stability property in part (4). The main mathematical difficulty comes from the fact that that the RD models (equations + boundary conditions) in question do not have, in general, homogeneous steady states. Thus, linearizing equations on a constant state, and solving for eigenvalues and eigenfunctions does not provide meaningful information, because the constant states are not solutions, the fact that sometimes is overlooked in biological literature. The growth of oscillations should be measured with respect to a proper steady state, which, in this case is non-homogeneous.
This leads us to the following problems that we address in this paper: existence of steady states and do exist, and sufficient conditions for stability.
There are local in time, unique, classical solutions to RD systems in the H'older space C 2+α,1+α/2 , as was proved by Ladyzhenskaja et al. [11] for a more general system of parabolic equations. Estimates on the max norm of the solution is needed to extend solutions for all times t > 0. This was done by Rothe [16] for activator-inhibitor systems with no-slip boundary conditions, which does not apply in our case. Moreover, the estimates obtained in [16] depend on the diffusion coefficients, which greatly complicates the asymptotic analysis. We note here, that the special structure of the activator-inhibitor equations does not allow the application of the theory of invariant regions of Chueh et al. [5], another well-known technique for the asymptotic analysis, see for example a book of Smoller [17].
Our approach is to rely on the maximum principle for the parabolic equations to obtain the uniform in (x,t) estimates of the solution. The application of this method places some restrictions on the coefficients of the RD system. The key point here is that the bounds are independent of the diffusion coefficient.
Then, we use the energy-types estimates for L 2 norms of the solution and it time derivative to identify the suitable stability condition that implies the exponential decay of the norm of the time derivative. Further analysis required to bound the spatial gradient of the solution in L 2 norm, with an upper bound, independent of time. The later fact allows us to extract a strongly convergent sequence u(x,t n ) with t n → ∞, whose limiting point is a steady state solution of the activatorinhibitor system. Here, u is a solution vector u = (u, v). Finally we show that such steady state, u s , is the limit of the u(x,t), at t → ∞. The stability condition mentioned above, for given source terms f (x), g(x), size of the domain Ω and the diffusion coefficient µ restricts the size of the initial data max Ω |u 0 |. This condition defines the basin of attraction for the steady state, i.e., there is a ball B in L ∞ , centered at zero, such the steady state u s belongs to B and for any initial data u 0 ∈ L ∞ ∩C 2+α , solution u(x,t) converges to that steady state.

RD MODELS OF ACTIVATOR-INHIBITOR TYPE
We consider an activator-inhibitor RD system where µ u , µ v are positive diffusion coefficients, a, c > 0 are rates of decay and b, d are positive reaction coefficients. This system was introduced by Gierer and Meinhardt [8]. The inhibitor enters the first equation through the factor (1 + v) −1 , where we added 1 avoid infinitely high rates when the level of the inhibitor v is small. This is a natural assumption in the experiments described in [20,4].
Here we denote by f (u, v) and g(u, v) the right-hand side of equations (1), (2). There is a single fixed point (stable node) In this case (0, 0) is a stable node, as the gradient matrix At other two points, the gradient matrix can be computed to be equal to where u is evaluated at a fixed point. Thus, we obtain the values of the trace and the determinant of the gradient matrix as These values determine the type of each fixed point. Substituting the values for u 1,2 from (3) we find that Det = ± (ac) 2 − 4dca 4 b −2 . This gives a saddle point, when the value of the determinant is negative. The remaining point is either a stable focus (a < c), a center (a = c), or an unstable focus (a > c). For the non-linear system (1), (2), the last two possibilities result in periodic motion. This type of motion is not observed in the experiment discussed in this paper, so we assume that a < c. The null clines for different values of the decay/reaction coefficients are sketched in Figure  1, which we use to distinguish the corresponding RD systems as type 1 and type 2. We will use system (1), (2) to model the dynamics of BMP4 and its inhibitor Noggin. To study the signaling waves of Nodal and Wnt we will use a source term f (x) in the activator equation that models the influence of BMP4 on production of Wnt: Here u is the concentration of Wnt and v the its inhibitor DKK. In the numerical simulations we will assume that f (x) is concentrated near the boundary of the domain which reflects the experimentally observed distribution of BMP. The activation of Nodal occurs through BMP → Wnt → Nodal signaling pathway. We will model this by a simplifying to BMP → Nodal signaling, and using equations (4), (5), with a different set of reaction coefficients, to model the dynamics of Nodal ant its inhibitor Lefty.
As we will show in by numerical simulation, the system of BMP4 and its inhibitor has a better fit into the reaction system of type 1, see Figure 1, while the systems for Wnt and Nodal are better described by type 2 dynamics.
The initial data correspond to the high initial concentration of BMP4 and low (zero) concentrations of other chemicals, in accordance with experiments described in [20].
To complete the model, we need to postulate boundary conditions for chemical concentrations. It should be emphasized that the boundary conditions are the integral part of the solution, that plays an important part in the way the dynamics proceeds. According the experimental set up the chemicals can diffuse from the domain of a cell sample, the intensity of this "leaking" being proportional to the difference between the boundary concentration of the chemical and the "background" concentration. For BMP4 it is reasonable to take the background concentration at the fixed levelū that equals to the initial concentration of BMP4. For other chemicals in question, the background concentration is zero. This is so-called "Newton's law of cooling." It is expressed as where n is the external, unit normal vector at ∂ Ω . Similarly, for the inhibitor, We note that an earlier works on mathematical modeling of hESC development such as [18], used ad hoc boundary conditions, not consistent with the experimental setting. An alternative way to deal with the boundary, is to embed the reaction domain into a larger domain where chemicals are only diffused, and postulate, for example, no-flux boundary conditions on the larger domain, as was done by in [4]. The reason being that, the precise form of the boundary conditions on the large domain should have minimal effect on the domain where the reaction takes place.

NUMERICAL SIMULATIONS
Simulations are performed for RD systems written in scaled variables, using the values of coefficients of the magnitude typically occurred in experimental studies, see Appendix for details. In particular, the computational domain is a disk of radius 1, and time t is measured in days. The values of parameters for system (1), (2) for BMP4/Noggin, and system (4), (5) for Wnt/DKK and Nodal/Lefy, used in the simulations, are listed in table 1. The parameters are selected is such a way that RD system for BMP4 is of type 1, and the systems for Wnt and Nodal are of type 2. Reaction coefficients for Nodal system differ by a factor of 0.4 from the corresponding coefficients for Wnt system, which means that the phase portraits for the reaction dynamics are identical in both cases.
The source term f (x) in (4) is set to be supported near the boundary of the disk: Finally, the background stateū is set to 3 for BMP4, as well as the initial values for BMP4. All other chemicals have zero initial values, and zero background states. We will show below that behaviors 1.-4., listed in the Introduction, are captured by the RD models described here. The exact timing of different phenomena described below does not necessarily correspond to experimentally observed values. That would require more precise estimation  (1), (2) and (4), (5), used in numerical simulations.
of the parameters of the model. Our main goal is to establish that qualitatively correct behavior is produced by the model.

3.1.
Terminal concentrations of proteins. The numerical simulations show that concentrations of all three proteins approach a radial, steady state profiles by t = 1 day. Figure 2 shows 2d and 3d plots of concentrations of the proteins at t = 3 days. BMP4 is concentrated at the boundary of the domain. Wnt peaks in the middle section, but takes comparable values at the center, and Nodal peaks at the center but somewhat extends to the middle section of the domain. The appearance of steady states was identified when the change between the successive iterations of the numerical solution became less than 10 −6 units, over a period of time of 1 day.
To illustrate the difference between three different diffusion-reaction regimes we map the radial cut of each protein and its inhibitor in the phase plane on Figure 3. Point B indicates the values at the boundary of the disk and point C represents the values at the center. The plots also show the stable fixed point, S, for the corresponding reaction dynamics. BMP4 starts at non-zero value at the boundary due to the influence from non-zero background stateb and then moves towards zero, according to the reaction, as one moves to the center of the disk. Point S is a stable focus for Wnt/DKK and Nodal/Lefty pairs, however the reaction coefficients are stronger in Wnt system than in Nodal system. This results in the radial profile of Wnt/DKK being "bent" in the direction of the reaction. This property results in maximum of Wnt to be located in the middle of the disk, whereas Nodal has maximum at the center.

3.2.
Signaling waves of Wnt and Nodal. Figure 4 shows the evolution of concentrations of Wnt and Nodal as functions of time. In both cases, the concentration first increase near the boundary, where f (x) is supported and then they move toward the center. This behavior corresponds to the "signaling waves" of the proteins described in Introduction.
3.3. Effect of inhibition of BMP4. Figure 5 shows the effect of inhibition of BMP4 ( f (x) = 0) at time indicated by the variable tcut on the shape of the final concentration at time t = 3 days.
The simulations show that there is a critical time t 0 with the property that if BMP4 is inhibited prior to t 0 the system converges to zero steady-state, but when BMP4 is inhibited after t 0 the system proceeds autonomously to a non-homogeneous, non-zero steady state. For Wnt dynamics t 0 is estimated to be between 0.001 and 0.005 days, and for Nodal, between 0.005 and 0.01 days, the difference is due to the difference in the magnitude of the reaction coefficients.
3.4. Dependence on parameters. For small variations of parameters given in Table 1 the numerical simulation produce qualitatively similar results, indicating that the problems are stable. That is, the terminal steady state concentrations are stable. This property is lost when the larger variations. We performed the numerical simulation of BMP/Noggin dynamics with large gap in diffusion coefficients, by selecting µ u = 1 µm 2 /sec and µ v = 55 µm 2 /sec, (instead of µ u = 11 µm 2 /sec, and µ v = 55 µm 2 /sec used previously), while keeping all other parameters. Figure 6 shows 2d and 3d plots of non-radial profile of BMP4 at time t = 3 days. The non-radial perturbations start to develop from a radially symmetric profile at the time about t = 0.2 days. Note that, due to the radial symmetry of equations, the problems has a unique, classical, radially symmetric solution if the initial data have this property, but the numerical solution deviate from it significantly.
A possible explanation of this phenomenon is that the problem has an unstable radially symmetric steady state to which the system moves from its initial values. Small deviations from radially symmetry due the numerical approximation lead to the growth of perturbations shown on the figure. That is, this is the case of Turing instability.
Interestingly, the instabilities appear to be restricted to the boundary, and further simulations (not shown here) produce a different number of peaks, with further variations in the diffusion coefficients. This non-homogeneous profile might, in principle, be associated with the formation of the outer ring of germs in a cell colony. Further investigation of a coupled BMP4-Wnt-Nodal system is needed to clarify if this behavior bears some significance in actual biological processes.

EXISTENCE AND STABILITY OF STEADY-STATE SOLUTIONS
One of the main features of the models considered in this paper and experimental works cited in the introduction is the formation non-homogeneous steady-states of the chemical concentrations. On the middle and bottom plots, S is a stable focus with counterclockwise rotation. Plots also show velocity fields of each reaction system.
Moreover, this steady-states appear to be stable, as none of the Turing-type instabilities is observed in the experimental setting, see Chhabra et al [4]. In this section we address the question of existence of stable steady-state solutions. The results that we prove below, apply to RD systems in the form: with the Direchlet boundary conditions This can be seen as a limiting case of the boundary conditions from earlier sections, when the rates of cooling h u , h v → +∞, i.e., there is high rate of transfer of chemicals to or from the background state. We will assume that decay and reaction coefficients are positive. We would like to compare the method we use with the well-known method of invariant regions by Chueh et al. [5], for establishing time asymptotic behavior of solutions of reaction-diffusion systems. The latter method applies to reaction-diffusion systems in the form where u ∈ R n , (x,t) ∈ Ω × [0, +∞), D is n × n, diagonal matrix with non-negative entries, and the vector source term F ∈ R n . The system is supplied with the zero-flux boundary conditions For the method to work the must have bounded invariant regions (see below) in order to establish bounds on supremum norm of |u(x,t)|. Next, the method relies on the fact that time asymptotic behavior of solutions of (9) can be compared with the solution of the system of ODEs: meaning that the limiting behavior of (9) is a homogeneous (constant in x) state.
If we look at the system (6), (7), we see that the right-hand side explicitly depends on x through functions f (x) and g(x). The boundary conditions (8) differ as well. Thus, in general, the steadystates of the problem (6), (7) and (8) are non homogeneous (non-constant).
We show now that system (6), (7) does not have invariant regions either, so that the uniform estimates must be obtained by other means.
A closed set S ⊂ R n is called an invariant set for (9) if for any t > 0, x ∈ Ω the solution u(x,t) ∈ S whenever the initial data u 0 (x) ∈ S, for every x ∈ Ω .
Suppose that S can be written as an intersection of "half-spaces": where G i are smooth functions. Theorem 14.14 from [17] gives a sufficient and necessary conditions for S be an invariant regions. (9), if and only if for every u ∈ ∂ S, (so that G i (u) = 0, for some i,):

Theorem. S is an invariant region for
(1) ∇G i (u) is the left eigenvector of D; (2) G i is quasi-convex at u; Recall that function G i (u) is called a quasi-convex function at u, if for any vector v such that Consider now system (6), (7). In this case D is a diagonal matrix µI, where I is 2 × 2 identity matrix. The first condition of the theorem implies that ∇G i (u) is proportional to vector (1, 0) or (0, 1), that is, the level sets of G i are either horizontal or vertical lines in (u, v) plane. Thus, an invariant region, if it exists, is a rectangle. The third condition of the theorem then implies that on a line {u = const.}, we must have and on a line {v = const.}, It can be seen from Figure 1 that it is impossible as non of these functions changes its sign.

Main theorem.
We will use the standard notation for the spaces of continuous, H'older continuous functions, as well as L p spaces. For definitions, we refer readers to [11]. Norms in L p (Ω ) space will denoted by u p , 1 ≤ p ≤ ∞. We let Q T = Ω × (0, T ) and Γ T = (Ω × {t = 0}) ∪ (∂ Ω × [0, T ]), for T > 0. Now we state our main result, that we will be proved below.
, and necessary compatibility conditions between the initial and boundary values hold. Then, there exists a unique classical solution of the system on Ω × [0, ∞). For any T > 0, u, v ∈ C 2+α,1+α/2 (Q T ), the following properties hold. (1) u(x,t), v(x,t) are non-negative and bounded above with a constant independent of x,t, and µ.
such that if an inequality

holds, where C(Ω ) is a constant from the Poincare's inequality, then, there exists a steadysolution (u s (x)
, v s (x)) of (6), (7), and (3) Under condition (10), there is a ball B = B( f , g, µ) ⊂ L ∞ (Ω ) such that if the initial datã u 0 ,ṽ 0 ∈ B and verify all other properties of the initial data stated above, then for the corresponding classical solutionũ(x,t),ṽ(x,t),

4.2.
Proof of main theorem. The proof is given below in a series of lemmas, where (u, v) is a local, classical solution of the problem. We will make a repeated use the of the maximum principle, that can be found, for example, the book by Evans [6].
Since f (x), g(x) ≥ 0, it follow from this lemma that u(x,t) and v(x,t) are non-negative functions. We proceed with uniform upper bounds. C > 0, depending on a, b, c, d, max u 0 , max v 0 , max f (x) and max g(x), but not µ, such that (11) max

Lemma 2. There is
where C is polynomial function of its arguments, with positive coefficients, and is independent of µ.
Proof. Let Z be a solution of the ordinary differential equation, As for v, we have, Let φ be a smooth, non-increasing function that will be chosen later. For function φ (v(x,t)) we Adding the last equation to (12) we get, for some α depending on a, b, c, d. Let W be a solution of with initial condition W (0) = 0, i.e., W (t) = α 1 max Q T ln(1 + v) (1 − e −at ) , with α 1 = α/a. Subtracting equation for W from (14), we obtain: Now, using the maximum principle (1) we obtain Therefore, for any (x,t) in the domain Q T ,û + φ −W ≤ max Ω u 0 , or, for some α 2 depending on a, b, c, d.
Using a maximum principle again we get By the elementary properties of function ln(1 + v), we find that max Q T v is bounded by a polynomial with positive coefficients in variables of max u 0 , max v 0 , f ∞ , g ∞ . The corresponding estimate for max Q T u follows from this and (15). Now, the global existence follows. Proof. We will use the following characterization of time maximal time of existence T of a local solution, from Rothe [16], theorem 1, page 111. It is proved there that if T < +∞ then the maxnorm over x of (u(x,t), v(x,t)) grows without bound as t approaches T . But this can not happen due to the estimates derived above in (11). Therefore, the contradiction leads us to conclude that the classical solution in fact exists for all times t > 0.
To show that the classical solution (u, v) of the reaction-diffusion system settles on a steady state it sufficient to show that the time derivative of the solution converges to zero. We will use an energy-type estimate to establish this fact. The proof makes use of the Poincare's inequality that we state for a reference below.

Lemma 4.
Let Ω be a bounded, connected, open subset of R n with a C 1 boundary ∂ Ω . Let 1 ≤ p < ∞. Then there exists a constant C, depending only on n, p and U , such that for any integrable function u with ∇u ∈ L p (Ω ) and zero trace on the boundary ∂ Ω , Proof can be found in Brezis book [2].

) is a polynomial function of its arguments, and is independent of µ, and T, and C(Ω ) is a constant from the Poincare's inequality. If the stability condition
holds, then Proof. Taking time derivative of both sides of the equation for u (1), we get, Multiply with u t and integrate over the domain to get, Using integration by parts, We will use Young's inequatlity, uniform bounds on u and v, and Poincaré's inequality (16) with p = 2 applied to the function ∂ t u (notice, that ∂ t u equals to zero on the boundary of the domain) to get the next result. The constant C(Ω , 2) from that inequality will be abbreviated to C(ω).
where c 0 has polynomial dependence on max u, max v. Similarly, for v, Now the statement of the lemma follows from the last inequality and uniform bounds from lemma 2 When uniform estimates from lemma 2 are substituted into function C s in lemma 5, condition (17) defines the range of L ∞ norms of admissible data f ∞ , g ∞ , u 0 ∞ , v 0 ∞ . We can state this in another way, by saying that if f and g are such that stability condition (10) holds with u 0 = 0 and v 0 = 0, then there is a ball B = B( f , g, µ) ⊂ L ∞ (Ω ), centered at zero such that the same condition (10) holds for any u 0 , v 0 ∈ B. In the remaining part of the proof we will assume this condition.
Next we will obtain bound on the gradients of ∇u, ∇v.
Lemma 6. For any t > 0, with a positive C, independent of time.
Proof. We multiply for u by u t and integrate by parts to get, Using boundary condition on u t and uniform bounds on u and v from lemma 2, with an appropriate c 0 . Similarly, for v, we get, Adding the last two inequalities we get, We integrate this inequality in time from 0 to t, use exponential decay estimate on u t 2 , v t 2 from previous lemma, together with uniform bounds on u and v to get: with some C > 0, as stated in the lemma.
Let u(x,t), v(x,t) be the solution of the section from previous section. Let t n be a non-decreasing sequence of time converging to +∞. Consider sequences of functions {u n (x) = u(x,t n )} and {v n (x) = v(x,t n )}. From the estimates of u, v and their gradients, it holds that there is C independent of n such that u n 2 ≤ C, v n 2 ≤ C, ∇u n 2 ≤ C, ∇v n 2 ≤ C.
We will need the following compactness result, the proof of which can be found in chapter 5 of Evans book [6].

Lemma 7.
Assume Ω is a bounded open subset of R n , and ∂ Ω is C 1 . Suppose 1 ≤ p < n. Then, there is a compact embedding Since Ω is a bounded set it follows this theorem that both sequences are pre-compact in L 2 (Ω ). This means that there is a subsequence of {t n }, that we still label by n, and two functions u, v ∈ W 1,2 (Ω ) such that lim n→∞ u n = u s , lim n→∞ v n = v s in L 2 norm. From the convergence in norm, it follows that a further subsequence can be extracted such that u n and v n converge to u s and v s almost everywhere in Ω . Moreover, since L 2 (Ω ) is a reflexive space and ∇u n , ∇v n are from bounded sets, there is still further subsequence such that ∇u n → ∇u s , ∇v n → ∇v s , weakly in L 2 (Ω ).
Notice, also, that as time derivatives of u, v are bounded, in L 2 norm. Now we pass to the limit in the equations.
Let ω(x) be a smooth test function, equal to zero on the boundary ∂ Ω . From the original reaction-diffusion system, considered at times t = t n we obtain the following integral relations: Passing to the limit in each term of these equations, using above compactness properties we obtain that µ ∇u s · ∇ω dx = −au s + bu 2 i.e, (u s , v s ) is a weak solution. As a pointwise limit of u(x,t n ), v(x,t n ), (u s , v s ) take boundary values u b and v b . By the well know regularity results for elliptic equations, it follows that u s , v s ∈ C 2+α (Ω ), and it is classical solutions of the same system. Now we prove the following.
Lemma 9. As t → ∞, u(x,t), v(x,t) of the system converges to u s (x), v s (x) in L 2 norm: Proof. Suppose that (u(x,t), v(x,t)) does not converge to (u s , v s ) in L 2 norm. Then, there is a sequence of times t n and ε > 0 such that Using the arguments of this section we conclude that there is another steady-state (ũ s ,ṽ s ) and Since (u s , v s ) and (ũ s ,ṽ s ) solve the same system of equations, subtraction corresponding equations, multiplying them byũ s − u s andṽ s − v s , and integrating over Ω , after simple manipulations we get where C s and C(Ω ) as in (17). Since µC(Ω ) −C s is positive, we conclude thatũ s = u s andṽ s = v s . This clearly contradicts statement (21) and the lemma is proved.
In the next theorem we show that the steady state (u s , v s ) is stable and "attracts" solutions of the reaction-diffusion system, with the same source terms f (x), g(x), provided that the solution of the latter satisfy the stability condition (10).

CONCLUSIONS
In this paper we address the mathematical modeling of recent experimental studies on selforganization of human embryonic stem cells during early stages embryo's development. Although several models based on reaction-diffusion equations were proposed in literature, those results are only partially satisfactory as they either contain a number of artificial assumptions on the reaction part of the model or use initial and boundary conditions that do not correspond to the experimental setup.
We showed that an Gierer-Meinhardt system of reaction-diffusion equations, with properly selected reaction coefficients and supplemented with Robin-type boundary conditions, qualitatively reproduces many of the experimental findings, thus identifying a proper mathematical framework. In this paper we only present numerical results for circular domains for brevity of presentation. Additional experiments with domains of irregular shapes, including non-convex domains, were presented by Bedekar [3]. All numerical simulations confirm very good qualitative agreement between our models and in vitro experiments. Moreover, the model produces various new phenomena for the reaction-diffusion system under investigation such as an interesting instability investigated numerically in section 3.4.
The second part of the paper is motivated by the numerical results obtained in the first part, and addresses the existence of of non-homogeneous steady state solutions and the asymptotically Parameters Values SI units µ bmp 11 µm 2 /sec µ nog 55 µm 2 /sec λ bmp 9 × 10 −4 1/sec k bmp 9 × 10 −4 1/sec λ nog 9 × 10 −4 1/sec k nog 9 × 10 −4 1/sec attract solutions of the reactions-diffusion system. In general, this is a hard mathematical problem, which we were able to address under certain simplifying assumptions about the system. With proper ramifications, the model considered in this paper can potentially lead to important scientific insights into the behavior of the biological system. In particular, instabilities outlined in numerical experiments performed here warrant careful analytical investigation. In addition, we can use experimental data to estimate parameters in the PDE model via a Bayesian approach and use the resulting realistic model to predict outcomes of experiments in domains of various sizes. We intend to carry out further investigation of the reaction-diffusion model presented here in subsequent papers.

APPENDIX
We use the following activator-inhibitor system for the dynamics of BMP/Nogin: whereũ,ṽ are some reference values for BMP4 and Noggin. The boundary conditions are whereū is the background value of BMP4, and H bmp , H nog are positive numbers. The initial conditions: u(x, 0) =b, v(x, 0) = 0, which correspond to a cell colony being treated with high concentration of BMP4. The typical magnitudes of the parameters are listed in Table 6. The colony size (radius of the disk) L = 500 µm, and a typical experiment takes up to 3 days (3τ, τ = 86400sec). The experimental data on the values of H bmp and H nog are not available. We set them to 1 (µm) −1 .
For the boundary and initial conditions, ratioū/ũ = 3. The scaling for Wnt/DKK and Nodal/Lefty RD systems are similar.
The numerical simulations are based on the forward Euler approximation of time derivatives with finite element methods, using piece-wise linear functions for the space discretization, and a suitable triangulation of the domain. Space and time partitions steps, (h, δ ) were set to h = 10 −3 and δ = 10 −6 , with δ = h 2 , to prevent numerical instabilities. The method was implemented using FreeFem++, see Heicht [9].