The genotype-phenotype relationship in multicellular pattern-generating models - the neglected role of pattern descriptors

Modelling phenotypic variation in biological pattern generation

The whole development process of higher organisms can be mathematically conceptualised as a recursive mapping - i.e. successive cell differentiations leading to a sequence of unfolding patterns at many different spatiotemporal scales, each pattern defining the context for further differentiation and thus for subsequent patterning processes. A deep understanding of what causes the phenotypic variation arising from such patterning processes cannot be claimed before we are able to recreate this variation theoretically by what we call causally cohesive genotype-phenotype models (cGP models) [1]. Unlike the broader class of mechanistic mathematical models describing how complex biological phenotypes arise from the interactions of lower-level systemic entities, cGP models are distinguished by linking together (cohering) the individual's genotype and its phenotype in a causal mathematical structure. cGP models thus allow the construction of genotype-phenotype maps, i.e. mappings predicting the phenotype associated with a given genotype based on what we know about the regulatory anatomy of a given biological system.

Irrespective of the level of biological resolution of a cGP model, genetic variation has to be represented as parametric variation. In a genotype-phenotype-map perspective one is thus interested in getting a clear understanding of the mappings between genotype parameter space and the generated phenotypic space. However, in the context of multicellular patterning models it is not trivial to establish this relation. The very concept of pattern implies that what matters is not the state of every single cell, but certain emergent qualities of the total cell aggregate which express relations between the states of subsets of the cells. Thus, in order to set up a genotype-phenotype map in such a spatiotemporal pattern setting one cannot just establish a mapping between domains in parameter space and certain properties of locally defined intracellular and extracellular state variables. One is actually forced to construct new descriptors of the emergent pattern features, express them in an abstract pattern descriptor space, and establish their relations to properties of the original model. In addition to leading to a genotype-phenotype map, this approach opens possibilities for additional validation of the model through prediction of higher level and empirically observable properties that are by no means part of the model's premise set, and which represent emergent features of the pattern that are interesting and relevant from the point of view of purpose and objectives of the original pattern-generating model.

There is today no general procedure for how to relate a set of patterns and their characteristic features to the functional relationships, parameter values and initial values of the original pattern-generating model. Here we propose a multivariate data modelling approach which is based on three major elements: (i) cost-effective computer simulations to probe the high-dimensional parameter space, (ii) informative ways to describe the model's graphical patterns quantitatively as points in a pattern descriptor (PD) space, and (iii) ways to establish two-way mappings between the PD space and the parameter space in terms of a statistically reliable and interpretable statistical prediction model (Figure 1).

We illustrate our approach on a simple mathematical model of pattern generation, and show that even in this case, traditional and intuitive methods prove inadequate to describe and understand its full potential of pattern formation, while our approach provides an overview and reveals unexpected features and their relationship to the model parameters.

The illustration model

A large number of cell differentiation processes involve the membrane-bound protein Notch which interacts with several transmembrane ligands in neighbouring cells [2–4]. Thus, understanding the diverse functions of Notch is of paramount importance in itself. One patterning mechanism in which Notch is involved is lateral inhibition: a cell-cell interaction whereby a cell developing a specific fate inhibits its neighbours from developing in the same way. Delta is one of the ligands binding to Notch; Delta in one cell binds to Notch in the cells in physical contact, leading to juxtacrine signalling [5]. The multicellular gene network regulating this process is called the neurogenic network [6]. However, as early as 1996 Collier et al. [7] presented what can be considered as a simplified multicellular cGP model for Delta-Notch activity in a one- or two-dimensional lattice of discrete cells. In this model the rate of change of Notch activity increases with the average activity level of its ligand Delta in neighbouring cells, while the rate of change of Delta activity decreases with increasing activity level of Notch in the same cell, both relationships being expressed in terms of sigmoidal response functions.

