‘Neighbourhood watch’ model: embryonic epiblast cells assess positional information in relation to their neighbours

In the late 1960s, Lewis Wolpert introduced the concept of ‘positional information’, asking the question of how cells within a morphogenetic field could adopt several cell-type identities in response to signalling cues from the embryo. The analogy of a French flag, with three colours: red, white and blue, was used to symbolise the cell types (Wolpert, 1968, 1969). Wolpert proposed that a gradient of a hypothetical ‘morphogen’ diffusing away from a local source and decaying with distance would be read by cells, which respond with discrete thresholds to adopt the various identities. He named this the ‘French Flag problem’.

Since then, several systems have been found in which a morphogen imparts positional information resulting in a defined morphological pattern. These include head and foot formation in Hydra (Schaller, 1973; Bode, 2011), patterning of the wing (Lecuit et al., 1996), leg and antenna imaginal discs of the fly (Postlethwait and Schneiderman, 1971) and limb regeneration in vertebrates (Kumar et al., 2007). Various mechanisms have been studied by which cells might interpret such morphogen gradients so that their positions are defined precisely and robustly. In cultured cells and explant systems (Gurdon et al., 1999, 1995) it has been shown that cells respond directly to morphogen concentration, in a manner most similar to that described by Wolpert (1969). In vertebrate neural tube patterning, the gradient of Shh is transformed into a dynamic profile of Gli (a transcription factor) to generate spatial patterns of downstream gene expression, suggesting that cells interpret positional information using intracellular regulatory networks, in which a temporal element is important (Dessaud et al., 2010; Cohen et al., 2013). The bicoid gradient, which sets up the anterior-posterior axis in fruit fly embryos, has been studied extensively (Driever and Nüsslein-Volhard, 1988a,b; Gregor et al., 2007a,b) and it has been suggested that spatial averaging across nuclei is one mechanism by which noise is reduced in the transduction of the bicoid signal (Gregor et al., 2007a).

All the above systems are relatively small (<100 cell diameters) (Wolpert, 1969) allowing stable gradients to be set up which span the entire field. However, some developing systems are much larger in size, begging the question of what mechanisms cells might use to assess their positions. An example of such a large system is the early chick embryo just before the onset of gastrulation. The embryo contains as many as 20,000-50,000 cells and is ∼3 mm in diameter. Within this large field, the primitive streak (the site of gastrulation) can arise at any point around the circumference. Any isolated fragment of this large embryo can initiate primitive streak formation; however, only one primitive streak forms, suggesting the existence of patterning events that coordinate cell behaviours across the whole field.

In these early embryos, the ‘pattern’ is established in the marginal zone, a ring-like region of extraembryonic tissue, lying just outside of the central disk-like area pellucida, where the embryo will arise. The primitive streak, the first indication of the future midline of the embryo, arises at one edge of the inner area pellucida, adjacent to the posterior part of the marginal zone, where the TGFβ-related signalling molecule cVG1 (also known as GDF3) is expressed. Previous studies have shown that positioning of the primitive streak requires ‘positive’ inducing signals by cVG1/NODAL from the posterior marginal zone near the site of streak formation, and that this is antagonised by BMP signalling, which is highest at the opposite (anterior) end of the blastoderm (Fig. S1A) (Shah et al., 1997; Streit et al., 1998; Bertocchini and Stern, 2012; Streit and Stern, 1999; Bertocchini and Stern, 2002; Skromne and Stern, 2002; Bertocchini et al., 2004). The distance between the two extremes of the marginal zone is ∼300 cell diameters. Previous studies have suggested that these signals are part of a ‘global positioning system’ to establish polarity in the chick embryo (Bertocchini and Stern, 2012; Arias et al., 2017) and therefore that the whole embryo is a coordinated system of positional information.

To find out how cells interpret morphogen concentrations to generate positional information, we designed two computational models to represent either a fixed gradient, read locally by cells, or a system in which cells compare themselves with their neighbours to determine their position in the field. Using a combination of embryological manipulations and computational modelling, we ask which of these two models can best account for the results of various manipulations in the spatial distribution, number and intensity of the inducing (cVG1/NODAL) and inhibitory (BMP) signals. We find that the ‘positional information’ that determines the site of primitive streak formation is explained better by a mechanism by which cells compare themselves to their neighbours rather than by a cell autonomous assessment of gradients. We name this the ‘neighbourhood watch’ model.

