A stepwise model of reaction-diffusion and positional information governs self-organized human peri-gastrulation-like patterning

How position-dependent cell fate acquisition occurs during embryogenesis is a central question in developmental biology. To study this process, we developed a defined, high-throughput assay to induce peri-gastrulation-associated patterning in geometrically confined human pluripotent stem cell (hPSC) colonies. We observed that, upon BMP4 treatment, phosphorylated SMAD1 (pSMAD1) activity in the colonies organized into a radial gradient. We developed a reaction-diffusion (RD)-based computational model and observed that the self-organization of pSMAD1 signaling was consistent with the RD principle. Consequent fate acquisition occurred as a function of both pSMAD1 signaling strength and duration of induction, consistent with the positional-information (PI) paradigm. We propose that the self-organized peri-gastrulation-like fate patterning in BMP4-treated geometrically confined hPSC colonies arises via a stepwise model of RD followed by PI. This two-step model predicted experimental responses to perturbations of key parameters such as colony size and BMP4 dose. Furthermore, it also predicted experimental conditions that resulted in RD-like periodic patterning in large hPSC colonies, and rescued peri-gastrulation-like patterning in colony sizes previously thought to be reticent to this behavior.


Suppleme ntary Tables
. B) Simplified mathematical representation of a generic dual inhibition model. C) Gradient formation of free BMP4 ligands as predicted by the model. D) Varying doses and colony sizes demonstrates the inability of the dual-inhibition model to generate a periodic Turing-like response in the free BMP4 distribution.

Background of the model
To test the assertion that a BMP4-NOGGIN RD system was self-organizing the pSMAD1 gradient, we set out to develop an RD model specific to our tissue geometry, initial, and boundary conditions to query testable predictions of signaling gradient formation that may arise in our geometrically confined hPSC colonies. Notably, the purpose of this model was not to attain detailed information of the response of the signal transduction within the differentiating cells, or identify the gene regulatory network that would allow differential induction of fates in response to different signaling levels. Instead, we set out to identify the simplest possible RD model that could produce a BMP signaling gradient for a 1000µm diameter colony when differentiated with 50ng/ml of BMP4 in the induction mediareminiscent of the pSMAD1 signaling patterns observed in our experimental data under those specific conditions. We then validated this model based on the predictions to perturbations of two key experimental parameters ( Fig. 4 main text); and then employed this model to made testable predictions of the differentiation behavior of the hPSC colonies beyond the conditions under which the model was build (Fig. 6, 7 -main text). The RD model described below is a simplified, and an idealized model. Nevertheless, this model generates predictions of morphogen (BMP4) distribution in response to perturbations of experimental conditions which is constantly shown to be consistent with both pSMAD1 gradient formation and associated hPSC differentiation.

Two-component Reaction-Diffusion system
We set out to develop a mathematical model of the concentration profiles of BMP4 and NOGGIN molecules as a function of space and time in our micro-patterned colonies, based on a reactiondiffusion (RD) system (Turing 1952;Murray 2008;Gierer & Meinhardt 1972). The partial differentiation equation (PDE) set can be described as follows:

Supplementary Materials and Methods: Model Description
Here, bmp, and nog are functions of both space, and time that represent the local concentrations of BMP4 and NOGGIN molecules at a particular point in our micro-patterned colonies. F (bmp, nog), and G (bmp, nog) represent the non-linear functions which describe the production rates of BMP4, and NOGGIN. The degradation rates of the molecules are given by dBMP, and dNOG; and DBMP, and DNOG represent the diffusivities of the molecules.
We assumed that the production terms of BMP4, and NOGGIN can be approximated by linear functions (close to the steady state). The nature of the production terms is such that as the values of F(bmp,nog), and G(bmp,nog) increase, the system transitions away from steady state, preventing convergence of the solutions (Murray 2008;Kondo & Miura 2010). Attempts to restrict the values for the morphogen near steady state, to enable convergence, have either used nonlinear functions that saturate at increasing values (e.g. the Hill function (Sick et al. 2006)) or have enforced a range in which the linear approximation of the reaction function is confined (Kondo & Miura 2010). Since we used linear production functions (2) in our model (1), we chose the latter strategy and restricted the reaction functions to a defined range (Kondo & Miura 2010). The production terms are represented by: Changing Variables: To circumvent the issue of intractability of the number of BMP4 and NOGGIN molecules in the circular region of interest modelled by our PDE solutions, we chose to change the variables bmp, and nog into normalized, dimensionless variables which we represent as bmp*, and nog*. The linearized PDE set from (1), and (2) together are of the following form: Taken together, the updated set of partial differential equations with the changed variables (bmp*, nog*, c*BMP, and c*NOG) is given below: * = * + * + * − * + ∇ 2 * (4) * = * + * + * − * + ∇ 2