On a two-dimensional lattice of n identical hexagonal cells we let D _k and N _k , k = 1,..., n, represent the activity level of the two proteins Delta and Notch in cell number k. In suitably scaled variables, assuming first-order degradation, and with the standard Hill function representing sigmoid stimulus-responses, the dimensionless model of Collier et al. is

(1a)

(1b)

in which μ is the ratio between the degradation rates for Delta and Notch, is the average value of Delta in the six neighbours of cell k, and S is the standard Hill function S(x, θ, p) = x ^p /(x ^p + θ ^p ). Thus, θ_D and θ_N are the thresholds for the action of Delta and Notch, respectively, and p_D and p_N are the corresponding steepness parameters (Additonal file 1: Table S-T1).

The positive feedback loops between any pair of neighbouring cells which these interactions generate, lead to differentiation: high Delta and low Notch activity in a cell favours low Delta and high Notch activity in neighbouring cells, and vice versa. On a one-dimensional string of discrete cells, the final, differentiated states are strongly dominated by an alternating pattern: apart from scattered exceptions in which two neighbouring cells both express either Delta or Notch, every second cell expresses Delta and suppresses Notch, and vice versa. On a two-dimensional lattice of square cells this behaviour is consistent with the well-known checkerboard pattern. Using standard linearisation methods, Collier et al. [7] showed that the system of regular, hexagonal cells possesses three basic 3-periodic patterns consisting of up to three cell types. In our simulations the basic regular patterns appear either by cells expressing Notch surrounded by cells in which Notch expression is inhibited, or vice versa.

However, two-dimensional lattices of real cells are usually rather irregular and do not have periodic patterns. Podgorski et al. [8] have combined the model of Collier et al. with models for apoptosis and differential adhesion to study the patterning process on irregular cell lattices. But even on a regular lattice of hexagonal cells irregular patterns arise because the system is frustrated [9]: it is impossible to satisfy the requirement of opposite expression of Delta and Notch between all pairs of neighbouring cells. If cell 1 expresses Notch, its neighbouring cells 2 and 3 should both express Delta, which causes frustration because cells 2 and 3 also may be neighbours. This results in a multitude of different patterns, sometimes with a high number of different Notch and Delta levels, depending on parameter values and initial conditions (Figure 2).

In fact, our simulations showed that the majority of patterns and pattern features are highly irregular. Initial trial and error simulations for different combinations of parameter values revealed a wide range of different kinds of patterns, and for certain parameter value combinations we observed a large number of different equilibrium levels. Trying to describe these irregularities as combinations of basic regular patterns did not seem to make much sense. Rather, our approach was to try to introduce new concepts and pattern descriptors to classify and describe the irregular pattern as patterns in their own right. It was also obvious that without a systematic, planned simulation design as outlined above, there could be no hope of achieving a sound description and classification of the pattern and reveal their relation to variations in parameter values.

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⁷ 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⁰, D⁰) lies in the attractor basin of just one patterned state. An important question is whether the patterns generated by all (N⁰, D⁰) 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).

Model analysis

With suitable boundary conditions, for example periodic, the model in Eqs. (1) has a homogeneous equilibrium (D*, N*) given by the unique solution of D* = 1-S(N*, θ_N, p_N), N* = S(D*, θ_D, p_D). This state is unstable unless the steepness of the sigmoid functions is very gentle [7, 17], and a perturbation of the homogeneous state will initiate a transient during which cells develop towards different final states, creating a spatial pattern. It is far from obvious which parameter value combinations generate which patterns and how patterns could be described and classified.

The patterns were obtained by integrating Eqs. (1) numerically from initial values D _k ⁰, N _k ⁰ on a 50 × 50 or 51 × 51 lattice of hexagonal cells with periodic boundary conditions. The lattice dimensions 51 × 51 were chosen to avoid certain shifts in the basic three-periodicity entirely due to the periodic boundary conditions (50 is not divisible by 3). The initial values for the individual cells represented small, random perturbations from the homogeneous steady state value D*. To generate clearly different initial value sets, we modified the maximum size M and sign s of the perturbations (see Additonal file 1: Table S-T1), if s > 0 (< 0), all perturbed values were larger (smaller) than the equilibrium value. Thus, including our parameterisation of the initial values, our model has a total of seven parameters which span the Parameter Space.