When a small pellet of cVG1-expressing cells (HEK293T cells transfected with a cVG1-expression construct) is grafted into the anterior marginal zone (the innermost extraembryonic epiblast, just outside the central embryonic area pellucida), it can initiate formation of an ectopic primitive streak that eventually develops into a full embryonic axis (Shah et al., 1997; Skromne and Stern, 2002). However, endogenous cVG1 mRNA is expressed as a crescent encompassing an arc of ∼60° in the posterior marginal zone (Fig. S1A). To mimic this distribution more closely, as well as to test the effects of a higher concentration of cVG1-inducing signal, we placed two cVG1-expressing cell pellets side-by-side in the anterior marginal zone, and assessed initiation of primitive streak formation by in situ hybridisation for BRACHYURY (cBRA; TBXT) after overnight culture (Fig. 1A-D). Only a single ectopic primitive streak was generated near the middle of the two cVG1 pellets (Fig. 1B); neither double nor thicker ectopic streaks were observed, similar to the effects of implanting a single pellet.

To provide a stronger and wider signal, we tested the effect of implanting four cVG1-expressing cell pellets side-by-side in the anterior marginal zone. Surprisingly, in the majority of cases (11/16 embryos), no ectopic primitive streak formed and no ectopic cBRA expression was seen (Fig. 1E,H,K). As application of the equivalent of a quad-dose of inducer spread over a fourfold wider area does not cause either more efficient or wider induction than a single dose, we speculated that ‘boundaries’ to the signalling domain may be required. To test this, we placed a control cell pellet (HEK293T cells transfected with pCAβ-GFP; see Materials and Methods) to split four cVG1-expressing cell pellets into two groups on either side. The incidence of ectopic streak formation doubled (Fig. 1F,I,K). If a boundary is indeed important, we might expect that, perhaps paradoxically, ectopic streak induction might increase if a pellet expressing the inhibitor BMP4 (rather than a control pellet) is used to interrupt the set of four cVG1-expressing cell pellets. This is indeed the case (Fig. 1G,J,K). Together, these results suggest that cells may sense variations in signal strength in relation to their neighbours, rather than measuring the absolute amount of local signal they receive, to determine the outcome of the inductive event.

The above experiments were carried out using pellets of transfected cells, as in previous studies (Bertocchini et al., 2004; Bertocchini and Stern, 2012; Shah et al., 1997; Skromne and Stern, 2002; Streit et al., 1998; Torlopp et al., 2014). One problem with this approach is that cells are likely to express other (unknown) factors that could influence the outcome of the signalling event. Another problem is that these pellets are relatively large (500-1000 cells). We therefore decided to substitute the use of cell pellets with protein-soaked microbeads (∼100 µm diameter). As neither VG1 nor NODAL are available as pure proteins, we decided to use ACTIVIN instead, which can induce axial structures and mesendodermal markers in chick epiblast (Mitrani et al., 1990; Stern et al., 1995). As shown in amphibian animal cap ectoderm explants (Green and Smith, 1990), ACTIVIN also acts through the SMAD2/3 pathway and generates finely graded responses of mesendoderm induction to different concentrations (Stern et al., 1995). BMP4-soaked beads were used as a source of inhibitory signal. First, we checked whether a single soaked bead can mimic the effects of a single cell pellet (Fig. S2). Grafting a bead soaked in ACTIVIN into the anterior marginal zone has the same effect as a cell pellet placed in the same position: it induces an ectopic cBRA-expressing primitive streak in adjacent epiblast (Fig. S2A-E). Conversely, placing a bead of the inhibitor BMP4 in the posterior marginal zone results in either displacement of the endogenous primitive streak to a more lateral position, or two primitive streaks, arising either side of the BMP4 bead (Fig. S2F-J). With a high concentration of BMP4 (50 ng/µl), primitive streak formation was abolished in about half of the embryos (Fig. S2J).