Dose dependence in BMP4 and NOGGIN production
We observed a BMPi dose dependent production of both BMP4 and NOGGIN (Fig S5). To incorporate this response into our model we chose the following expressions for aBMP, and aNOG.

Initial conditions of BMP4 and NOGGIN distributions in micro-patterned colonies: BMP4
BMP4 is added in the differentiation medium, and presented to the colonies at a uniform dose. Therefore, we considered the initial concentration of BMP4 to be a constant value given by BMPi.

NOGGIN
The initial conditions for NOGGIN are more nuanced. Although NOGGIN is produced in response to BMP signaling in the cardinal 'Activator-Inhibitor' paradigm, BMP4 inhibition in the initial condition (at time t=0), can be achieved by a variety of different molecules (e.g. FOLLISTATIN (FST), CHORDIN, GDF3, and CERBERUS-Like (CERL) among others) in addition to NOGGIN (Wu & Hill 2009). Since we observed elevated basal expression of BMP signaling inhibitors like FST, GDF3, and CERL relative to NOGGIN in hPSCs during basal culture conditions ( Fig. S4B-G), we opted to consider the spatial profile of a 'generic BMP inhibitor' as the initial condition for the RD paradigm.
To identify the specific spatial profile of a generic BMP inhibitor, we developed a simplified model of a passive diffusion-driven profile that would arise in a circular hPSC colony where each cell is a source of the secreted molecule. Over the course of the formation of a confluent hPSC colony, we assumed that the expression profile of a 'generic BMP inhibitor' would reach a steady state.
To approximate this steady state spatial profile, we considered each cell (a point source of the inhibitor), evenly distributed within the colony (Fig. S14A), and assumed an infinite sink at a large distance from the colony (Fig. S14B). Simulation of a steady state response revealed a spatial profile that could broadly be approximated as an elliptical paraboloid (Fig. S14C-D). Accordingly, we considered the initial condition of NOGGIN, which at t=0 can be replaced by the effective contribution of all BMP inhibitors being expressed by the hPSCs in the micro-patterned colony, to be an elliptical paraboloid function (Fig. S14E). Notably the peak concentration chosen for the NOGGIN initial condition is arbitraryand the patterning of free BMP4 distribution was found to be robust to the choice of the peak value.

Boundary conditions for the BMP4-NOGGIN reaction-diffusion system in micro-patterned colonies
We assumed that the cells at the radial edge of the micro-patterned colonies were always subjected to the same concentration of BMP4 that is in the bulk medium, which parallels the 'edge-sensing' model that has been proposed by previous studies (Etoc et al. 2016). For the case of NOGGIN, we assumed a no flux boundary condition.

Initial Conditions Boundary Conditions
The diffusivity values for both NOGGIN, and BMP4 are in realistic ranges (Raspopovic et al. 2014;Sick et al. 2006;Inomata et al. 2013). For instance, Inomata et al. calculated the diffusivity of NOGGIN in Xenopus embryos to be 37±6.6 μm 2 /s, which is close to the predicted diffusivity of Noggin in our system. However, the exact values of the diffusivities in our system have not been measured.
The value for bNOG was chosen to be zero since NOGGIN does not have any receptors and is therefore, unable to repress its own expression.  Fig. S15, we show the response of the predicted spatial profile of free BMP4 molecules within the geometrically-confined hPSC colony to varying the above parameters to provide a sense of the sensitivity of the model output to the model parameters. The sensitivity of the model to perturbing the mesh definition is shown in Fig. S16.  Figure S15: Response of predicted gradient formation to perturbation of model parameters.