Simulations and descriptive analysis of patterns

As a start we make some useful distinctions. For each choice of parameter values and initial values the model yields a stable solution, from which we generate a two-dimensional image in which each cell is given a shade of gray defined by its Notch level (white: N = 0, black: N = 1). The image exhibits certain pattern characteristics, in short, a certain pattern. A pattern is a conspicuous, emergent feature of an image (could be periodic or non-periodic). The pattern characteristics, e.g. Whiteness, are qualities. If the overall Notch level is low, we characterise the image by attaching the quality Whiteness to it. Using sensory methods or dedicated mathematical filters, the degree of whiteness is given by the value of the scalar variable Whiteness for this particular image.

The most interesting ranges for the individual model parameters were sought by trial-and-error in an initial exploration (1). Each parameter combination was tested with a new set of random perturbations for the simulation. A second, full factorial screening design 2⁷ with centre-points (2) was then chosen as a compromise between computational load and probability of finding important patterns (see Additonal file 1: Figure S-F1, Table S-T2). The multivariate profiling (see Additonal file 1: Table S-T3) and data analysis (cross-validated PLS regression of the relationships between the design parameters and about 80 automatically generated image analysis descriptors) showed that several kinds of pattern were indeed generated by the simulations (see Additonal file 1: Figure S-F2). Based on preliminary results, the design was expanded in a couple of stages to ensure that the relevant parameter domains were adequately spanned.

To ensure a more informative and interpretable characterisation, we decided to apply human visual evaluation of an interesting subset of the solutions in the main experiment (3). Since this is a more laborious and expensive measuring principle, a reduced statistical design was chosen for this third stage. From the original screening design, a fractional factorial design 2^7-2 was selected by combining a high or low value of each of the seven parameters, confounding the least interesting main effects with higher-order interactions of other factors (see Additonal file 1: Table S-T4). Six of the 32 chosen parameter combinations resulted in a return to the homogeneous equilibrium, and were replaced by three centre points in two replicates with different initial conditions.

For each of the 32 experimental conditions, the full 50 × 50 equilibrium lattice of Notch was printed on paper in black-and-white. The sensory descriptive analysis, performed by a panel of eleven assessors working as trained sensory judges at the Norwegian Food Research Institute, had two distinct steps, carried out according to the ISO method convention [13], as e.g. used in [14]. (i) The assessors and the researchers agreed on twelve descriptors reflecting characteristic pattern features (see Additonal file 1: Table S-T5). The assessors were then trained to respond similarly with respect to each descriptor. (ii) Each image was judged with respect to each of the twelve descriptors on a scale from 1 to 9 by the eleven individual assessors, with the images anonymised and presented in random order. Finally, PLSR with split-half cross-validation was used in various ways to find reproducible patterns of co-variation between the twelve sensory descriptors and the seven design factors and their interactions (Figure 4). Graphical inspection of the PLS score plot (Figure 4a) as well as cluster analysis led us to group the solutions into four distinct classes (see Additonal file 1: Figure S-F3, Table S-T6).

In order to pursue in detail the effect of changing the parameter θ_D from the low values giving class I through class III (θ_D ≈ 0.7) up to θ_D = 0.9, a final sensory analysis experiment (4) was performed on simulation results for 36 values of θ_D and constant values of the other parameters, this time using a 51 × 51 lattice and a fixed set of randomly chosen initial values (Figure 7). To focus on the local fine-structure of the images, the panel and the project leaders agreed upon a slightly revised set of twelve descriptors of which two are reported here (see Additonal file 1: Table S-T7). The eleven assessors and the general evaluation conditions were the same as before.