Next, we mirrored the experiments carried out with two or more cell pellets but using soaked beads (Fig. 2). After grafting a single ACTIVIN protein-soaked bead flanked by two control beads, 43% of embryos (6/14) had ectopic cBRA expression (Fig. 2B,G,K). When three ACTIVIN beads were grafted in a row to expose a wide domain to the inducing signal, the majority of embryos (78%, 7/9) showed no ectopic cBRA expression (Fig. 2C,H,K). When boundaries to the signalling domain were generated either by introducing a BSA-soaked control bead (Fig. 2D,I) or a BMP4-soaked bead (Fig. 2E,J) among the ACTIVIN beads, the proportion of embryos with ectopic cBRA expression was restored, to 40% (4/10) and 50% (6/12), respectively (Fig. 2K). In both cases, the proportion of embryos with two ectopic sites of cBRA expression was significantly higher than after grafting three ACTIVIN beads (P=0.0260 and 0.0507, respectively; Boschloo's test). Therefore, as with experiments using cell pellets, these results suggest that cells may sense inducing signals relative to their neighbours, rather than the absolute local amount of inducing signal.

To distinguish between the two alternative mechanisms of how cells might sense their positions (absolute local morphogen concentration or comparison of local signal strength in relation to their neighbourhood), two mathematical models were designed, one for each of these mechanisms, to make experimentally-testable predictions (for details see Materials and Methods). We model the marginal zone as a one-dimensional ring of cells (Fig. 3A). Positional information is provided by the balance between an inducer (SMAD2/3 activation in response to a VG1/ACTIVIN/NODAL-type signal) and an inhibitor (SMAD1/5/8 in response to a BMP signal) within each cell (Fig. 3B). Model A proposes that each cell independently assesses the concentration of morphogens (inducer versus inhibitor) it receives: when a threshold is exceeded, the cell is triggered to start primitive streak formation. Model B proposes that cells communicate with their neighbours to assess how the streak-inducing signal changes in space: each cell in the ring compares itself with the average signal strength in its neighbourhood to determine whether or not to initiate streak formation (Fig. 3B).

As an initial test of the model comparison method, we asked whether there exist parameter values allowing both models to replicate the experimental results shown in Fig. 2. We automated the search for parameter values using Bayesian computation, which scores values with a ‘likelihood function’ (Fig. S3). This function quantifies how well the predicted number and position of ectopic streaks match experimental results on a cell-by-cell basis. As an additional check of the automated parameter search, all parameter values found were tested for their ability to predict the initiation site of ectopic primitive streaks. If a set of parameter values allows for the model prediction to match the target result (here, the experimental results in Fig. 2), we define this set of parameters as ‘successful’. ‘Success’ is defined for each embryonic manipulation, as each manipulation results in different expected number and locations of ectopic streaks. The success rate of a model is defined as the proportion of predictions that are successful over the five manipulations shown in Fig. 3C-G. Although many parameter values yielded the same model success rates (Fig. S4), the likelihood function allowed further discrimination. Fig. 3C-G illustrates the output of each simulation when run with the single set of parameter values that provides both the greatest success rate for each model and the highest likelihood score. Even when the parameters were chosen to favour Model A, no set of parameter values was found that allowed Model A to replicate both experimental results in Fig. 3C and D simultaneously (the consequences of placing one or three beads of inducer in the anterior marginal zone). In contrast, Model B successfully predicts that broadening the domain of ectopic inducer reduces the chance of initiating ectopic streak formation (Fig. 3D).

The two models also differ in their ability to portray the effects of placing a bead of inhibitor between two beads of inducer (Fig. 3E-G). Model A predicts that the presence of the inhibitor will reduce the likelihood of ectopic streaks (Fig. 3E,F). However, Model B correctly predicts that only a low dose of inhibitor increases the chance of forming an ectopic streak (Fig. 3F,G). The same results were obtained irrespective of whether the sources of inducer and inhibitor were of smaller diameter (Fig. 3C-G, to simulate microbeads as in Fig. 2) or of larger size (Fig. S5, simulating a cell pellet as in Fig. 1).

We sought a single set of bead parameters that would allow both models to mimic the experimental findings (Fig. 3H). However, choosing a single set of bead parameters could act as a constraint, giving an advantage to one of the models. Therefore, we also performed the parameter inference to allow bead parameters to vary for each model independently (Fig. 3I). It is important the model success rates are the same, regardless of whether a single set of parameters is chosen to fit both models, or whether parameter values are optimised for each model separately (Fig. 3J) – strikingly, Model B always outperforms Model A.

In both models, cells measure their position by assessing the relative strength of the intracellular downstream effectors of the inducers (VG1/NODAL/ACTIVIN) and inhibitors (BMP). Therefore, decreasing the streak-inhibiting signal alone should induce ectopic primitive streak formation. In this case, both models predict this outcome (Fig. 4A,B).

To test these predictions experimentally, we used dorsomorphin, an inhibitor of BMP signalling (Yu et al., 2008). A dorsomorphin-soaked bead was grafted in the anterior marginal zone (Fig. 4A). After overnight culture, an ectopic primitive streak (with cBRA expression) was seen to arise close to the bead (Fig. 4C,D; P=0.0195; Boschloo's test). This result is consistent with a previous study showing that a graft of a cell pellet expressing the BMP antagonist CHORDIN in the area pellucida induces an ectopic streak (Streit et al., 1998). When embryos that had been grafted with a dorsomorphin-bead were examined 6 h after the graft, ectopic expression of cVG1 mRNA in the area pellucida (cVG1 expression is an early target of VG1/NODAL signalling; Skromne and Stern, 2002; Torlopp et al., 2014) was found in the vicinity of the bead (Fig. 4E,F; P=0.0005; Boschloo's test).

A more counterintuitive prediction arises when the strength of inhibition by BMP is increased in a region that normally expresses high levels of BMP (Fig. 5A). The two models predict different outcomes: Model A predicts that increasing BMP signalling in the anterior marginal zone will reduce the chance of ectopic streak formation (middle, Fig. 5B), whereas Model B predicts that introducing a bead of inhibitor will increase the streak-inducing values in an area adjacent to the bead (bottom, Fig. 5B). However Model B also suggests that this effect will be small, perhaps insufficient to result in formation of a mature ectopic primitive streak.

In embryological experiments in which a BMP4 bead was grafted into the anterior marginal zone, no cBRA expression or streak formation was observed after overnight incubation (Fig. 5C). After short incubation (4.5 h), however, cVG1 expression was observed in cells surrounding the grafted BMP4 bead in the anterior marginal zone and slightly in the adjacent area pellucida (Fig. 5D). cVG1-expression was absent from cells directly overlying the bead (Fig. 5F) (see also Arias et al., 2017). In addition, the ectopic expression was very weak, only detectable after prolonged chromogenic development of the in situ hybridisation (Fig. 5D,F). This ectopic expression of cVG1 in the anterior marginal zone was transient: it was seen at 4.5 h and disappeared by 6 h, remaining mostly in the lower layer of the area opaca (germ wall; Fig. 5E,G). In conclusion, this experimental result conforms with the predictions of Model B but not those of Model A.

We have seen that an increase in streak-inhibiting signal can result in paradoxical induction of cVG1, which is only predicted by Model B. However, no ectopic cBRA expression is observed. If it is indeed the case that cells assess their position in comparison with their neighbours (Model B), rather than measuring the absolute local levels of inducer and inhibitor, then introducing a sub-threshold amount of inducer flanked by low amounts of inhibitor would both deepen and steepen the gradient and might therefore be expected, perhaps paradoxically, to generate a new streak. Model A, in contrast, might predict that neither concentration is high enough locally to affect cell fates, resulting in a failure of ectopic streak formation. To simulate this, we explored parameter values for both models that could generate this result (Fig. 6). We found that only Model B can predict the initiation of an ectopic streak (Fig. 6D-F). No parameters were found that allowed Model A to produce the same result (Fig. 6D-F).

Next, we tested this prediction experimentally. We began by establishing the minimum threshold of ACTIVIN concentration for primitive streak induction; 2.5 ng/μl of ACTIVIN does not induce cBRA (Fig. S2D). When two BMP4-beads (6.25 ng/μl) were separated by a control bead, no ectopic primitive streak formed (0/9) (Fig. 6A,G). When an ACTIVIN-bead (2.5 ng/μl) was flanked by control beads (‘C-A-C’), 97% of embryos showed no ectopic primitive streak (n=37) (Fig. 6B,H). We then tested the predictions of the model experimentally: when a sub-threshold ACTIVIN bead was flanked by BMP4 beads (‘B-A-B’), cBRA expression was seen in 12.5% of cases (n=56; higher than ‘C-A-C’; P=0.0678; Boschloo's test) (Fig. 6C,I). However, a higher concentration of BMP4 (12.5 ng/μl) in the neighbouring beads reduced the proportion of embryos with an ectopic site of cBRA expression (to 9%; n=22), suggesting that at this concentration the total amount of inhibitor may overcome the small amount of inducer emitted by the sub-threshold ACTIVIN-bead. In conclusion, therefore, only Model B correctly predicts the counterintuitive results of this experiment.

Taken together (Fig. 7) our results favour a model by which cells assess their status (in terms of whether or not they will constitute a primitive-streak-initiating centre) in relation to the relative amounts of inducing and inhibiting signals they experience and also in relation to the status of their neighbours, rather than by direct readout of the local concentration of a morphogen that diffuses freely across the entire embryo.

Here, we propose a ‘neighbourhood watch’ model to explain how cells interpret positional information to determine the site of gastrulation. Our present results, both from computational modelling and experiments, favour the idea that cells do not read the concentrations of inducer and inhibitor (‘SMAD-value’) locally and cell autonomously, but rather interpret their own SMAD-value in relation to that of their neighbours. Moreover, the results suggest that the distance over which such comparisons take place is greater than just the immediately neighbouring cell on either side. In our ‘neighbourhood watch’ model, a neighbourhood size of 100-130 cells is predicted to satisfy experimental observations, based upon the parameter values estimated by the Bayesian inference algorithm.

In previous studies, multiple mechanisms have been uncovered by which cells interpret morphogen gradients. Can these other mechanisms explain our results? A key check when answering this question is to ask whether an alternative mechanism can explain the lack of ectopic streak and cBRA expression when an inducing signal is applied ectopically as a broad domain (Fig. 2). The first possible mechanism is that cells respond directly to morphogen concentration in a graded manner, as studied in explants of Xenopus embryos with a bead graft (Gurdon et al., 1995). Another study using cultured blastula cells not only supports this but also suggests that interaction with neighbouring cells is not required for the interpretation of morphogen concentration (Gurdon et al., 1999). However, this mechanism cannot explain our result of why a broad domain of inducer paradoxically reduces ectopic cBRA expression. A second possible mechanism of morphogen interpretation is one in which cells transform the signal concentration into the intracellular activity of a transcription factor, generating dynamic gene expression patterns with regulatory networks as shown for neural tube patterning (Cohen et al., 2013). Although this mechanism explains well the precision of different thresholds for interpreting morphogen concentrations based on duration and level (strength) of signals, it cannot explain our experimental observations, especially because we find that a broad domain exposed to inducing signal, without changing the duration of signals, reduced the incidence of ectopic cBRA expression. These considerations make it more likely that interactions between neighbouring cells are needed to position the primitive streak. A recent paper proposes that a neighbourhood comparison of signal strength [called the ‘spatial fold change (SFC)’ model] is required to position the determination front to regulate somite size in the zebrafish trunk and tail bud (Simsek and Özbudak, 2018), another example of a large developing field undergoing patterning. This suggests that a mechanism involving neighbourhood comparison for the interpretation of positional information may be used by different systems, especially if they are of large size.

In the ‘neighbourhood watch’ model in this study, as well as in the SFC model (Simsek and Özbudak, 2018), cells adopt a relative or normalised value to be evaluated, rather than the absolute morphogen concentration to assess their position. A relative value can provide a stable response of cells to signals, promoting robustness and precision in signal interpretation. Interestingly, a recent in vitro study suggests that cells sense relative signal intensity in the TGFβ/SMAD pathway as a fold-change value relative to the background to compensate for cellular noise (Frick et al., 2017).

How do cells communicate with their neighbours? In other words, by what mechanism could cells assess their environment? In the wing imaginal disc of Drosophila embryos, the TGFβ-related protein Decapentaplegic (Dpp) acts as a morphogen, conveying positional information that results in positioning the wing veins and other features of the wing. Signal-receiving cells have been shown to extend thin and long filopodia, called cytonemes, which extend several cell diameters to the proximity of Dpp-producing cells (Miller et al., 1995; Ramírez-Weber and Kornberg, 1999; Roy et al., 2011). It is worth noting that the existence of filopodia extending very large distances (connecting the invaginating archenteron with the future oral ectoderm at the opposite end of the embryo) was observed by Gustafson and Wolpert in studies of gastrulation in the sea urchin in 1961 (Gustafson and Wolpert, 1961) – this was one of the studies that initiated thinking on pattern formation. Similar structures have been observed in chick embryos during somite development (Sagar et al., 2015) but have not yet been sought at earlier stages. Another important question is: by what mechanism do cells sense relative signals compared with their neighbours? In our simulations, we mimic how each cell encodes the relative strength of inductive (SMAD2/3 activation by VG1/NODAL/ACTIVIN signals) and inhibiting (SMAD1/5/8 activation by BMP) cues they receive as the ratio between them. This is based on the proposal (Candia et al., 1997) that these two classes of SMADs (SMAD2/3 versus SMAD1/5/8) compete for binding to the ‘co-SMAD’, SMAD4. This could take place in both models – one possible mechanism to provide information about the status of neighbouring cells could involve hypothetical intermediate messengers conveying information about this state. Neighbourhood information could also be transmitted via a positive-feedback mechanism, for example a cell sensing higher levels of BMP would be stimulated to produce more BMP protein (Jones et al., 1992; Metz et al., 1998; Re'em-Kalma et al., 1995; Schulte-Merker et al., 1997).

One question is whether the mechanism proposed here (involving only local cell interactions and no long-range diffusion) is a feature unique to very large fields (several mm), where meaningful positional information conveyed by diffusion alone is likely to be impossible (Crick, 1970). There do appear to be several examples where diffusion of informative morphogens is key, such as initial patterning of the Drosophila blastoderm (Driever and Nüsslein-Volhard, 1988a; Gregor et al., 2007b) and mesoderm induction by activin in Xenopus animal caps (Gurdon et al., 1994, 1995; McDowell et al., 1997). However, in the chick embryo, the anterior-posterior distance between the two extremes of this ring should span ∼300 cell lengths (in reality the marginal zone has a thickness of about 120 µm, corresponding to ∼10 cells – here, we represent it as being one-cell-thick). As argued by Crick, it appears to be unlikely that this geometry can support the formation or maintenance of long-range gradients of morphogens generated by free diffusion (Crick, 1970). It therefore appears likely that positional information can be imparted by a variety of different mechanisms, perhaps according to the size and characteristics of the field to be patterned. It will be interesting to perform experiments comparable to those in this paper in a system such as anterior-posterior patterning of the chick limb, which is also large at early stages (HH18-20) and involves a localised signalling region (the Zone of Polarizing Activity) (Riddle et al., 1993; Tickle et al., 1975).

Here, we propose that positional information (when interpreted by a collection of cells) defines the location of the signalling centre (NODAL-expressing) that initiates primitive streak formation (Bertocchini and Stern, 2002). Initiation of a streak can be seen as the event that defines embryonic polarity. Our experiments and the associated models were designed to ask questions about how cells within the marginal zone assess their positions around the circumference of this signalling region, and thereafter determine the site next to which (in the area pellucida) the primitive streak will start to form. However, it is important to realise that in the embryo, the downstream consequence of these processes is not only a spot of cBRA expression, but rather a true ‘streak’, gradually extending towards the centre of the embryo. It has been shown previously that this elongation involves a process of cell polarisation and intercalation affecting the same site in the area pellucida where cells receive the inducing signals from the marginal zone (and which itself expresses cVG1 and NODAL) (Rozbicki et al., 2015; Voiculescu et al., 2007, 2014). Here, we observe cases where cBRA is induced but this is not followed by formation of an elongated primitive streak. For example, this result is seen when three beads are placed in the anterior marginal zone (A-B-A). One possible reason for this is that the embryos were not incubated for long enough to allow the intercalation to take place, but it is also possible that signals other than cVG1 and inhibition of BMP are required. Indeed it appears that non-canonical (planar cell polarity) WNT signalling drives intercalation (Voiculescu et al., 2007) within the area pellucida. Whatever mechanisms operate in the normal embryo to determine the site of primitive streak formation must somehow coordinate these signalling events to generate the full structure.

Taken together, we provide evidence that in a large system with two opposing gradients, cells assess their position in the field by measuring their location based on the relative concentrations of the inducing (cVG1/NODAL) and inhibitory (BMP) signals, and this is refined by taking cues from their local environment to assess the rate of change of these signals locally. However, the gradients are unlikely to involve long-range diffusion of two morphogens. Regulation of their strength is likely to involve other mechanisms resulting in gradients of transcription and therefore rates of production of the factors.

Fertilised White Leghorn hen eggs (Henry Stewart, UK) were incubated for 2-4 h to obtain EGK X-XI embryos, which were then harvested in Pannett-Compton saline (Pannett and Compton, 1924). After setting up for modified New culture (New, 1955; Stern and Ireland, 1981), the cell pellets or beads were grafted as required for each experiment, and the embryos cultured for the desired length of time before fixation in formaldehyde. Whole-mount in situ hybridisation was conducted as previously described (Stern, 1998; Streit and Stern, 2001). The probes used were: cVG1 (Shah et al., 1997), cBRA (Kispert et al., 1995) and BMP4 (Liem et al., 1995). Stained embryos were imaged under an Olympus SZH10 stereomicroscope with a QImaging Retiga 2000R camera. Some embryos were sectioned at 10 μm using a Zeiss Microm Type HM315 microtome.

HEK293T cells were seeded at 5×10⁵ cells/well in a six-well dish and incubated for 2 days (or 1×10⁶ cells/well for transfection on the next day) at 37°C in a total of 2 ml 10% fetal bovine serum (FBS) DMEM (growth medium)/well. On the day of transfection, the growth medium was changed to 1 ml/well of 5% FBS DMEM (transfection medium) at least 30 min before transfection. Transfection was carried out using polyethylenimine (PEI) as reported previously (Papanayotou et al., 2013). Briefly, 3 μl PEI (1 mg/ml) was added for every 1 μg of DNA transfected, in a total volume of 150 (for 0.5-2 μg)-200 µl (for 3-6 μg) DMEM in a sterile Eppendorf. Then 2 µg DNA were transfected/well (containing 6 μl PEI/well). Expression plasmids were the previously described DMVg1 (myc-tagged chimeric Vg1 containing the pro-domain of Dorsalin; Shah et al., 1997), pMT23 (murine BMP4; Dickinson et al., 1990) and pCAβ-IRES-GFP (as a control). The latter was also used to estimate transfection efficiency. Transfection mixtures were vortexed and then left for 10 min at room temperature for the PEI/DNA to complex. The transfection mixture was then added dropwise to the confluent monolayers of cells and incubated overnight at 37°C. The next day cells were checked for transfection efficiency of the GFP plasmid; typically, efficiency ranged from 60-90%. Cells were washed three times with 1× PBS, trypsinised and resuspended in a total of 1.5 ml growth medium and put into a sterile Eppendorf. The cell concentration was estimated in a haemocytometer. A bulk cell suspension of the transfected cells was made in the growth medium, so that each drop contained 500 cells in a total of 20 μl growth medium. Hanging drops were formed by placing the 20 μl aliquots on the lid of a 6 cm cell culture dish, the bottom of which was filled with 5 ml of sterile PBS or water to create a humidified atmosphere. After placing several such aliquots well-spaced in a circle, the lid was inverted and placed over the bottom of the dish, creating a mini culture chamber, to allow the cells to coalesce into pellets without adhering to the plastic. Culture dishes were incubated for 36-48 h at 37°C for the formation of pellets ranging in size from 500-1000 cells and used for grafts as required.

Recombinant human BMP4 (R&D Systems, 312-BP) was delivered using Affigel Blue beads (Bio-Rad, 1537302); recombinant human ACTIVIN A (R&D Systems; 338-AC) was delivered using Heparin-Acrylic beads (Sigma-Aldrich, H5236) and dorsomorphin hydrochloride (Tocris, 3093) was loaded onto AG1X2-formate beads. In each case the beads were incubated overnight at 4°C in the desired concentration of protein or chemical. Beads were washed in Pannett-Compton saline immediately before use.

The marginal zone was modelled as a one-dimensional ring of cells, comprising 600 cells in total [based on the assumption that the embryo at this stage contains 20,000-50,000 cells (Bertocchini and Stern, 2012) and on electron microscopy data (Lee et al., 2020; Voiculescu et al., 2007) for estimates of cell size and the radius of the marginal zone]. Proxies for streak-inducer and -inhibitor concentrations were assigned to each cell i, represented as V_i and B_i, respectively, with i=1,…,600 (Fig. 3A).

Before the addition of beads, streak-inducer and -inhibitor levels were inferred from a combination of RNA-seq reads (Lee et al., 2020) and in situ hybridisation of cVG1 and BMP4 (Fig. S1A), respectively, at approximately stage EG&K XII. To mimic these patterns, we use a Gaussian function to model the inducer levels based on the observed strong expression of cVG1 posteriorly, whereas inhibitor levels are modelled with a parabolic function to reflect the shallow, anterior-to-posterior gradient of cBMP4 (Fig. S1B). The placement of a bead is modelled as having an additive (or subtractive) effect on local protein concentration. The added values are constant for the width of the bead, and then decrease exponentially in space. Therefore, placement of a bead invokes four parameters (Fig. S5A): the position of the centre of the bead, the width of the bead, the concentration of the bead (relating to magnitude of the added values, see Fig. S5B) and the rate of decay of the added compound in space (i.e. the ‘spread’ parameter of the exponential distribution, see Fig. S5C).

For each cell to make its decision to initiate streak formation, we define the relationship between the amounts of SMAD2/3 (as a proxy for amount of inducer received) and SMAD1/5/8 (as a proxy for amount of inhibitor received) within the cells. This is based on the fact that inducing TGFβ-related signals (VG1/ACTIVIN/NODAL) act by phosphorylation of SMAD2/3, whereas inhibitory TGFβ-related signals (BMPs) phosphorylate SMAD1/5/8 – cells have been proposed to evaluate the relative strength of signals through competition of binding of these two classes of SMADs to the ‘co-SMAD’, SMAD4 (Candia et al., 1997). Inducing and inhibitory SMADs compete to form complexes with a fixed limited amount of SMAD4. The inducer- and inhibitor-linked SMAD complexes then move to the nucleus and regulate expression of different target genes.

We define two models for how cells interpret the SMAD-value to make the decision to initiate a primitive streak.

Both Models A and B have as parameters a threshold value (α or β) and protein concentration scalings (a_V and a_B). In addition, Model B requires the size of the neighbourhood (n) to be defined as a parameter.

For the final stage of the modelling process, we asked whether there exists a set of parameters allowing each model to replicate a target result. As both models invoke many parameters, resulting in a large and high-dimensional parameter space, we automated the search with a MCMC Bayesian computation algorithm. Parameter values were scored using a likelihood function, which quantifies how well model predictions match a target result. The target result was defined based upon an experimental result (Fig. 3 and Fig. S5) or a new possible theory (Figs 4–6).

For the parameter search, we fixed the expected width of the streak initiating domain, as well as the positions and widths of the beads. We allowed the concentration and spread parameters of the beads to vary (denoted c and s) in addition to all model parameters (α, β, a_V, a_B, n). Uniform prior distributions were defined for all parameters except the protein concentration scalings, a_V and a_B. For these parameters we defined b_V=log₁₀a_V and b_B=log₁₀a_B, which were then uniformly distributed. We defined biologically plausible ranges within which parameters were allowed to vary (shown in Fig. S7).

The posterior distributions of the parameters were obtained via the MCMC Bayesian computation in the pyDREAM package (Shockley et al., 2018), which implements a DREAM_(ZS) algorithm (Laloy and Vrugt, 2012). The algorithm was run using five Markov chains for a minimum of 5000 iterations per chain, and convergence was tested using the Gelman–Rubin statistic (Gelman and Rubin, 1992; Brooks and Gelman, 1998). The posterior distributions are shown in Fig. S7. An approximate neighbourhood size can be inferred from the posterior distribution of the parameter n (defined in Eqn 4), which peaked between 50-65 cells for all experiments.

The Bayesian computation algorithm maximises the likelihood (Eqn 8), quickly and efficiently finding sets of parameter values minimizing the distance between the target result (D_i) and the model result (f_i). Specifically, the likelihood function is defined so as to strongly favour sets of parameters where D_i and f_i have the same sign (i.e. both above zero or both below zero). Occasionally this means that parameter values obtained by the algorithm gave model values close to, but not exceeding, the threshold and therefore did not predict ectopic streaks as required by the target result. Therefore, all parameter values found using the Bayesian computation algorithm were checked to ensure that ectopic streaks were predicted in locations dictated by the target result. This was done by verifying that at least one cell exceeded the threshold to produce an ectopic streak in the expected location (i.e. the location of a bead). Thus, if parameter values for a given model allowed the prediction of correct ectopic streak placement, these values were deemed to give ‘success’ for a specific experimental design. The parameter values used in the plots in Figs 3–6 and Fig. S3 were chosen to maximise both the success rate and the likelihood. We have verified that there is a positive correlation between the success rate and the likelihood score (Fig. S4). All parameter values are given in Table S1.

The parameter search was performed for each group of experimental designs comprising Figs 1, 2, 4/5 and 6. Ideally, the parameter search must output a single set of bead parameters, allowing both models to approximate the target results as closely as possible (Fig. 3H). However, this acts as a restriction that might limit the ability of either model to replicate the target result. Therefore, the parameter search was also performed with all parameters varying for both models independently, removing this restriction (Fig. 3I). We verified that seeking a single set of bead parameters did not reduce the ability of either model to replicate the target result (Fig. 3J).

All statistical tests were performed in the R programming environment. The one-sided Boschloo's test was used.