### General procedure

We combined various explorative principles in order to maximise insight from the computer simulations. Factorial statistical design plans ensured that the parameter space was spanned systematically. Each solution image was characterised by a range of dedicated, quantitative descriptors, and the results interpreted by multivariate data modelling.

We analysed the high-dimensional, non-linear Delta-Notch model on a 2D lattice of hexagonal cells by this approach. For each set of parameter values the system of differential equations was integrated in time until a stable state was reached. The spatial distribution of the stable-state Notch concentrations were plotted as black and white 2D images of the cell lattice for easy visual inspection.

The exploratory process proceeded through four distinct stages (Figure 3). Together, these four stages were intended to yield maximum new insight into the behavioural repertoire of the model with minimum effort and minimum number of additional assumptions.

(1) The interesting ranges for the model parameters were determined in a preliminary exploration. In this trial-and-error based phase the parameters were changed one at a time, and the solution images inspected informally, to determine values to be used as "high" and "low" for each parameter in the subsequent designs.

(2) Secondly we defined a 2^{7} full-factorial design and expansions of this in order to explore the parameter space domain found interesting in phase (1). The equilibrium solutions images were submitted to computerised image analysis by standard mathematical filters of various kinds, like gray-tone density and spatial frequency histograms. This profile of many, but *per se* meaningless image descriptors could be related to the seven known model parameter values by well-established multivariate methods from chemometrics, primarily conventional reduced-rank Partial Least Squares Regression (PLSR) [10] with optimal rank determined by cross-validation as described in [11, 12]. This permitted us to identify parameter combinations which could be ignored in the main explorative experiment because they had little or no effect on the solution patterns.

(3) In the main experiment we employed a more informative, but also more laborious solution profiling. It consisted of a fractional factorial design created in order to span the range of all parameters and their combinations that were found to have a clear effect in (2), by means of a minimum number of numerical simulation solutions. Due to the insight obtained from our previous screening, our new design could be reduced to only 32 selected combinations of a high or low value of each of the seven model parameters, with six non-informative (no patterning) solutions replaced by two replicates of an intermediate parameter combination with independent random starting value sets, each assessed independently in three sensory parallels (see Additional file 1).

Human sensory descriptive analysis was used for image description. It provides systematic, inter-subjective description in intelligible terms with high repeatability. Such sensory descriptive analysis is standard procedure in food science [13] and also used for describing and comparing e.g. perfumes or the quality of MRI images [14]. Here we used it for assessing printed images of outcomes from the model of cell differentiation. A panel of professional human assessors first developed adequate terminology, and was then trained in describing and comparing a set of images of steady-state solutions from the simulations with respect to these descriptor terms [13, 15]. During this process each assessor's subjective "private language" was replaced by a common "inter-subjective" set of terms which was agreed upon by all persons involved, both with respect to what it means qualitatively and to how it should be reported quantitatively. The sensory analyses are detailed elsewhere [16].

The resulting tables of profile data were modelled by multivariate "soft modelling" based on reduced-rank Partial Least Squares (PLS) regression [10–12]. The main covariation patterns among the observed descriptors were identified and related to the chosen model parameter values in a subspace-model defined by statistically significant covariance eigenvectors. Cross-validation was used to distinguish between phenomena with predictive validity and apparent noise effects. This made it possible to find and interpret the main patterns of co-variation while ignoring minor details that otherwise would have lead to over-parametrisation.

The 32 image profiles could thereby be mapped onto points in the PD space, and PLSR was used for further analysis of the descriptors and their relation to the model parameters. Figure 2 shows examples of pattern classes thus identified.

(4) Finally, we pursued unexpected pattern types discovered in the main experiment (3) in more detail in a final follow-up experiment, involving a dense sampling of the model parameter *θ*_{D} for a few fixed combinations of the other parameters. The solutions were characterised by both sensory descriptive analysis and by a computational method by which the frequency of the observed Notch concentration levels were recorded at each parameter value. (See Additonal file 1 for documentation of the research process according to Figure 3.) In the following we summarise the main results.

### Pattern Descriptor Space and Pattern Descriptor Model

The parameter ranges found to be of interest in the preliminary exploration (1) (see Additonal file 1: Table S-T2) were used for designing the subsequent two experiments (2) and (3). The image analysis filters employed in the extensive screening experiment (2) (see Additonal file 1: Table S-T3) revealed systematic co-variation patterns in the image analysis profile data. But these patterns were difficult to interpret. The main, sensory-based experiment (3) employed a modified version of the reduced design (see Additonal file 1: Table S-T4), and consisted in using a trained sensory panel to define verbal pattern descriptors and then quantify the patterns in each image in terms of these descriptors. Based on input from the investigators the sensory panel defined 12 descriptors which they considered sufficient and by which they subsequently quantified the patterns (see Additonal file 1: Table S-T5). This allowed us to develop a PD model (Figure 4) by bi-linear PLSR regression in latent variables, relating the chosen 12 sensory descriptors (see Additonal file 1: Table S-T5), to the seven model parameters via 3 latent variables, over the 32 images. Perturbations in the model parameters *θ*_{D}, *θ*_{N}, p_{D} and p_{N}, alone and in combinations, clearly affected the solutions in systematic ways, while variations in *μ* (the ratio between the degradation rates of Delta and Notch) and in the initial conditions had fairly small effects. Based on split-half cross-validation [12], the PD model correctly predicted 88% of the variance in the 12 sensory descriptors from the model parameters. Conversely, 85% of the variance in the selected non-linear model parameter combinations was correctly predicted from the pattern descriptors.

