Continuum limits of pattern formation in hexagonal-cell monolayers

Intercellular signalling is key in determining cell fate. In closely packed tissues such as epithelia, juxtacrine signalling is thought to be a mechanism for the generation of fine-grained spatial patterns in cell differentiation commonly observed in early development. Theoretical studies of such signalling processes have shown that negative feedback between receptor activation and ligand production is a robust mechanism for fine-grained pattern generation and that cell shape is an important factor in the resulting pattern type. It has previously been assumed that such patterns can be analysed only with discrete models since significant variation occurs over a lengthscale concomitant with an individual cell; however, considering a generic juxtacrine signalling model in square cells, in O’Dea and King (Math Biosci 231(2):172–185 2011), a systematic method for the derivation of a continuum model capturing such phenomena due to variations in a model parameter associated with signalling feedback strength was presented. Here, we extend this work to derive continuum models of the more complex fine-grained patterning in hexagonal cells, constructing individual models for the generation of patterns from the homogeneous state and for the transition between patterning modes. In addition, by considering patterning behaviour under the influence of simultaneous variation of feedback parameters, we construct a more general continuum representation, capturing the emergence of the patterning bifurcation structure. Comparison with the steady-state and dynamic behaviour of the underlying discrete system is made; in particular, we consider pattern-generating travelling waves and the competition between various stable patterning modes, through which we highlight an important deficiency in the ability of continuum representations to accommodate certain dynamics associated with discrete systems.

waves and the competition between various stable patterning modes, through which we highlight an important deficiency in the ability of continuum representations to accommodate certain dynamics associated with discrete systems.