We found that the PD model clustered the 32 graphical patterns in four main classes, here named I, II, III, IV (see Additonal file 1: Table S-T6). In a cross-validation experiment in which each image was treated as unknown against a PLS discriminant analysis model estimated from the other 31 images, the set of pattern descriptors could be used for successful classification of each of the 32 images into the four classes. The clustering was confirmed by independent hierarchical cluster analysis of the patterns based on their Euclidean distances in PD space, and appears to be quite robust (see Additonal file 1: Figure S-F3).

Images belonging to one and the same class share obvious visual pattern characteristics that distinguish them from images from the other classes, and can be associated with different domains in parameter space (see Additonal file 1: Table S-T6, Figure S-F4). For example, low values of *θ*_{D} are a clear indication that the resulting pattern will fall into class I. The images in class I are likely to be high in the descriptors *Sharpness*, *Contrast* and *PatternBlack*, and low in *Whiteness* and *PatternWhite*, implying primarily black and white cells, the former outnumbering the latter and participating in the curls which dominate the pattern in this class (Figure 2A).

### Parameter-dependent differentiation phenomena

We subsequently (Phase 4) checked that the clustering is not just an artefact stemming from the low number of patterns in the fractional factorial design by running the fourth and final set of simulations to pursue details. We varied the Delta threshold parameter *θ*_{D} in small steps to traverse the PD space along a straight line, illustrated in the biplot of Figure 4a from the previous simulation experiment, from class I (*θ*_{D} = 0.1, characterised e.g. by sensory descriptor *Sharpness*) via class III (*θ*_{D} = 0.7, characterised by *Whiteness* and *MultiShade*) and up to *θ*_{D} = 0.9. Some new sensory descriptors were added (see Additonal file 1: Table S-T7).

From these results, good predictive PLS regression models were developed between the sensory descriptive profile and the model parameters within the parameter range calibrated for, as exemplified for *θ*_{D} in Figure 5. Moreover, rather than spreading evenly in PD space, the image properties jump between a few clearly separated domains in the 3-dimensional PLS component space, with few intermediates (Figure 6).

As *θ*_{D} increased, the estimated Notch activity passed through a series of complicated, abrupt changes resembling bifurcations (Figure 7a), showing that the statement of Collier *et al*. [7] that the patterns are insensitive to the precise values of the model parameters is only partially true. Ordinary bifurcation diagrams show the number of coexisting stable states of which the system may occupy just one at a time, but Figure 7a shows the number of stable Notch levels actually expressed simultaneously over the lattice for each value of *θ*_{D}. That is why this is not a bifurcation diagram in the common sense. We propose to call it a *differentiation diagram*, where *θ*_{D} is a *differentiation parameter* regulating the number of simultaneous levels, or the degree of differentiation. This being said, the diagram shows a large number of "bifurcation" values of *θ*_{D} for which the number of simultaneous Notch levels changes. The number of levels is related to the complexity of the image: a high number of levels come with a complex and highly irregular pattern. A high number of equilibrium levels, e.g. at either side of *θ*_{D} = 0.41 and 0.62, correlates with long integration time (Figure 7a). Both can be seen as a consequence of high *frustration* [9], which implies that the system needs extensive fine-tuning and either takes a long time to settle in a stable, patterned state or is chaotic. Despite the fact that no descriptor was designed with analysis of differentiation in mind, the descriptors characterise the patterns in a way that at least partly reflects the complicated differentiation pattern for varying *θ*_{D} (Figure 7b).

### Initial state variations, repeatability and boundary effects

For a given set of parameter values the total number of stable states in a large network is probably very large. However, each initial state (*N*^{0}, *D*^{0}) lies in the attractor basin of just one patterned state. An important question is whether the patterns generated by all (*N*^{0}, *D*^{0}) for a given set of parameter values lead to more or less the same state in PD space. For a realistic model of a biological system we would expect this to be true, as biological systems are almost by definition structurally stable and thus insensitive to internal and external stochastic fluctuations. Thus, contrary to stable points in phase space, we would expect points in PD space to have large attractor basins.

Once the clustering into four classes was discovered during the multivariate analysis of the sensory data, it was straightforward to test the robustness and the visual appearance of these pattern types. Varying initial conditions within the imposed limits (see Methods) changed the details of the patterns without appreciably altering their descriptor profiles. This is illustrated by the proximity of the two independent random initial condition replicates in class IV (Figure 4a, open triangles Δ, ▽). The precision of the sensory descriptive analysis is also illustrated in the same figure by the closeness of the three independent sensory parallels within each of these two initialisation replicates.

With periodic boundary conditions, which are commonly used to mimic an infinite domain, we expected that the boundary effects would be negligible as regards the pattern characteristics. This was confirmed by simulations in which the original 50 × 50 lattice was embedded in a larger (75 × 75) lattice which was run with the same set of random initial perturbations as in the 50 × 50 lattice. The final patterns in the 50 × 50 lattice in the two cases belonged to the same class, and in most cases were almost identical (data not shown). Some very large-scale spatial patterns appear as boundaries between domains with a three-periodic pattern, but with a phase shift along the boundaries (Figure 3B). In an infinite-sized lattice we expect these patterns would fail to appear. By increasing the pattern size to 51 × 51 we found that some of these changed to a 3-periodic pattern all over the lattice without phase-shift boundaries (51 is divisible by 3).