Introduction
The derivation of continuum models which represent underlying discrete phenomena is emerging as an important part of mathematical biology: integration between subcellular, cellular and tissue-level behaviour is crucial to understanding tissue growth and mechanics with self-evident application to, for instance, in vitro tissue engineering or the understanding of tumour growth and invasion.
It is well known that cell-signalling mechanisms regulate differentiation, cell-fate determination and, ultimately, tissue and organ development. Such regulation is mediated by the production, transport and binding of intercellular signalling molecules, which may be free to diffuse throughout the tissue or may be anchored in the cell membrane. In the latter case, cell-signalling molecules may bind only to directly adjacent cells; such a juxtacrine signalling mechanism is therefore of particular interest in closely packed structures such as epithelia.
Lateral inhibition is a juxtacrine pattern-forming mechanism employed by developing tissues to create fine-grained patterns of cell differentiation, in which adjacent or nearby cells diverge to achieve differing cell fates. This mechanism is controlled by a negative feedback loop: receipt of inhibition reduces the ability of a cell to inhibit others, leading to the amplification of differences between cells. This mechanism is evolutionarily conserved and is observed in insects, nematode worms and vertebrates, in all of which the transmembrane proteins Notch and Delta (or their homologues) have been identified as mediators of the interaction (see Collier et al. 1996 and biological references therein). Well studied examples of such fine-grained patterns include compound eye and nervous tissue development in insects, such as the fruit fly, Drosophila (Haddon 1998;Appel et al. 2001;Carthew 2007). Within the context of regenerative medicine, Delta-Notch signalling has been shown to regulate cell fate in stem cell clusters (Lowell et al. 2000). While other ligand-receptor interactionmediated cell signalling mechanisms have been characterised (e.g. the binding of cyclic AMP to Dictyostelium cells (Martiel and Goldbeter 1987;Dallon and Othmer 1997) and Transforming Growth Factor-α and Epidermal Growth Factor binding in keratinocytes (Clark et al. 1985;Coffey et al. 1987), here we consider the well-studied Delta-Notch signalling interaction, which provides an ideal model system to illustrate our methodology.
Inspired by the microscopic, contact-dependent nature of the juxtacrine signalling mechanism and the short-range patterns often observed in tissue development, many authors have used a discrete mathematical formulation to investigate the pattern-forming potential of such signalling mechanisms. The first such study was presented by Collier et al. (1996), in which Delta-Notch binding was considered. Lateral inhibition was shown robustly to produce fine-grained patterns of the kind observed in early development, provided that the feedback strength was sufficiently strong. The model was formulated in terms of ordinary differential equations (ODEs) representing Delta and Notch activity on individual cells. Many other relatively recent studies have considered a discrete representation of juxtacrine signalling. Owen and Sherratt (1998) analysed a more complicated model, considering explicitly the numbers of ligand and free and bound receptors on each cell. Lateral induction (positive feedback between ligand-receptor binding and subsequent ligand production) was accommodated and the range over which juxtacrine signals may be transmitted was studied; Wearing and Sherratt (2001) performed a comprehensive nonlinear analysis of this model, highlighting that linear analysis alone is unable to predict the model's behaviour for large numbers of cells. Webb and Owen (2004) extended this model, considering the dynamics of ligand and free and bound receptors in systems of varying geometry (strings and square or hexagonal arrays), showing that lateral inhibition can produce patterns with a lengthscale of many cell diameters and that cell shape is a crucial determining factor in the patterns produced.
Such discrete approaches necessitate intensive numerical study, especially when linear analysis cannot predict the patterns formed. Extensive study of pattern formation has also been undertaken within continuum formulations. Continuous reactiondiffusion models have formed the basis of many models of biological pattern formation since the study of Turing (1952) in which it was shown that reaction and diffusion of chemicals can produce heterogeneous distributions of chemical concentration that consequently determine cell fate (we remark that Turing (1952) employed both discrete and continuous analyses, however). In these models, spatial patterns of morphogens are assumed to induce cells to differentiate. Examples include Kauffman et al. (1978), in which segmentation of Drosophila was considered, and Varea et al. (1997), who considered a Turing system on a growing domain to model the formation of skin patterns in fish. An alternative continuum approach is known as mechanochemical modelling in which the patterns in biological tissue are dictated by mechanical laws applied to cells and their environment, reflecting, for instance, tissue deformation or cell migration (see, for example, Murray et al. 1988). Here, pattern formation and morphogenesis take place simultaneously and the system may therefore adjust to external disturbances, an important feature of embryonic pattern formation.
Discrete models are able to reflect the inherently discrete nature of cell population behaviour, capturing explicitly the interactions between individual cells, cell movement or short-range patterning. Appropriate continuum models of such phenomena facilitate their incorporation into tissue-scale modelling and, in addition, may admit analytic progress or simpler numerical analysis. For these reasons, multiscale (or homogenisation) techniques have been employed to derive continuum models directly from underlying discrete systems, enabling some of these discrete effects to be incorporated into tissue-scale models in a mathematically precise way. The method of multiple scale expansions for partial differential equations is well-developed (see, for example, Kevorkian and Cole 1996) and widely used to derive models for a variety of physical and biological problems. In a biological context, such techniques have been employed by (e.g.) Turner et al. (2004) and Fozard et al. (2009) to represent, within a continuum formulation, the collective motion of adherent epithelial cells.
It has typically been assumed that the analysis of cell signalling mechanisms demands a discrete approach, especially when considering contact-dependent juxtacrine signalling processes (Wearing et al. 2000;Plahte and Øyehaug 2007) and short-range patterning (Roussel and Roussel 2004); however, in our previous work (O'Dea and King 2011), we have shown how a multiscale method may be employed to analyse the fine-grained patterns of period two (in each coordinate direction) generated by the discrete model of Collier et al. (1996) in an array of square cells within a continuum formulation. The resulting reduced model, which takes the form of a semilinear partial differential equation (PDE), demonstrates that on the macroscale the interaction between adjacent cells manifests itself as a diffusive process despite the short-range variation in signal concentration. Travelling wave analysis showed excellent quantitative agreement between the continuum formulation and the underlying discrete model.
Square cells were considered in O'Dea and King (2011) as an initial geometry with which to illustrate the methodology. However, Webb and Owen (2004) have demonstrated that cell shape is a determining factor in the patterns produced by cell signalling systems; furthermore, relevant biological structures such as the simple squamous epithelium are not necessarily well-approximated by square cells, and can display hexagonal morphology of striking regularity (e.g. squamous epithelial cells; the hexagonal configuration is, of course, one of minimal surface energy, as seen in soap films, for example). An example of such cell shapes is given in Fig. 1. Inspired by this, in this paper, we extend our previous work to consider the fine-grained patterning behaviour in an array of hexagonal cells. Such investigations are therefore biologically relevant as well as revealing important mathematical insight; while in the cell signalling systems under consideration, cells may not form exactly hexagonal arrays, our caricature enables the derivation of a continuum model amenable to analysis which, nevertheless, reflects microscale complexity of biological relevance. We note further that, Webb and Owen (2004) show that small random perturbations of cell geometry do not prohibit the formation of regular patterns; we therefore expect the qualitative features of our results to apply to the non-uniform case.
In one-dimensional strings (or arrays of square cells), fine-grained patterning corresponds to patterns with a period of two cell lengths (in each coordinate direction), whilst hexagonal cells imply patterns of period three in each direction; though still fine-grained in nature, the increased number of possible patterning modes results in a significant increase in mathematical complexity. Concentrating on patterns of period three cell lengths, we construct individual continuum models for the generation of patterns from the homogeneous state and the transition between the various patterning modes in response to variation of a parameter associated with feedback strength. In addition, via a two-parameter expansion, we construct a more general continuum representation of this system, capturing the emergence of the multiple bifurcation structure from a single pitchfork-like bifurcation. The resulting continuum models allow representation of the generation of microscale patterns in a tissue in response to macroscale variation in cell signalling (e.g. that induced by tissue-level chemical or physical stimulation). Via comparison with both the steady-state and the dynamic behaviour of the underlying discrete system, we show that the continuum models faithfully represent the fine-grained patterning behaviour displayed by the discrete signalling model. Quantitative comparison is made by analysing the travelling-wave behaviour; specifically, the speed of a 'linearly-selected' pattern-generating wave invading an unstable patterned state is considered, providing an indication of the range of applicability of our continuum formulations. In addition, we consider in detail the competition between various stable patterning modes, revealing important insight into the ability of continuum representations to accommodate such discrete dynamics.
In Sect. 2, the model of Collier et al. (1996) is summarised and its period three patterning behaviour examined. In Sect. 3, a continuum formulation for period three pattern generation by Delta-Notch signalling in hexagonal cells based on this underlying discrete system is derived. Solution of this model is compared to numerical simulations of the discrete system in Sect. 4, illustrating the generation of patterns under spatial variation of feedback strength, and pattern competition and travelling-wave behaviour in the case of uniform feedback. In Sect. 5 our findings are summarised, together with suggestions for future research.

Formulation
In Collier et al. (1996), the feedback between the binding of a membrane-bound signalling protein, Delta, to its receptor, Notch, and subsequent Delta expression was considered. The crucial aspect of the feedback loop is that elevated Delta expression in a cell downregulates Delta expression in its neighbours via the receptor, Notch. This mechanism, known as lateral inhibition, is a fundamental cell-fate control mechanism (Mitsiadis et al. 1999), creating fine-grained patterns in developing tissues which determine subsequent cell development. The key postulate of the Delta-Notch signalling model is that the level of activated Notch in a cell determines its fate: low levels lead to the adoption of default (primary) fate, whilst high levels relegate the cell to the secondary fate. In the specific case of nervous tissue in Drosophila, the primary fate corresponds to the adoption of a neural phenotype, the secondary fate being the maintenance of the epidermal phenotype (Lehmann et al. 1983;Campos-Ortega 1993). The model is formulated in terms of Delta and Notch "activity"; the details of the signalling pathways, as well as cell division, are neglected for simplicity.
The model comprises a pair of ordinary differential equations that govern the levels of Delta (d j ) and Notch (n j ) activity in each cell ( j). In dimensionless terms these are (Collier et al. 1996) where dots denote differentiation with respect to time and λ is the ratio of the decay rates of Delta and Notch activity. We remark that a single subscript is used to denote each cell for simplicity; however, generalisation to higher spatial dimensions is trivial.
In (1), f (d j ) and g(n j ) are feedback functions representing the coupling between adjacent cells and the inhibitory effect of Delta-Notch binding, respectively and d j denotes the mean level of Delta activity in the N cells surrounding cell j: The positive parameters k, h, a and b determine the feedback strength. Detail of the biological meaning of the exponents together with example molecules is given in Webb and Owen (2004). In vivo, Delta-Notch binding and the resulting Delta production are dependent on the cell's biochemical and biophysical environment. Experimental evidence suggests that environmental inhomogeneities are significant, which provides motivation for the consideration of microscale patterning in response to macroscale variation in tissue stimulation, leading us to include spatial variation in the parameters associated with the feedback functions f and g. For constant feedback strength Collier et al. (1996), Plahte (2001) and Plahte and Øyehaug (2007) remark that a striking feature of the Delta-Notch model is the robustness of the fine-grained pattern. 1 Webb and Owen (2004) have demonstrated that cell shape is an important factor in the patterns produced in cell-signalling models; furthermore, certain epithelial cells are not necessarily well-approximated by square cells (see Fig. 1). In view of these considerations, below, we extend our previous work (O'Dea and King 2011) to investigate the emergence of fine-grained patterning modes in hexagonal cells in response to variation in the parameters associated with the feedback strength.

Patterning and stability
By linear and bifurcation analysis we now determine the parameter regimes for which fine-grained spatial patterning will be produced by this model. In arrays of square cells, fine-grained patterning corresponds to patterns with a period of two cell lengths in each direction (and the periodic unit therefore contains four cells). In contrast, the shortest wavelength pattern which fits onto a discrete hexagonal mesh is of period three in each direction (the periodic unit containing three cells), corresponding to the dominant fine-grained patterning mode discussed in Sect. 2.1. In each case, in what follows, we will refer to these as patterns of period two and three, respectively. Due to the local nature of the juxtacrine signalling mechanism, considering patterns of Delta and Notch activity in such a periodic unit provides insights into more extensive arrays and the results obtained are crucial to the subsequent multiscale analysis. The system (1) reduces to a system of six coupled equations to which the solutions are either homogeneous or patterned with a wavelength of three cells. We remark that, in view of (2), in the period-three regime the connectivity of the periodic network is identical both in arrays of hexagonal cells and in one-dimensional strings. For the purposes of the following patterning and stability analysis, we therefore persist with the single index j to indicate each of the three cells in the periodic unit.
Via linear analysis of (1) it is straightforward to show that the homogeneous state is unstable and period three patterns exist if f (d * )g (n * ) < −2, where (n * , d * ) are the homogeneous steady states satisfying d * = g(n * ), n * = f (g(n * )) and denotes differentiation (Collier et al. 1996). Furthermore, the period three pattern is the fastest growing mode from the uniform state. Similarly, period three steady states n * j , d * j ( j = 1, 2, 3) become unstable to period three perturbations if f (d * j )g (n * j ) < −2, where d * j denotes the weighted sum of neighbouring steady-state values defined by (2). The transition from parameter regimes in which the homogeneous steady state is stable (so that adjacent cells exhibit identical levels of Notch activation and thereby attain identical cell fates) to that in which stable heterogeneous states emerge leads to the generation of fine-grained patterns of Delta and Notch upregulation (so that nearby cells diverge to different fates). We remark that since it is postulated that the level of Notch activity determines its fate, in the following, we will use the terms "up-" or "down-regulated" to refer to the level of Notch activation in each cell.
Bifurcation analysis indicates that the system admits three distinct pattern types of period three: (i) a ratio of two identically up-regulated cells to one down-regulated, (ii) a ratio of two identically down-regulated cells to one up-regulated and (iii) one up-regulated cell, one down-regulated and one intermediate; in the following, we shall refer to these patterns as type (i), (ii) and (iii) for clarity. Figure 2 shows typical numerical simulations displaying period three patterns of type (i) and (iii).
In the three cell system under consideration the various pattern permutations mean that types (i) and (ii) comprise three patterns each, whilst type (iii) contains six pattern configurations. Figure 3 shows the transition from one patterning mode to the next under variation of the exponents k and h which govern the feedback strength in the signalling model (similar bifurcation behaviour is observed under variation of k or h alone), indicating that there is a distinct range of parameter space in which each of the three types of pattern is stable, and that type (i) and (ii) patterns are stable for a significantly larger portion of parameter space than those of type (iii). We highlight that 'stability' as indicated in Fig. 3 refers to stability of period three patterns within the three-cell periodic unit. The various pattern configurations within each pattern type are highlighted on the relevant curves; Fig. 3b, c shows the configurations in detail. The solution branches indicate the level of Notch activation in each cell. For instance, in the type (i) patterning regime, referring to Fig. 3a, if the level of Notch activation in a cell follows the lower solution branch marked with circles, the remaining pair of cells follow the upper branch (marked with diamonds). Similarly, in the type (ii) regime, if the level of activation in a cell tracks the upper branch (filled circles), the remaining cells track the lower branch (marked with squares); in between these parameter ranges, type (iii) patterns are observed. Lastly, we remark that the bifurcation structure is significantly more complex than that associated with the generation of period two patterns (in each coordinate direction) in square cells (O'Dea and King 2011). It is important to note, however, that O'Dea and King (2011) considered the specific case of checkerboard patterns, for which the analysis corresponds exactly to that in a one-dimensional string of cells. In the case for which the periodic unit contains four distinct levels of upregulation (only two exist for checkerboard patterns) the behaviour is likely to be more complex. Figure 3a indicates that type (i) patterns branch from the homogeneous steady state subcritically; patterns containing three distinct levels of upregulation (type (iii) patterns) and type (ii) patterns are created at subsequent bifurcations labelled C a -C d . The bifurcation points C a , C c and C b , C d form solution pairs: the Notch activity in two of the three cells diverges from the supercritical bifurcations at C b and C c ; the remaining cell increases or decreases from its bifurcation point value C a or C d via the relevant solution branch. The pitchfork bifurcations at C b and C c imply that patterns with three distinct levels of Notch activation are generated from type (i) or type (ii) patterns by the local, symmetric divergence of two of the three cells. Figure 4a illustrates how this bifurcation structure changes under variation of the inhibitory feedback function Bifurcations at which patterns are created from the homogeneous state and between different pattern types are denoted C * and C a -C d , respectively. The homogeneous steady state is denoted n * ; for type (i) and (ii) patterns, up-and down-regulated cells are denoted '+' and '−'; type (iii) is denoted (1:1:1). Panes b and c show in more detail the behaviour near the branch points C b and C c at which the transition from types (i) and (ii) to type (iii) patterns occurs. The different type (iii) steady states are labelled 1, 2, 3 to denote increasing levels of Notch upregulation in each cell. Asterisks denote type (iii) solutions with associated type (i) or (ii) patterns parameter b, demonstrating that the five bifurcations shown in Fig. 3 collapse onto a single patterning bifurcation from the homogeneous steady state as b is reduced. Additionally shown are the resulting pattern-forming bifurcations obtained under variation of k, h for values of b at which the bifurcations collapse onto a subcritical bifurcation (Fig. 4b) and a supercritical pitchfork-like bifurcation (Fig. 4c). Here, only stable type (ii) and unstable type (i) patterns exist. The patterns shown in Fig. 2 are robust, being generated from random initial data as well as appropriate periodic or near-periodic initial conditions for a wide range of parameter values. We remark that longer-range patterns (of period greater than three) may be generated in this array of hexagonal cells by appropriate choice of periodic initial conditions and domain size. We note that these patterns are stable to periodic or aperiodic perturbations but that their domain of attraction is small and therefore are only observed for suitable initial data. Numerical simulations of (1) indicate that regions of these longer-range patterns form stable patterned distributions when competed with period-three patterns and, additionally, are able to invade the unstable homogeneous steady-state (invasive or competitive behaviour is generated via initial data comprising a region of stable longer-range patterning adjacent to the unstable homogeneous state or stable period-three pattern in the remainder of the domain). However, the introduction of random noise to periodic initial data, which is not itself a steady state of (1), seems always to result in the emergence of the dominant period-three patterning mode; similarly, an initial state comprising a region of unstable period-three pattern adjacent to a region of stable longer-range pattern, leads to the invasion of the unstable state with the corresponding stable period-three, rather than long-range, pattern (though the region of stable longer-range pattern remains). Figure 5a shows an illustrative stable pattern configuration obtained by competition of stable period-three and period-five patterns; Fig. 5b shows the modulated travelling wave of period-five pattern invading the unstable homogeneous state. Lastly, we note that the range of patterns formed is limited by the hexagonal mesh; in particular, patterns of even period in each coordinate direction are prohibited (though we remark that striped patterns of arbitrary period are, of course, permitted).
To summarise, the signalling model of Collier et al. (1996) displays robust finegrained patterning, consistent with those patterns observed in relevant biological signalling systems. In square cells, consideration of such patterns leads to the investigation of checkerboard-type patterns with period of two cell lengths in each coordinate direction. In hexagonal cells, the fine-grained patterning mode corresponds to patterns of period three in each direction. Collier et al. (1996) identified these patterns as the fastest-growing mode a via linear analysis. The results contained in this section show that the generation of fine-grained patterns in hexagonal cells is significantly more complex than their checkerboard counterparts in square cells. We have isolated the three possible types of fine-grained pattern, demonstrated their robustness via numerical simulation and shown how the bifurcation structure varies under variation of certain feedback parameters. Such information will prove crucial for the following multiscale Numerical simulations of (1) performed in a 60 × 60 hexagonal mesh illustrating the behaviour of longer-range patterns. a A stable configuration comprising a stable period three pattern of type (iii) (left) adjacent to a stable period five (right) Notch activation pattern. b The modulated wave of period-five Notch activation travelling through the hexagonal mesh ( j, l) and invading the unstable homogeneous steady state at successive times t = 0, 10, 25 along the line l = 30. Initial conditions comprise the domain split equally between the unstable homogeneous state ( j < 30) and stable period-five pattern ( j 30). All parameter values as in Fig. 2b analyses. In addition, via numerical simulations of (1) we have shown that patterns of wavelength greater than three may be produced by this model, though these patterns have small domains of attraction and are therefore only observed for suitable initial data, further exemplifying the dominance of the fine-grained model.
In the following sections, we employ a multiscale analysis to construct continuum models of the period-three patterning behaviour for parameter values near the patterning bifurcation points investigated above.

Model formulation
In this section, we employ a multiscale method to derive continuum models based upon the discrete system (1) that capture the period-three patterning phenomena described in Sect. 2.2. Figures 2 and 3 imply that introducing appropriate spatial variation of one (or more) of the model parameters will induce fine-grained patterning in certain regions of the domain, or transition from one patterning mode to the next; in the following, we construct continuum models which capture fine-grained pattern generation in response to macroscale tissue stimulation, corresponding to spatial variation in feedback strength. (Such formulations include spatially-uniform feedback as a special case, analysis of which will be especially instructive when considering pattern competition and travelling-wave behaviour.) In a biological context, such parameter variation might correspond to differences in the sensitivity of the cells to Delta-Notch binding, leading to the adoption of a particular programme of gene activation by a subset of the cell population according to spatial position. Such environmental inhomogeneities are a common feature of biological systems, such as in embryonic tissue growth. We note, however, that the parameter values employed in the remainder of this work are chosen to illustrate the interesting patterning behaviour exhibited by the model and investigated in Sect. 2.2 and are not motivated by specific biology.
To capture the different facets of the system (1) in the period-three regime within appropriate continuum representations, we consider in turn (a) the patterning bifurcation from the homogeneous state, (b) the transition between different period three patterning regimes and, (c) the behaviour for parameter values close to the set for which the separate bifurcation points collide leaving a single pitchfork bifurcation (see Figs. 3 and 4a, c). As remarked earlier, since in the parameter regime b > b * , illustrated in Fig. 3a, the bifurcation from the homogeneous state generating type (i) patterns is sub-critical, the asymptotic behaviour derived in case (a) will not be reflected in the observed solutions to (1). However, the simple model analysed here is ideally suited to illustrate our methodology; the analysis holds for alternative signalling models whose stability properties imply that the results are more immediately applicable. To construct a continuum model capturing period three patterning phenomena in response to spatial variations in feedback strength, a homogenisation process is required. To preserve the local periodicity, we assume that the variation of the spatially non-uniform parameter is slow compared to the variations of Delta and Notch activity: that is, we construct two-scale models which capture fine-grained periodthree pattering patterning phenomena. This separation of scales (between 'fast' and 'slow' variation) allows the integration of fine-grained (microscale) patterning within a macroscale framework.
The continuum model is derived as follows. Considering an array of hexagonal cells, we denote the distance between cell centres by δ 1, introduce slowly varying continuum variablesx = δ j,y = δl and represent the levels of Delta and Notch activity in the multiple-scales form: d j,l = d ( j, l,x,y, t), n j,l = n ( j, l,x,y, t), in which j, l represent the fast, andx,y the slow, spatial scales. The labelling scheme for hexagonal cells, together with the periodic repeating unit for the period three pattern is illustrated in Fig. 6. Additionally, since adjacent cells differ significantly in their Delta and Notch activity the patterning regimes, we write: and we will refer to the cells corresponding to Eqs. (3)-(5) as types a-c. We remark that since the discrete labelling scheme underpins our multiscale approach, we are unable, as it stands, to accommodate cell division, which would require a global relabelling; our analysis applies on the timescale of Delta-Notch mediated cell fate determination, which is significantly exceeded by that of the cell cycle (Hartenstein and Posakony 1900). Transforming to an orthogonal Cartesian coordinate system (x, y) = (x −y/2, √ 3y/2), expanding in Taylor series and exploiting the periodicity, the spatial coupling term d i for each cell type i may be written: Spatial variation in biochemical or biophysical conditions within the domain (leading to differences in feedback strength) is modelled by introducing slow spatial variation to the parameters k and h. For simplicity, we assume k = h and keep the remaining parameters fixed; we assume k = h = k(x, y). We expand around the bifurcation point under consideration (C * , C a -C d ; see Figs. 3 and 4) via: where i = a, b, c indicates expansions associated with each cell type, k = k * denotes the relevant bifurcation point, at which n i * , d i * are the steady states associated with each cell type, and ε 1 (δ), ε 2 (δ) 1. Additionally, we rescale time according to 1 will be chosen such that we analyse the equations at the timescale on which spatial coupling first appears.
We now pause to define some notation which will be of use in the following sections. We define linear operators L and M as follows: together with the averages d i p for each cell type i at the pth asymptotic order:

Emergence of patterns from the homogeneous steady state
In this section, we consider the generation of fine-grained patterns from the homogeneous steady state at the bifurcation point C * in the regime b > b * as illustrated in Fig. 3a. Employing the linear operators (10), the scalings ε 1 = ε 2 = ε 3 = δ 2 , the averages (11) and noting that at C * , (n i * , d i * ) = (n * , d * ) for i = a, b, c, the equations governing Delta activity in each cell type i at each order are: and Notch is governed by: The functions F p ,F p and G p ,Ĝ p are the O(δ 2 p ) perturbations to g, f and g , f , expanded around the bifurcation point value k * , and defined by: these being evaluated at k = h = k * , d = d * and n = n * , and in view of the expansion (7) depend on the spatial coordinates x and y.
The O(δ 2 ) linear system is of rank four and therefore has two levels of degeneracy, reflecting the additional degree of freedom implied by the analysis of a periodic unit of three cells (in comparison to that observed when analysing checkerboard patterns) in which the different cell types may be freely interchanged. Noting that f (d * )g (n * ) = −2 (see Sect. 2.2) and n i * = n * for i = a, b, c, the combinations: where i = b, c and p = 1, 2, . . . , are identically zero. These correspond to the two eigenvectors (with zero eigenvalue) of the linear system, based on which we make the ansatz: where A and B are to be determined and, from (13), (16) we obtain: C = g (n * )F 1 (d * ) + F 1 (n * ) /3. We remark that through appropriate choice of A, B, all the types of possible period-three pattern discussed in Sect. 2.2 may be represented by this ansatz, including the additional pattern configurations obtained under cyclic permutation of the periodic unit. Employing the linear combinations (19) to eliminate the O(δ 2 ) perturbations n i 2 , d i 2 , equations (12)-(17) may be expressed as the following pair of coupled partial differential equations for A(x, τ ), B(x, τ ): where , γ , α 1 and β 1 are defined by: We reiterate that, in view of Eqs. (7) and (18), γ = γ (x, y) and β 1 = β 1 (x, y), while and α 1 are constants. The corresponding Notch activity is given by Eq. (13). As in the case of square cells (see O'Dea and King 2011), Eqs. (21) and (22) demonstrate that, on the macroscale, the juxtacrine signalling interaction manifests itself as linear diffusion, with effective diffusivity 1/ , despite the short-range patterning under consideration.
This continuum model differs from that derived in O'Dea and King (2011) (in which a checkerboard patterning regime in square cells was considered) by the quadratic (rather than cubic) form of the nonlinearity, due to the asymptotic scalings appropriate to reflect the solution behaviour near C * . Additionally, we have expressed our equations in a more general form to enable the wider array of period three modes to be captured.
We note that the form of the continuum model is not crucially dependent on the details of the underlying signalling system; additional interaction terms (such as those of the form n j d j or n j d j as employed by Owen and Sherratt (1998) in a more complex juxtacrine signalling model, accommodating ligand, free receptors and bound receptors) do not materially affect the diffusive interaction or the quadratic nonlinearity. As noted above, since the bifurcation at C * is subcritical, the behaviour captured by Eqs. (21) and (22) will not be reflected in the observed solution behaviour; however, the analysis presented here will apply to similar systems whose stability does not conform to that shown in Fig. 3. For constant parameter values, the uniform steady states of (21) and (22) are determined by solving a pair of coupled quadratic equations. (We highlight that these are "uniform" in the sense that they are constant for their associated cell type, adjacent cells attaining different steady states.) The signs and magnitudes of the parameters γ, α 1 and β 1 in (21), (22) determine the existence of such solutions, reflecting the stability properties of the underlying discrete system. In the case investigated here, multiple steady-state solutions to (21), (22) exist on both sides of the bifurcation point, reflecting its inability to differentiate between stable and unstable states; in contrast, a supercritical transition, such as the pitchfork bifurcations at C b , C c , would imply a reduction in the number of real steady states as the bifurcation point is traversed.

Transitions between patterned states
In this section, we derive a continuum model to capture the transition between period three patterns of types (i) and (iii) or types (ii) and (iii) bifurcation points labelled C a -C d , as illustrated in Fig. 3. We note that, in contrast to the subcritical behaviour analysed in Sect. 3.1.1, in this case, the transitions occur at supercritical pitchfork-type bifurcations; the continuum model we shall derive is therefore immediately applicable to the dynamics of the underlying discrete system.
Adopting the notation employed in Sect. 3.1.1, the appropriate linear combinations in this case are: wherein i = b, c and n i * , d i * now denote the levels of Notch and Delta activity in each cell type at the bifurcation point between pattern types, again denoted k * in each case.
Considering the bifurcation points C b , C d at which the pattern has a ratio of two upregulated cells to one downregulated cell (see Fig. 3), it may be observed that the downregulated state is an order of magnitude smaller than the upregulated state and we therefore scale the downregulated state with δ. Figure 3 indicates that at the bifurcations C a , C c , the steady states of Notch activity are of comparable order; however, the corresponding up-and downregulated steady states of Delta activity display an order of magnitude disparity (not shown) and a corresponding rescaling must therefore be chosen when considering these bifurcations. We remark that this change in scalings is also revealed from the asymptotic analysis.
In the more general case for which the downregulated steady-state value of Notch activation is of the same order as the level of upregulation, the scaling employed in Sect. 3.1.1 is appropriate. Considering a bifurcation of the form C b , C d (with steady states of comparable order), the calculation follows the same method as that outlined in Sect. 3.1.1 and we obtain: in which , α 1 , γ and β 1 are defined by (23) (26) and (21) having identical form. The parameters α-χ are known functions of both up-and down-regulated steady states: n i * , d i * , d i * (i = a, c), but are too cumbersome to include here; these are included in Appendix A for completeness. We now consider in detail the behaviour at C b , C d under the configuration (n a * , n a * , n c * ), where n c * = δn c * and the steady-states n a * andn c * are O(1). In this case, at first order, the level of Notch activation in cells of type a and b diverges at the pitchfork bifurcation C b , the activity in type c cells remaining constant. Such behaviour corresponds to choosing A = 0 in Eq. (20). Here, the appropriate scaling is ε 1 = ε 3 = δ 2 and ε 2 = δ so that spatial and temporal variation enters at O(δ 3 ). The remaining two configurations obtained under cyclic permutation may be analysed in a similar way and correspond to B = 0 and A = −B in Eq. (20) in each case. We note that an equivalent derivation provides an equation governing the transition between pattern types (ii) and (iii) at C a , C c , but is omitted.
The equations at O(1), O(δ) and O(δ 2 ) are of similar form to (12)-(17); the equations governing Notch and Delta activity in each cell type are as follows: wherein i = a, b, c and F 1 ,F 1 , G 1 andĜ 1 are defined by (18) where is defined by Eq. (23) and α 2 and β 2 (x, y) are defined: Equation (34) is invariant under the transformation B → −B, reflecting the symmetric divergence of activity in cells of type a and b at the supercritical pitchfork bifurcation C b ; the level of activity in cells of type c remains constant at O(δ). The continuum model (36) is therefore identical in form to the reaction-diffusion equation derived in O'Dea and King (2011) when considering a period-two patterning regime in square cells. We reiterate that the patterning mode under consideration and the associated scalings, rather than the details of the underlying signalling model, result in the form of the partial differential equation (34).
For constant β 2 (i.e.k = constant), the uniform steady states of (34) are B = (0, ± √ −β 2 /α 2 ); 'uniform' here means that they are spatially-constant in cells of type a. For k < k * we find β 2 > 0 and only the trivial solution B = 0 exists, corresponding to the maintenance of the type (i) heterogeneous steady state: the periodic unit remains in the state n ∼ (n a * , n a * , δn c * ). As the feedback strength increases, β 2 changes sign: k k * , β 2 0; in this regime non-trivial roots exist and the behaviour in the periodic unit is n ∼ (n a * , n a * , 0) + δ(± √ −β 2 /α 2 , ∓ √ −β 2 /α 2 ,n c * ) corresponding to the emergence of type (iii) patterns. Considering only the local value of the feedback strength β 2 = β 2 (x, y), we therefore expect the generation of patterned states to reflect the pitchfork bifurcation at k * . For brevity, hereafter we will refer to this patterning behaviour as that influenced by "local feedback strength".
We note that a similar analysis may be performed at the bifurcation points C a , C c yielding a corresponding result.

Two parameter expansion
The above multiscale analyses have produced continuum models, based upon an underlying discrete patterning system, for the emergence of (unstable) patterns from the homogeneous steady state (Sect. 3.1.1) and the (stable) transitions between patterns of types (i) and (ii) and types (ii) and (iii) (Sect. 3.1.2) in parameter regimes for which these patterning bifurcations are widely spaced. Inspection of the model equations derived in Sects. 3.1.1 and 3.1.2 reveals that the generation of, and transition between, period-three patterns is dependent on α 1 (for clarity, we emphasise that α 1 = α 1 (n a * , d a * , k * , b), n a * and d a * being the steady states in cells of type a at the bifurcation point k * for a specific value of b). The multiscale analyses presented in Sects. 3.1.1 and 3.1.2 correspond to the bifurcation structure illustrated in Fig. 3 for which α 1 O(1). In the case α 1 = 0 (b ≈ 1.927), the five bifurcations collapse onto one pitchfork-like patterning bifurcation generating stable type (ii) and unstable type (i) patterns from the homogeneous steady state (see Fig. 4c). Here, we analyse the behaviour near the bifurcation at α 1 = 0; inspired by Fig. 4 we consider in addition to the expansions (7)-(9), in which (d i * , n i * ) = (d * , n * ) denote the homogeneous steady-state values of Delta and Notch activity at b * , k * , the bifurcation point values at which α 1 = 0. In this regime, the model (21) and (22) breaks down, providing only the trivial steady state; instead, a continuum model for period three patterns being generated from the homogeneous steady state is obtained by choosing ε 1 = ε 3 = ε 4 = δ 2 and ε 2 = δ and spatial and temporal variation enters at O(δ 3 ) as in Sect. 3.1.2. The equations at each order are identical to (28)-(33) with f and its derivatives evaluated at d i * = d * ; here, F 1 ,F 1 , G 1 ,Ĝ 1 are the O(δ 2 ) perturbations to f, g and f , g expanded around b * , k * . As in the previous analysis of pattern emergence from the homogeneous steady state, the linear combinations (19) (32) (20) with C = 0, (28)-(33) may then be expressed as the following pair of coupled partial differential equations for A(x, τ ), B(x, τ ): where α 2 and β 2 are defined by (36) and evaluated at α 1 = 0, d a * = d * . The corresponding Notch activity is given by Eq. (13). The analyses presented in Sects. 3.1.1 and 3.1.2 describe the behaviour near the individual bifurcation points C * or C a -C d ; above, via the isolation of the combination α 1 which dictates the bifurcation structure of the system (except in the special case for which the downregulation of Notch activity in one cell is extreme), we have derived a more general formulation which captures the generation of the five bifurcation points illustrated in Fig. 3. Furthermore, the model corresponding to the structure illustrated in Fig. 4c may be recovered by settingb = 0.
We remark that the nonlinearities in the model (38), (39) are of a different form to those in the models (21), (22) and (26), (27), obtained by analysing pattern generation in the regime α 1 O(1) for steady states of comparable order. However, by considering the matching between the inner (in parameter space) formulation (38) and (39) and the outer formulations (21), (22) and (26), (27), 2 it can be shown that these models are consistent in the limitsb → ∞, b → b * . This calculation is summarised in Appendix B for completeness.

Summary
We have employed a multiscale analysis to derive continuum models based upon the discrete Delta-Notch signalling model that describe both the emergence of period three patterns from the homogeneous steady state and the transition between two types of period-three pattern in appropriate parameter regimes.
In the former case, the bifurcation from the homogeneous state generating type (i) patterns is sub-critical, and so the asymptotic behaviour derived will not be reflected in the observed solutions to (1). However, the analysis contained herein serves to illustrate our methodology; similar continuum models may be derived for alternative signalling systems whose stability properties imply that the results are more immediately applicable. On considering the transition between patterned states, two distinct regimes were identified. In the case for which the heterogeneous steady states of Delta and Notch activation are of comparable order, the generation of patterns from the homogeneous state and the transition between patterning modes are governed by models of a similar form. In each case we obtain a pair of coupled reaction-diffusion equations with quadratic nonlinearity: Eqs.  (2011), in which a cubic nonlinearity is obtained. This disparity is due to the asymptotic scalings required to capture solution behaviour near the bifurcation points in the multiscale analysis. An exception is the transition between patterning modes in parameter regimes for which the downregulation in one of the three cells is extreme (the transition at the pitchfork bifurcations C a or C d shown in Fig. 3). Here, since the level of activation in the downregulated cell remains constant (to first order), the resulting continuum model (34) is identical in form to that derived in O' Dea and King (2011) in which the symmetric divergence of a periodic unit of two cells at a pitchfork bifurcation was considered.
The model Eqs. (21), (22) and (26), (27) are expressed in a more general form than that in O' Dea and King (2011), leading to a system of PDEs in place of the single PDE obtained in O' Dea and King (2011). This reflects the more complex pattern forming behaviour exhibited in a hexagonal mesh. In the case of checkerboard patterns, the pattern-forming bifurcation is of pitchfork type and the possible pattern configurations are equivalent up to a phase shift of a single cell. The asymptotic behaviour is therefore captured by one PDE (with cubic nonlinearity) and the Notch activation in alternate cells (n 1 , n 2 , say) diverges symmetrically from the bifurcation point; the alternative pattern therefore displays the corresponding, but inverse, behaviour, viz: ±n 1 = ∓n 2 . In the case of patterns of period three, the various pattern configurations associated with each mode mean that such a symmetry is not sufficient to capture the entire behaviour. A second PDE is required to provide an additional degree of freedom and these patterns may then be captured via suitable choice of the slowly-varying functions A and B. In (34), however, due to the specifics of the problem at hand, we are led to consider the case for which A = C = 0 and the resulting single equation models the transition from type (i) patterns in the configuration (n a * , n a * , δn c * ) to corresponding patterns of type (iii) only. The return to a single PDE and cubic nonlinearity obtained in this case reflects that the third cell in the periodic unit remains constant to first order, the remaining two cells diverging symmetrically, as is the case for checkerboard patterns. Figure 4 reveals that the bifurcation structure is controlled by the feedback parameter b, reflected in the appearance of the combination α 1 in the asymptotic analysis; the five bifurcations illustrated in Fig. 3 collapse onto a single supercritical pitchfork-like bifurcation at a critical value b * (corresponding to α 1 = 0). Inspired by this, we performed a two-parameter expansion, constructing a continuum model for the unfolding bifurcation at k * , b * . In this case, the second order equations become degenerate and the resulting reaction-diffusion equation has a cubic nonlinearity as in O' Dea and King (2011), reflecting that the behaviour in two of the three cells is identical, tracking the upper solution branch of the pitchfork. The remaining cell follows the lower branch; the alternative patterns are obtained by inverting this behaviour. However, a system of PDEs is nevertheless obtained, reflecting the additional possible pattern configurations. By considering matching between this formulation and the preceding analyses in appropriate limits, we confirmed that the various continuum limits are consistent and that the model derived for the behaviour near b * contains the other models as limit cases.

Steady states
As noted above, the continuum model (21), (22) is unable to predict the observed behaviour of the discrete model (1) since the period three patterns branch subcritically from the homogeneous steady state at C . In place of time-dependent numerical simulations of (1), we therefore compare the unstable steady states of (1) near C with those of (21) and (22) to demonstrate that our multiscale analysis does in fact reflect  Figure 7a shows how the steady state levels of Notch activity compare to the bifurcation behaviour of the nonlinear model (1) as the parameter k is varied; Fig. 7b shows this in more detail. Clearly, the steady states of the continuum model form a good approximation to the (unstable) underlying discrete behaviour close to the bifurcation point; further from the bifurcation point, the transcritical-type behaviour predicted by the continuum limit diverges from the nonlinear model.

Simulations
To illustrate the dynamic behaviour of our continuum formulations, we now present numerical solutions of (34) on the domain 0 x L , 0 y L and compare these to corresponding numerical simulations of the discrete system (1). These simulations capture the transition between type (i) and type (iii) patterns. We choosê k(x, y) = 2x/L − 1 to demonstrate how spatial variation in feedback strength dictates the emergence of fine-grained patterning in the domain. As noted in Sect. 3.1.2, the uniform steady steady states of (34) are (0, ± √ −β 2 /α 2 ). Considering only the local value of the feedback strength β 2 = β 2 (x, y), we therefore expect the generation of patterned states to reflect the pitchfork bifurcation at k * . In addition to the comparison with the underlying discrete model, in the following, we demonstrate how solutions to (34) with spatially-varying feedback differ from the pitchfork bifurcation behaviour that arises from considering only this 'local feedback' strength. Simulations of (38) and (39) capture the emergence of stable type (ii) patterns from the homogeneous steady state; while the specifics of the pattern transition represented by (38) and (39) are significantly different to those described by (34), the behaviour of the relevant variables at first order is qualitatively similar (two of the cells following the upper solution branch of the pitchfork-like bifurcation, the third cell tracking the lower) and simulations are therefore omitted for brevity. Equation (34) was solved by discretising on a uniform spatial grid and employing the initial value problem solver ode15s in MATLAB; in the following simulations we choose the grid spacing x, y as x = y = 10 −4 . We note, however, that the grid spacing is not prohibited from exceeding significantly the diameter of a single cell since the cell-scale variation has been systematically averaged out. Since we expect macroscale variation over the domain due to the chosen form ofk(x, y), we impose no-flux conditions at x = 0, L in place of a macroscale periodic domain imposed in the discrete simulations shown in Fig. 2.
In Fig. 8a, the steady-state pattern in a 19×19 section of a small 2D array of hexagonal cells is shown, obtained from simulation of the discrete system (1) indicating the transition from one patterning mode to the next. So that patterning effects are evident on this small domain, we choose δ = 0.5, L = 10.5. In Fig. 8b the steady-state value of the level of Notch activation in cells of type a and b as predicted by (34) are compared to a corresponding simulations of (1) in an array of 900 × 900 hexagonal cells (L = 9, δ = 0.01), demonstrating that a smooth transition to from type (i) to type (iii) patterning is obtained. In contrast, considering only the local feedback strength (see Sect. 3.1.2) predicts a pitchfork bifurcation from spatial homogeneity to patterning; far from the bifurcation k * , the local feedback approximation shows good quantitative agreement with the numerical simulations.

Travelling waves
In this discrete model, spatial patterns in Delta and Notch activation can be generated by a modulated travelling wave, behind which a regular pattern forms. Comparison between the wave speed of pattern invasion predicted by our continuum models and that obtained from simulation of the discrete system provides a quantitative description of the accuracy of the continuum approximation and an indication of its range of applicability in parameter space. Furthermore, such phenomena are also observed in vivo; for example, the morphogenetic furrow in the retina of Drosophila, behind which a regular pattern of photoreceptors surrounded by support and pigment cells is created (Owen 2002;Plahte and Øyehaug 2007).
Travelling-wave solutions to (34) may be constructed by considering a type (iii) pattern-forming wave invading the unstable type (i) state on an infinite domain. As noted in O' Dea and King (2011), while linear propagation front travelling-wave solutions of the form B(x + vt) may be constructed, these correspond to 'pushed' wave fronts. Equations (1) and (34) describe 'pulled' wave fronts (Hadeler and Rothe 1975;Stokes 1976;Rothe 1981;Plahte 2001), whose minimum wavespeed may be calculated by analysing the behaviour at the leading edge of the wave via: The wavenumber in each coordinate direction k = (k x , k y ) and the minimum speed v are determined from the dispersion relation obtained from (34) together with the corresponding repeated root condition. Considering waves travelling in the direction y = 0, we obtain v = 2 √ β 2 / .
Pattern-forming travelling waves may be induced in simulations of the discrete and continuous models through choice of suitable initial data; the following simulation corresponds to an initial state consisting of stable type (iii) pattern in the region 0 < y < L , 0 < x < L/30 and unstable type (i) in the remainder. We choose L = 100, δ = 0.05 and, as in Sect. 4.2, no-flux conditions are applied at x = 0, L in place of periodicity. Figure 9a shows the position of the midpoint (defined as the halfmaximal value of Delta activation in cells of type a) of a modulated travelling wave of a type (iii) pattern invading a region of unstable type (i) pattern obtained from numerical simulation of (1) and (34), indicating that the numerical solution of the homogenised model faithfully reproduces the behaviour of the underlying discrete system. While providing a good qualitative approximation, the 'linearly-selected' wavespeed deviates from this observed behaviour; analysing travelling waves closer to the bifurcation point reduces this deviation: Fig. 9b indicates how the speed of the invading wave predicted by the discrete and continuous systems varies as a function of δ. In these simulations, the wavespeed was calculated in the central part of the domain to minimise boundary effects. Clearly, the minimum wavespeed predicted by the continuum model provides an excellent approximation to the underlying discrete behaviour, even for parameter values relatively far from the bifurcation point; this approximation improves as the bifurcation point is approached.

Pattern competition
In each of the three patterning regimes depicted in Fig. 3, more than one period-three heterogeneous state is stable: there are three possible type (i) and type (ii) patterns, and six type (iii) pattern configurations. Since, in each case, the pattern configurations are equivalent (up to a suitable permutation of the periodic unit) and equally  (1) and (34). The domain is of side L = 100, (δ = 0.05) and all parameters are as in Fig. 8 except, in b, δ as stated stable, a priori, neither has the propensity to invade the other: In Fig. 10, we demonstrate how the system (1) evolves from an initial state containing more than one type of stable pattern configuration. Figure 10a illustrates that two type (i) patterns may coexist stably, with a smooth transition region in between. This behaviour is observed irrespective of the proportion of the domain occupied. In contrast, Fig. 10b indicates the evolution of the system from an initial state comprising one stable type (iii) pattern enclosing another, showing that these patterns may invade leading to a state where only the pattern which initially occupies the majority of the domain exists. In this regime, the multiscale analysis outlined in Sect. 3.1.2 applies and comparison with the homogenised model may be performed; as demonstrated in O'Dea and King (2011) the initially (approximately) square enclosed region becomes circular as it is invaded as expected from the corresponding continuum system (figure omitted for concision). Referring to Fig. 3b, the patterns shown in 10(b) are type (iii) patterns in (3,2,1) and (2,3,1) configuration which lie on the same solution curve emanating from the associated type (i) pattern (+, +, −) configuration. The remaining two matching configurations (2,1,3) and (3,1,2) which branch from the bifurcation point C b also display similar invasive behaviour; all other competing pattern configurations coexist stably in the domain with smooth transitions between patterning states (cf. Fig. 10a; omitted for concision).
The non-invasive behaviour displayed by the discrete model is termed "pinning" and illustrates an interesting feature of such patterning systems in which multiple modes are stable: a minority of patterning modes display invasive behaviour, the remainder are able to pin, forming discrete regions of stable pattern. Such invasive behaviour is exhibited only by 'compatible' patterning modes, such as those highlighted in Fig. 3b, c. These correspond to the symmetrical pairs of pattern configurations, reflecting the pitchfork structure of the pattern-generating bifurcation; the remaining configurations are not invasive. In view of this, the longer the characteristic wavelength of the pattern under consideration, the more pathological this pinning behaviour will 1) along the line y = 1 from an initial state comprising a (3, 2, 1) configuration type (iii) pattern enclosing a (2, 3, 1) configuration (k = k * + δ 2 ). Here, the continuum model solutions are indicated by circles. In b the Delta evolution in the remaining cell is omitted for clarity (the extreme Delta upregulation in type c cells renders the graphs unclear). Except as stated, all model parameters are as in Fig. 2 be, since permutation of the periodic unit implies a larger number of possible configurations. Pinning of patterns can therefore be expected to occur even for parameter values very near the bifurcation point, for which behaviour more appropriate to a continuum description would be anticipated. Far from the bifurcation point, strong discrete effects may prevent pattern invasion, and pinned pattern fronts are generally observed even for 'compatible' patterns. Pinning is exhibited in the related juxtacrine signalling model of Owen and Sherratt (1998) and will be discussed in a subsequent publication. The corresponding continuum limit is unable to capture such behaviour, its description being restricted to the relatively few compatible pattern configurations emanating from the bifurcation point; continuum models of patterning systems with a wide range of patterning modes are therefore severely restricted in the range of patterning behaviour they are able to capture.

Discussion
In this paper, we have examined thoroughly the range of fine-grained (period-three) patterns formed in an array of hexagonal cells which exhibit intercellular signalling behaviour governed by the generic Delta-Notch signalling model of Collier et al. (1996). This model takes the form of a pair of ODEs, representing the activation of juxtacrine (membrane-bound) signalling molecules and their receptors in each cell; a key assumption is that receipt of inhibition reduces the ability of a cell to inhibit its neighbours. For sufficiently strong feedback, this negative feedback loop is able to amplify differences in Notch activation between cells, the level of which determines the cell's fate (low levels lead to the progression of cells to the primary fate, whilst high levels relegate the cell to the secondary fate) and is therefore an important cell fate regulator in densely-packed tissues.
We employ the method outlined in O' Dea and King (2011) to derive continuum models based upon the underlying discrete model which are capable of describing the fine-grained patterning which arises under particular assumptions on the model parameters associated with signalling feedback strength. It has previously been thought that the short spatial scale of oscillation of these patterns necessitates analysis via discrete methods; however, by considering period-two (checkerboard) patterning in square cells, in O' Dea and King (2011) we demonstrated that such microscale variation may be accommodated within models derived via the straightforward assumption of slow variation and judicious choice of variables, exploiting the local periodicity displayed in specific parameter regimes. The resulting continuum model demonstrates that this interaction may be represented as a diffusive process even in a fine-grained patterning regime.
Inspired by the fine-grained patterns observed in early development of some tissues, and the discrete analysis of Webb and Owen (2004) which demonstrated that cell shape has a strong influence on the patterns produced in cell signalling models and the morphology of certain epithelial cell types, we extend our previous work to analyse fine-grained patterns in an array of hexagonal cells. In contrast to the checkerboard patterns analysed in O'Dea and King (2011) (which, although the periodic unit contains four cells, are equivalent to period-two patterns in a one-dimensional string), here, fine-grained patterning corresponds to a wavelength of three cell lengths in each direction (and the periodic unit therefore contains three cells). While still fine-grained in nature, the increased pattern wavelength leads to a far richer pattern-generating bifurcation structure than that investigated in O'Dea and King (2011) and we are led to consider separate continuum models for the generation of patterns from the homogeneous state and the transition between patterning modes. Via a two-parameter expansion, we additionally derive a more general formulation capturing the emergence of the bifurcation structure.
The continuum model derived in O'Dea and King (2011) took the form of a single reaction-diffusion equation with cubic nonlinearity (algebraic formulae providing the three remaining solutions), the juxtacrine interaction manifesting itself as linear diffusion. In this case, due to the symmetry of the two possible pattern configurations, the entire checkerboard patterning behaviour of the system may be represented by one such equation; here, the form of the continuum model is determined by the patterning regime (and associated scalings) that we wish to capture. The crucial difference between the continuum formulations obtained in square and hexagonal geometries is that, since the various period-three pattern configurations in hexagons conform to a more complex symmetry than checkerboard patterns, an additional degree of freedom is necessary, resulting in a pair of coupled reaction-diffusion equations.
In the case for which the bifurcation points are widely spaced (α 1 O(1)), the nonlinearities in the systems of PDEs representing both generation of patterns from the homogeneous steady state and transition between patterning modes are quadratic rather than cubic. A notable exception is Eq. (34) which applies to the specific case where downregulation in one of the cells is extreme; here, we obtain a model of identical form to that derived in O' Dea and King (2011), reflecting the fact that the third cell in the periodic unit remains constant to first order, the remaining two cells diverging symmetrically as is the case for checkerboard patterns. A cubic nonlinearity is also obtained in the parameter regime for which 0 α 1 1 in which the five bifurcation points shown in Fig. 3 collide to leave a single pitchfork bifurcation generating type (i) and (iii) patterns from the homogeneous state; however, in this case a system of PDEs is nevertheless required to capture the various possible pattern configurations.
We remark that the continuum formulations are general: the form of the nonlinearity is influenced by the required asymptotic scalings rather than the specific details of the signalling model under consideration. Ongoing analysis of period two patterning phenomena in a more complex juxtacrine signalling framework, modelling explicitly the number of ligand molecules and free and bound receptors on each cell results in a continuum model of the form presented in O' Dea and King (2011). It is important to note that due to the bifurcation structure of the specific signalling model under consideration (see Fig. 3), the continuum representation (21) and (22) of pattern generation from the homogeneous steady state in the regime α 1 O(1) is not reflected in the observed solution behaviour; however, such a model will apply in general to similar signalling systems whose stability properties imply that this analysis is more immediately applicable.
The ability of our formulation to capture the underlying discrete dynamics is illustrated via comparison between the nonlinear steady states of the continuum and discrete models (in the case of Eqs. (21) and (22) for the reasons outlined above) and between time-dependent numerical solution of the continuum model and simulations of the corresponding discrete system in a variety of scenarios. For concision we present time-dependent simulations of the transition model (34) only; corresponding behaviour is observed in solutions to the remaining continuum equations (26), (27), (38) and (39) in appropriate parameter regimes.
We demonstrate firstly that the steady-states of (21) and (22) provide a good approximation to the bifurcation behaviour of the nonlinear model (1) for parameter values close to the bifurcation point. Further from the bifurcation point, the nonlinear behaviour rapidly diverges from the transcritical-type behaviour predicted by the continuum limit. Secondly, we consider the time-dependent transition between patterns of types (i) and (iii). We demonstrate that under slow variation of parameters associated with feedback strength, the continuum model faithfully reproduces the transition between patterning types demonstrated by the underlying discrete system. We also compare the speed of propagation of travelling waves of stable type (iii) patterns invading a region of unstable type (i) patterns, demonstrating excellent qualitative agreement between numerical solution of, and the speed of a linearly-selected travelling-wave solution to, our homogenised model and the behaviour of the underlying discrete system, even for parameter values far from the bifurcation point.
Lastly, we investigate the competition between the multiple patterning modes. We show that pairs of 'compatible' type (iii) patterns exhibit invasive behaviour, the system evolving to contain the pattern configuration initially occupying the larger part of the domain. These compatible modes are those pattern configurations which reflect the symmetry of the underlying bifurcation structure. Competition between the remaining patterning modes results in pinned behaviour in which the system supports multiple stable patterning regions with a smooth transition region in between. Such behaviour is not reflected in the continuum description of the patterning system, which highlights a limitation of the ability of continuum models of the type presented herein to capture patterning behaviour in systems with a wide range of stable modes. We remark that such shortcomings are likely to be more pathological in longer-range patterning systems in which a wider range of pattern configurations exist for each patterning mode. In addition, the dynamics of such discrete systems may be strongly influenced by the underlying lattice geometry (Cahn et al. 1998), a feature not captured in the continuum representations presented here. A subsequent publication will investigate lattice-induced anisotropy in discrete diffusion and cell signalling systems. We remark that the stable "pinned" patterning behaviour highlighted here may be interpreted biologically in terms of a stem cell niche, as it corresponds to the constraining of cells of a distinct phenotype to a distinct, fixed region. As such, further investigations into pinning phenomena in discrete systems may have implications cell differentiation research.
The aim of this work was to develop continuum representations of the discrete patterning phenomena arising from models of intercellular signalling. Such an approach enables the integration of important microscale behaviour into tissue-scale models. The specific intercellular signalling model analysed in this work was chosen in large part for illustrative purposes, its relative simplicity admitting analysis; however, the mathematical techniques employed here may be applied to more elaborate spatially-discrete signalling models with little difficulty, enabling integration of discrete cell-signalling mechanisms into tissue-scale models in a wide variety of scenarios. Our methodology is, however, limited in two important ways. Firstly, it is unable to account for large microscale stochastic noise in cell signalling activity, due to its reliance on a separation of scales between cell signalling variation and feedback strength. Our analysis therefore applies only in those cases for which the noise is sufficiently small to be averaged out on the macroscale, or has no affect on the overall stability of patterns (see Rudge and Burrage 2008). Secondly, although our methodology can be extended to reflect cell growth or domain deformation (see, e.g., Fozard et al. 2009), we are unable, as it stands, to accommodate cell division, which may require global relabelling of cells. Both of these phenomena may be modelled with ease within an individual-based model and their inclusion in continuum models poses interesting challenges for future work. Other biologically-relevant extensions include investigations of continuum models capable of capturing multi-mode patterning states. Additionally, the patterns typically observed in biological systems have approximate short-range periodicity with defects; capturing such phenomena will provide many interesting mathematical challenges and is integral to the development of relevant continuum models of cell signalling behaviour.
homogeneous case,d * =d * andD † = D † = 0. When considering an outer expansion around the heterogeneous state C a in the configurationd * = (d a * ,d a * ,d c * ), we setD † = (0, 0,D † ) and similarly for D † . Corresponding forms are chosen for the remaining patterning configurations at C a or expansion around the heterogeneous state at C c . Matching to the homogeneous outer solution, in the limitsb → ∞ and b → b * and examining the equations at each order, the nonlinearities in (21), (22) are replaced by terms of the following form: The coefficients α † and β † are known functions of the expansion of α 1 and β 1 in the limit b → b * , D † ,d * and d * ; similarly, D † is a known function of the model parameters. These are omitted for brevity. The corresponding nonlinearities in (38), (39) have an identical form.
In the case of the heterogeneous steady state the nonlinear term in (26) is identical to (46); the nonlinearity in (27) has the form: in which a 1 -f 1 are functions of the expansions of α 1 , β 1 ,d a * ,d c * and d * independent ofD † = a 2 d c † 1 − b 2 d a † 1 + g 2 (in which a 2 , b 2 and g 2 are similarly dependent on the steady states and model parameters) which are too long to include here. Corresponding nonlinear terms are obtained from (38) and (39) in the limitb → ∞.