Abstract
We hypothesize that a population of migrating cells can form patterns when changes in local strains owing to relative cell motions induce changes in cell motility. That the mechanism originates in competing rates of motion distinguishes it from mechanisms involving strain energy gradients, e.g. those generated by surface energy effects or eigenstrains among cells, and diffusion–reaction mechanisms involving chemical signalling factors. The theory is tested by its ability to reproduce the morphological characteristics of enamel in the mouse incisor. Dental enamel is formed during amelogenesis by a population of ameloblasts that move about laterally within an expanding curved sheet, subject to continuously evolving spatial and temporal gradients in strain. Discretecell simulations of this process compute the changing strain environment of all cells and predict cell trajectories by invoking simple rules for the motion of an individual cell in response to its strain environment. The rules balance a tendency for cells to enhance relative sliding motion against a tendency to maintain uniform cell–cell separation. The simulations account for observed waviness in the enamel microstructure, the speed and shape of the ‘commencement front’ that separates domains of migrating secretorystage ameloblasts from those that are not yet migrating, the initiation and sustainment of layered, fractureresistant decussation patterns (crossplied microstructure) and the transition from decussating inner enamel to nondecussating outer enamel. All these characteristics can be correctly predicted with the use of a single scalar adjustable parameter.
1. Introduction
Turing demonstrated that diverse spatial patterns can form in the concentrations of chemical factors when two or more different molecules diffuse through a body at different rates, while their concentrations are simultaneously changing owing to chemical reactions [1]. The reaction–diffusion phenomenon has been a cornerstone of theories of pattern formation in nature by chemical signalling [2].
Others have addressed an alternative possibility that gradients in strain energy (either stress in a volume or surface energy) generated by cells as they migrate or proliferate might be the symmetrybreaking mechanism underlying pattern formation. Many experiments have shown that strain is a factor in cell migration, proliferation and differentiation [3–11]. When incorporated into physical theories, such effects can account for various shapes and patterns arising in organogenesis [12]. In the ‘differential adhesion hypothesis’, cells of different type are sorted spatially if their celladhesion energies differ [13–15]. Static analysis of the related phenomenon of cell surface tension can account for some patterns among cells of single type [16]. Network structures similar to trabecular bone can be generated if cell growth is correlated with stress variations around local fluctuations in material density [17]. As cell alignment is influenced by the stiffness of material to which the cell adheres, alignment patterns can be induced by a feedback mechanism between cells and secreted extracellular matrix (ECM) [18]. Similarly, interactions of the force dipoles exerted by cells interacting with their environment can account for simple linear arrangements [19]. Growth can be correlated with changes in local stresses; in conjunction with surface tension, this phenomenon can account for the development of nodules and branchlike structures [17,20]. Mechanical buckling of layers of cells on a soft substrate can introduce periodic patterning in organs [21]. Thus energy gradients among cells offer alternative patternforming mechanisms to the reaction–diffusion phenomenon.
This paper offers a third mechanism for pattern formation, in which patterns are formed neither by gradients in chemical signals nor energy gradients. The patterning arises from changes in cell motility cued by relative cell motions. Straincued motility derives from the cell's living capacity for initiating and sustaining motion, rather than its passive or inanimate movement under an external driving force or energy gradient.
We hypothesize that patterns can form in cell populations if the migration of an individual cell follows rules that depend on the local cell density and the velocity of the cell relative to its contacting neighbours. The rules balance a tendency for cells to enhance relative sliding motion against a tendency to maintain uniform cell–cell separation. If the balance tips in favour of sliding, a homogeneous population will spontaneously move into layered patterns, which are stable once formed.
We test the hypothesis by simulating the migration of populations of 10^{5} secretorystage ameloblast cells in the mouse incisor. The simulations predict the migratory paths of each individual cell and thus patterns in the morphology of the enamel they form. Because the migratory paths of ameloblasts are permanently recorded in mineralized enamel as distinct ‘prisms’ or ‘rods’ [22–28], prediction and nature can be compared in more detail for amelogenesis than for any other cell migration that occurs in a natural context. Exploiting only a single scalar parameter that is not known a priori (with other parameters in the theory calibrated independently), the theory reproduces several characteristics of the morphology of mouse incisor enamel.
2. Idealization of secretorystage ameloblasts
When secretorystage ameloblasts form enamel, they migrate from the dentine–enamel junction (DEJ) in a curved, expanding sheet. The histological section of figure 1 illustrates the process in a human molar; the geometry and some details differ in the mouse incisor, but this image illustrates generic aspects of amelogenesis pertinent to this work.
Ameloblasts are elongated in shape, with length approximately 100 μm and diameter 3.5 μm in the mouse incisor [22,30]. The cells remain approximately invariant in shape during secretorystage migration, at least over shorttime frames, and oriented with their long axis normal to the population sheet and approximately parallel to their neighbours (figure 2a). Secretorystage ameloblasts do not possess lamellipodia or filipodia.
As cells migrate, Tomes’ process, a feature at the tail of the cell, secretes the amelogenin scaffold that is the template for a prism (figure 2a). The prisms, therefore, record the paths followed by Tomes' processes for all cells, which parallel the motion of the cells' centre of mass. The centre of mass of cell i is denoted by x_{i}. The population contains N cells. The surface containing all centres of mass will be termed the growth front, Ω.
The sheet of cells migrates away from the current growth front at speed v_{f}. Simultaneously, individual cells can execute lateral motion (in the plane of the growth front), during which slip between contacting neighbours is assumed to occur approximately uniformly along their zone of contact (figure 2).
A cell's migratory trajectory is represented by a onedimensional locus in threedimensional space, x_{i}(t), where t denotes time. The polar angle γ and azimuthal angle θ define the orientation of the trajectory relative to the growth front (figure 2c).
Even though only the single point x_{i} is used to represent an individual cell, the set of all points {x_{i}, i = 1, … , N} contains information about cell shape: the cell crosssections can be estimated by the Voronoi construction (see the electronic supplementary material, appendix SA2). To pursue the kinetic theory, only the shapes of cell crosssections parallel to the growth front need be computed. The Voronoi construction also determines the set of contacting neighbours of any cell, which can change arbitrarily during a simulation owing to the cells’ lateral motions.
The ameloblast sheet remains continuous, with minor exceptions (e.g. enamel tufts [31]), which are neglected. Furthermore, matrix is secreted only very near the lower end of the cell, where Tomes' process is located, so that along almost their entire lengths cells contact directly, with no intervening ECM [22]. Last, ameloblasts remain close to fixed in number once they commence their secretorystage migration, with limited cell death or proliferation. Therefore, stretching or compression of the population parallel to the growth front implies changes in cell density (the number of cells per unit area) and changes in cells’ crosssectional areas and shapes (see the electronic supplementary material, appendix SA3).
3. Amelogenesis and morphology in the mouse incisor
3.1. The commencement front
In mammalian teeth, the first migration of secretorystage ameloblasts occurs at locations where the DEJ has maximum curvature. Other cells commence migration in turn at progressively more remote locations. Thus a ‘commencement front’ propagates across the DEJ separating the domain of ameloblasts that have already begun to secrete amelogenin and migrate away from the DEJ and the domain of ameloblasts that are assembled on the DEJ, but remain stationary and nonsecretory. As a consequence of the progressive commencement of migration, the first enamel formed is wedgeshaped, tapering to zero thickness at the current position of the commencement front (figure 1).
The fact that the first enamel expands across the DEJ as a wedge implies that cells are affected by a shear strain at the very moment they mature from stationary and nonsecretory to migrating and secretory. The hypothesis was introduced in Cox [29] that the shear strain propagates as a pulse through the cell population and that its arrival at a cell triggers the cell's transition to the secretory stage. The speed of the shear strain pulse is determined in the kinetic theory by the rules of response of individual cells to local strain variations. The measured velocity of the commencement front therefore offers a calibration point for the kinetic theory (see below).
The incisor in a rodent erupts continuously during the life of the animal (figure 3). In the mature animal, the entire spatial pattern formed by the different stages of the odontogenic process achieves an approximately steadystate configuration, translating in a selfsimilar manner back along the erupting tooth towards the tooth organ. The commencement front takes the form of an arc traversing the tooth and progressing back along the tooth, with a velocity equal in magnitude to the velocity of eruption, but opposite in sign.
3.2. Decussation
Prisms in most mammalian enamel tend to be organized in layers of different orientation, a patterning known as decussation (the forming of X patterns; [22]). Decussation is key to fracture toughness and therefore survival. Its presence implies patterning in the motions of migrating cells. In the mouse incisor, cells arrange themselves into layers a single cell thick, transverse to the long axis of the incisor. Cells within a layer move in parallel, but alternate layers migrate to the right or left, forming the X pattern (figure 4). Thus alternate layers of cells shear past one another during migration. The shear strain rates have the very low order of magnitude 10^{−5} s^{−1} (see the electronic supplementary material, appendix SA4).
Towards the end of amelogenesis in mammals, a transition occurs from socalled ‘inner’ enamel, marked by decussation, to ‘outer’ enamel, where prisms lie parallel to one another (figure 4), implying that each cell is now migrating with a fixed set of neighbours. This migration mode will be termed ‘invariant ordering’. One test of the theory is its ability to reproduce decussation and the transition to invariant ordering.
Ameloblast trajectories (prisms) follow wavy loci, with the longest observed wavelength typically comparable with the final thickness of the enamel [33]. The waviness is another aspect of morphology accountable by the kinetic theory [29].
4. Formulation of simulations
4.1. Governing differential equations
We use a discretecell representation of ameloblast migration introduced in prior work [29]. Simulations begin with a cell population assembled on a prescribed DEJ, which is a doubly curved, convexoutward twodimensional surface in threedimensional space (figure 5). The radii of curvature ρ_{x} and ρ_{y} (coordinate system of figures 2 and 3) are determined from images of teeth. The cells are assembled in an approximately hexagonal closepacked array, consistent with observations (C. E. Smith 2012, personal communication).
Cells begin to migrate when passed by the commencement front, which moves at velocity v_{p}. The initial velocity assigned to each cell has magnitude v_{f}, the rate at which the growth front expands away from the DEJ, and direction normal to the DEJ. Thus, the population is entirely homogeneous when migration begins—the initial velocities contain no pattern.
As migration of a cell proceeds, it acquires a lateral velocity component within the plane of the growth front. The lateral velocity changes continuously and must be computed; it reflects the outcome of changes in a cell's motility owing to its changing strain environment. Cell trajectories are recorded and constitute the primary output of a simulation.
Calculations track the ameloblast population through a series of growth front surfaces, Ω_{m}, m = 0, 1, … , separated by the time interval δt (see the electronic supplementary material, appendix SB1). Each growth front is a known surface, derived from the known shape of the DEJ, the motion of the commencement front and the front velocity v_{f}, which are all prescribed in advance. The perpendicular distance of the growth front from the DEJ, at any point equals the thickness of the accumulated enamel, v_{f}(t – t_{0}), where t_{0} is the time at which the commencement front passed that point (see figure 5 and electronic supplementary material, appendices SA2 and SB1). Where the commencement front has already passed, growth fronts are separated by δs = v_{f}δt.
The problem consists of a system of coupled nonlinear differential equations to determine the set of migration trajectories {x_{i}(t), i = 1, … , N}, with one equation for each cell. A cell's velocity v_{i} = dx_{i}/dt varies in proportion to a vector stimulus function S_{i}: 4.1where s is the distance advanced by the growth front, and S_{i} is a function of the relative positions of the contacting neighbours of the ith cell, x_{j} – x_{i}, j = 1, … , n_{i}, with n_{i} the current number of neighbours, and possibly the relative velocities v_{j} – v_{i}.
Possible changes to the set of contacting neighbours of cell i are critical to the value of S_{i}. Therefore, as migration proceeds, the map of contacting neighbours is recomputed using the Voronoi procedure at each new growth front, so that arbitrary changes in cell–cell contact can be permitted. In the current theory, by the definition of S_{i}, a cell is influenced only by its contacting neighbours.
While S_{i} is expressed for computational convenience in terms of cell positions and velocities, this is equivalent to stating a dependence of S_{i} on the local strain and strainrate environment of cell i. Local strain fields in the plane of the growth front can be computed from the set of all cell positions (see the electronic supplementary material, appendix SA3). Tracking the evolving local strain environment of individual cells within the migrating population is crucial to predicting patterning due to straincued motility.
The formulation of the theory is simple, but numerical implementation presents some challenges (see the electronic supplementary material, appendix SB).
4.2. Constraints, cell–cell adhesion and history effects
Cells are constrained to lie on the current growth front, forming a smooth twodimensional curved sheet, rather than being able to spread into a threedimensional assembly. This mimics the universal observation for secretorystage ameloblasts, which always remain on the smooth surface of the enamel formed so far [30]. If this is not enforced, the cells’ response to stimulus usually leads to nonphysical instability (loss of integrity of the growth front), suggesting that the constraining effect of the enamel surface plays an important role in amelogenesis.
Limits are imposed to the lateral velocity and acceleration of any cell (see the electronic supplementary material, appendix SA7). The limits are reached only by a few outliers and do not affect patterning. However, without imposed limits, outliers can cause numerical instability.
Optional elaborations include a historydependent constraint on changes in the velocity of an individual cell and a tendency for pairs of cells to move in parallel if they have been in contact for a long time (cell adhesion) (electronic supplementary material, appendices SA5 and SA6). These phenomena affect pattern stability but do not drive pattern formation.
4.3. Mass inertia and stress
The theory does not consider either mass inertia or mechanical force or stress. Kinetic energy associated with cell motion is approximately 15 orders of magnitude lower than that associated with strain (the cells move very slowly) [29].
Computation of stress is unnecessary. The key concept, patterning due to straincued motility, arises as a kinetic instability and can be described entirely in terms of cell motions and rates of motion. Further remarks on force and stress appear in the electronic supplementary material, appendix SD.
4.4. Solution of an inverse problem
Given sufficient knowledge of the internal workings of a single cell, stimulus rules S_{i} could be prescribed a priori and used in a ‘bottomup’ approach to define the parameters of a multicell simulation and thence predict the outcome of organogenesis. An alternative approach, followed here, solves the inverse problem: begin with the observed characteristics of a mature organ as input data and ask the question, what singlecell response function when assigned to each cell in a large population will lead to the population reproducing the observations? Obtaining a solution to the inverse problem proves that straincued motility is a feasible mechanism for pattern formation.
5. Correlation of wavy microstructure and the ‘commencement front’
In prior work [29], a stimulus S_{i} was introduced that depended linearly on fluctuations in the separations d_{ij} = x_{j} – x_{i} of cell i from each of its contacting neighbours: 5.1aand 5.1bwhere (d_{ij} – d)/d is the average strain between the two cells' centres; is a unit vector from cell i to neighbour j; and w_{ij} is a weight factor that is proportional to the computed contact area between cell i and neighbour j (see the electronic supplementary material, appendix SA). (A simple example of the evolution of the neighbour map and w_{ij} appears in the electronic supplementary material, appendix SE.) The linearity of equation (5.1) was justified by reanalysing experimental observations of how changes in the cytoskeleton of aortic cells correlate with the magnitude of applied strain, as well as by modelbased interpretation of enamel waviness [29,34].
The tendency of the rule of equation (5.1) to restore a preferred cell centretocentre separation combined with ratebased ‘inertia’ generates wave propagation in the population (see the electronic supplementary material movie, ‘pulse 9.7.10’; [29]). The factor in equation (4.1) governs the rate of propagation of strain disturbances, v_{p}.
All constants in the linearstimulus theory reduce by scaling invariances to a single adjustable parameter, the ratio v_{p}/v_{f} (see the electronic supplementary material, appendix SA2). Other parameters in the theory describe the geometry of the DEJ and are therefore known a priori.
The observed wavelength of wavy prism loci in the mouse incisor is λ = 100–170 μm [29]. The predicted wavelength is determined by the ratio of rates, v_{p}/v_{f}, and the width, W, of the tooth measured along the DEJ. Numerical solutions yield 5.2for the linearstimulus rule of equation (5.1) [29].
If the zone of odontogenesis in a mouse incisor (figure 3) remains stationary relative to the tooth organ, the whole zone, including the commencement front for secretorystage ameloblast migration, must move with speed equal to the observed rate of eruption of the incisor. Following the hypothesis that the commencement front is governed by a propagating strain pulse [29], the rate of eruption should equal v_{p}. Using eruption data for a mature wildtype mouse [35], one infers v_{p} = 250 ± 40 μm d^{−1}.
The observed migration velocity for ameloblasts in the mouse incisor is v_{f} ≈ 15 ± 2 μm d^{−1} (C. E. Smith 2010, personal communication). Thus v_{p}/v_{f} ≈ 16 ± 2. The DEJ has an arc length W = 0.6–0.8 mm (S. N. White & D. E. Smith 2012, personal communications; [36]). Thus equation (5.2) predicts a wavelength λ = 60–90 μm, which is close to the measured wavelength.
The hypothesis that the commencement of secretory migration is triggered by a strain pulse correlates waviness and the speed of the commencement front to within a factor of 1.5 in the mouse incisor and even better in the human molar ([37] and work to be published).
6. Simulations of the mouse incisor
A segment of the mouse incisor is represented by a curved rectangular strip, with a commencement front propagating along its length (figure 5). Simulations considered up to 100 000 cells, occupying an area 0.25–0.35 mm wide and 2–3 mm long. Choosing a long simulated domain allows patterns to develop fully after initial transient behaviour. The simulations are subscale, the simulated width being approximately one half the width of the incisor, but tests for simulations of different size showed that this does not affect pattern formation. For most simulations, but not all, the initial hexagonal array of cells is oriented so that rows of cells with separation d lie across the incisor. All cells in such a row are reached simultaneously by a straight commencement front propagating along the incisor.
6.1. Shearstimulus conjecture
Studies based on the linearstimulus rule of equation (5.1) establish the importance of the rate ratio, v_{p}/v_{f}, relate the dominant wavelength to the response of an individual cell, and show that the wavelength is almost independent of the initial curvature of the DEJ [29]. However, other interesting features of the patterns formed during amelogenesis are not obtained. Cells responding to the linearstimulus law exhibit invariant ordering if they migrate from a DEJ possessing only modest curvature: their trajectories are approximately parallel, and each cell retains the same contacting neighbours throughout the migration. If the curvature increases, the theory predicts the onset of ordering instability whereupon invariance is lost, but specific patterns found in nature do not arise [29].
The stimulus rule was therefore modified to allow differentiated treatment of shear strains. The relative velocity of cell i and its neighbour j, v_{ij} = v_{j} – v_{i}, was decomposed into components acting along the vector between the two cells’ centres and perpendicular to . Denote the latter . Nonzero implies that either the interiors of cells i and j experience nonzero shear strain rate or their contacting surfaces slip past one another (see the electronic supplementary material, figure SA1). The sum of these effects accounts for ; their relative proportions need not be determined to test the theory.
Rules were sought that reproduce decussation. Some trial rules are summarized in the electronic supplementary material, appendix SA1. The two most instructive, one failed and one successful, are as follows.
6.1.1. Easy slip model
The shearing motion of decussation suggests that a mechanism for its creation might be planes of easy slip, analogous to dislocation slip planes in crystal plasticity. The pertinent slip for ameloblasts would be in their lateral motion across the growth front.
One trial therefore reduced the stimulus S_{ij} exerted by neighbour j on cell i so that the pair can move relatively freely past one another when is large. This modification reduces resistance to slip by specifying that stimulus no longer tends to restore the equilibrium separation between two cells if their relative motion is a sliding motion; a cell's passive response to shear is modified to facilitate its tendency to slide. The modification generates qualitative changes in the patterns formed by a population, but local ordering of cells into groups is transitory and the stable and regulated structure of figure 4 is not achieved.
6.1.2. Stimulated shear model
The successful trial added a second stimulus factor that actively promotes further sliding motion of cell i past a neighbour j when is already nonzero: 6.1awhere w_{ij} is the weight factor of equation (5.1) and 6.1bwith c_{1} and c_{2} dimensionless scalar parameters. The shear stimulus owing to the jth neighbour acts in the direction of (see the electronic supplementary material, figure SA1). The net shear stimulus on cell i is the vector sum of the stimuli from all neighbours and thus indicates the net shear velocity gradient around an individual cell.
Equation (6.1) describes more than easy slip that accommodates an applied shear strain passively; it implies that cells exert themselves to accelerate past one another, with possibly rapid motion (as far as ameloblasts can be rapid!) when stimulated. Thus the cells exhibit enhanced motility under shear stimulus.
Parameter c_{1} sets the magnitude of the shear stimulus. The exponential factor (parameter c_{2}) limits the stimulus owing to large sliding velocities. Solutions that successfully mimic nature are presented below for c_{1} > 0, but solutions for c_{1} < 0 are also interesting.
6.2. Calibration of the theory
Simulations are controlled by a small number of parameters that relate to geometry, the response of a cell to strain stimulus, or constraints on cell motions.
The geometrical parameters are the two radii of curvature of the DEJ, ρ_{x} and ρ_{y}, and the preferred separation of cells, d. For the mouse incisor: ρ_{x} = 6.7 mm, ρ_{y} = 0.3 mm and d = 3.5 μm.
The response of a cell to stimulus depends on three parameters: the ratio v_{p}/v_{f} and the parameters c_{1} and c_{2} of equation (6.1). From data for the mouse incisor (§5), v_{f} = 13 μm d^{−1} and v_{p} = 250 ± 40 μm d^{−1} [35]. The unknown parameters related to strain stimulus are c_{1} and c_{2}. Simulations show that only c_{1} has a significant effect on symmetry breaking and the initiation of patterns.
The parameters related to constraints on cell motion describe the strength of history effects or cell–cell adhesion, or limit the maximum lateral velocity or lateral acceleration. Their values, listed in the electronic supplementary material, appendix SA7, affect only the stability of patterns.
6.3. Decussation
The shearstimulus rule of equation (6.1) successfully creates decussation: cells spontaneously form layers with the direction of the lateral velocity alternating from layer to layer. At the beginning of the migration, the shearstimulus effect on cell motility, governed by c_{1}, tends to amplify any small local shear gradient, potentially leading to large amplitude velocity excursions. Local shear gradients can arise, for example, because of the curvature of the DEJ (initial cell velocities are set normal to the DEJ), the motion of the commencement front or small fluctuations in the velocities of individual cells. The cells tendency to relax to the preferred density, at a rate governed by v_{p}/v_{f}, counters the disrupting effect of the shear stimulus, tending to keep cells in their initial hexagonal pattern, with invariant neighbour sets. The outcome of this competition is settled very quickly (figures 6 and 7); for fixed c_{1} and v_{p}/v_{f}, decussation (or precursors of decussation) either occurs within a few microns of the DEJ or not at all.
The layered patterning has no explicit origin in equation (6.1), which is applied isotropically (without orientational bias) to each individual cell in the population and its neighbours; nor in the initial conditions, which define a spatially homogeneous cell population, with cells close to equally spaced at their preferred separation and each cell migrating at the same speed perpendicularly away from the DEJ. Patterning arises because the stimulated shear response generates a kinetic instability in the cell population and a tendency to favour layered ordering.
A typical simulation shows a transient regime, characterized by incoherent attempts for local cell groups to develop patterns, emerging into a universal stable pattern of decussation that spans the width of the tooth. The transient regime is confined spatially to the domain where cells first begin to migrate in the simulation (details in the electronic supplementary material, appendix SB3). Stable decussation appears as a steadystate solution in the wake of the commencement front as it leaves the transient domain, and, once established, continues indefinitely along the tooth. In the steadystate region, the initial cell motions and the selforganization of cells into decussation develop a system that translates invariantly along the incisor.
Typical steadystate decussation is illustrated in figure 6. These discretecell section images show the Voronoi polygons representing the crosssections of individual cells in a patch of approximately 3400 cells extracted from a simulation of 97 000 cells. The commencement front is moving towards the top of the figure; cells above the line are not yet migrating while those below have migrated a distance out of the page that is proportional to their distance from the commencement front. Each cell is colourcoded to show the azimuthal and polar angles θ and γ of its trajectory (figure 2c). For established decussation, γ will be termed the ‘decussation angle’ (figure 4). Cells that are still stationary are assigned angle γ = 0, because by assumption they initially move normal to the front. Figure 6a shows that lateral motion large enough to be observable in histological images (say, γ > 5°) sets in approximately 50 μm behind the commencement front, when the front has advanced approximately 3 μm from the DEJ. Prior to that, the cells have already begun to establish themselves in decussating layers: small lateral motions create systematic variations in θ (figure 6b). (The apparent twodomain pattern in the values of θ above the commencement front is not significant; it arises from the treatment of indeterminate values when γ = 0.)
In the example of figure 6, decussation nucleates at two positions across the front, with a phase incompatibility between the two: observe that rows in the circled region in figure 6b have colours that differ from their extensions into the region of the furthest advanced patterning. This patterning error is relieved by further development, which can be seen in the electronic supplementary material, movie ‘incipient decussation (vor 0.06 48 292 defectfree 5.28.12).wmv’ (see the electronic supplementary material, appendix SF).
The decussation pattern formed on a plane transverse to the tooth after substantial migration of the cell population away from the DEJ is illustrated in the trajectory slice of figure 7, which should be compared directly with figure 4. The electronic supplementary material, movie ‘decussating cells (y.2 b 4.22.12).wmv’ highlights the lateral motions of a few cells in steadystate decussation and their repeated swapping of neighbours.
The simulations reproduce nature in that decussation begins very close to the DEJ (3 μm), layers with alternating direction of lateral motion are a singlecell thick and the layers are oriented transverse to the long axis of the incisor (but see also the electronic supplementary material, appendix SB9).
None of these characteristics is explicit in equation (6.1). The shear response applies isotropically to any cell until symmetry is broken and the layer thickness is not specified: the cells achieve singlelayer thickness autonomously.
When decussation is established, the angle θ for a given cell takes a value very near either 0° or 180°, while, for the case of figures 6 and 7, the angle γ attains an average value of approximately 40°, about which systematic oscillations arise.
The plot for γ (figure 6a) shows waves of variation in the orientation of a cell's trajectory. The waves correlate with distortions in the cell shapes derived from the Voronoi procedure as cells swap neighbours, which are visible in both figure 6a,b. Far from a free edge, distortions in cell shapes tend to occur along straight lines, as an entire row of cells switches simultaneously to new sets of neighbours. Near a free edge, distortion patterns become sinuous; as discussed below, the commencement front is delayed near the free edge, which causes delays in neighbour switching.
The left boundary of figure 6a,b is a boundary of the population. Disorderly cell behaviour arises at the boundary in the simulations, because the shear stimulus drives cells out of the population, whereupon their motions are confused.
6.4. Predicted shape of the commencement front
The predicted development of lateral cell motions is delayed near the boundary (left edge) in figure 6. Cells relax strains more easily near a free edge, because they are not entirely surrounded by neighbours. This delays the buildup of large shear velocities and thus decussation.
Edge retardation is also seen in the time derivative of von Mises’ equivalent strain (figure 8a), which is approximately proportional to the magnitude of shear velocities (see the electronic supplementary material, appendix SA3). The commencement front was stipulated to be straight and is therefore not the source of edge retardation.
The predictions of figure 8a have an interesting implication for the shape of the commencement front. If, instead of assuming a straight commencement front, migration was assumed to be triggered by the attainment of a critical value of the shear strain, the commencement front should follow the C shape of figure 8a. This is indeed the case in nature (figure 8b).
6.5. Roles of c_{1} and v_{p}/v_{f}: transitions in and out of decussation
Whether decussation or invariant ordering arises depends on the values of c_{1} and v_{p}/v_{f}. In a population of 10^{5} cells moving independently and interacting in a strongly nonlinear manner, the transition from invariant ordering to decussation is not easily predicted by analysis. Therefore, it was studied by inspecting the outcome of numerous simulations, assessed away from any transient regime (see the electronic supplementary material, appendix SB4). The numerical survey showed that for the rules of equations (5.1) and (6.1), decussation only occurs if c_{1} > 0.5 (figure 9). On the other hand, if the shear stimulus is too strong (typically c_{1} > 1.0), orderly decussation gives way to noisy patterns. Therefore, relevant values of c_{1} lie in the fairly narrow interval 0.5 ≤ c_{1} ≤ 1.0.
Recall that v_{p}/v_{f} = 16 (§5), with uncertainty of approximately 50 per cent. Given 0.5 ≤ c_{1} ≤ 1 and 10 ≤ v_{p}/v_{f} ≤ 20, figure 9 predicts that any combination of v_{p}/v_{f} and c_{1} that produces decussation will lie not far from the critical condition. The possibility exists of crossing the critical condition at some time during amelogenesis.
The possibility of a transition is strengthened by the fact that the wavelength observed in the enamel of most mammalian species increases substantially from the DEJ to the outer enamel surface, typically by a factor of 1.5–3 [33]. This observation is at least partially correlated with changes in migration velocity [37]. An increase by a factor more than 2 is observed in the mouse incisor [37]. The implied decrease in v_{p}/v_{f} is substantial on the transition map of figure 9.
However, decreasing v_{p}/v_{f} alone proves insufficient to change a pattern once decussation is already established in a simulation; and the pattern is also maintained if the shearstimulus parameter c_{1} is decreased well below the critical value needed for its nucleation. The conditions for maintaining decussation are milder than those for its nucleation. However, if c_{1} is decreased to values of approximately 0.2 in midsimulation, a spontaneous switch from decussation to invariant ordering does occur. In the simulation of figure 7, the parameter change was made over the whole simulation simultaneously, resulting in rapid reversion to invariant ordering (motion along parallel paths with fixed neighbours).
A transition similar in appearance to figure 7 is common in mammalian enamel.
Other than its role in determining the transition from invariant ordering to decussation, the parameter v_{p}/v_{f} governs the rate of propagation of strain disturbances across the population (see the electronic supplementary material, movie ‘pulse 9.7.10’; [29]). The other main influence of parameter c_{1} is in determining the decussation angle, as follows.
6.6. Predictions of the decussation angle γ
The decussation angle γ depends mainly on the shearstimulus parameter c_{1}. When c_{1} is varied mildly with decussation already established, decussation is preserved but with changing γ (figure 10).
Values of γ range from approximately 45° down to zero as c_{1} decreases from approximately 1 to 0.2 (the threshold for reversion to invariant ordering). However, if c_{1} is restricted to the values necessary to nucleate decussation, i.e. c_{1} ≥ 0.5, predicted values of γ are in the range 35°–55° (figure 10). In the mouse incisor, γ ≈ 35° (figure 4).
With c_{1} > 0.5 as required to nucleate decussation, γ should lie in the range 35°–55° early in migration and remain in this range unless c_{1} decreases. A decrease in c_{1} is of course inferred to occur, at least approaching the transition from inner to outer enamel. Measured variations of γ during migration, such as those visible in figure 4, could be a source of information about time variations in c_{1}.
7. Further aspects of the theory
ESM appendices contain:

— numerical methods, including the idealized shape of the DEJ, the definition of incremented growth fronts, the integration of differential equations, neighbour mapping, the definition and computation of strains, the definition of stimulus functions, history and cell–cell adhesion models, noiselimiting bounds and visual rendering of output;

— details particular to the mouse incisor, including model calibration, illustrations of the variety in predicted ordering patterns, illustrations of stability and transient behaviour, mesh tests, prediction of the decussation angle γ and discussion of the sources and magnitude of uncertainty in predictions;

— discussion of the relation of the present theory to a recent theory of network formation, in which a leading factor is also the ratio of cell migration rate to cell relaxation rate;

— remarks on mechanical forces and mechanical stress;

— analysis of the kinetics of decussation in a simplified limiting case, illustrating the changes in neighbour separations and relative velocities as rows of cells slide past one another, and the accompanying changes in neighbour sets and variations of the area of cell–cell contact. The limiting case also relates the computational model to the sine–Gordon equation; and

— movies.
8. Concluding remarks
Patterning via straincued motility is a new concept for pattern formation among migrating cells. Kinetic effects can break symmetry in a cell population that is initially homogeneous and governed by rules that do not explicitly prescribe any pattern. In the case studied, ameloblasts in the mouse incisor spontaneously form layered patterns. Once the decussation pattern is formed over the width of the tooth, it tends to be stabilized by cell–cell adhesion and history effects. Pattern stability in the presence of perturbation promotes robust organogenesis.
An inverse problem is solved and therefore the solution cannot be assured to be unique; a modified motility rule might also work. However, the existence of one solution proves that patterning via straincued motility is a viable mechanism for creating enamel morphology.
The rules for cell behaviour must contain a linearstimulus response to reproduce the wavy character seen in the enamel of the mouse and other animals. They must also contain a shearstimulus response in which sliding motion is enhanced to be able to generate decussation. Enhanced shear response implies living action, powered by the cells’ capacity for work, not inanimate response to energy gradients.
A discretecell model, in contrast to a continuum model, allows realistic depiction of symmetry breaking that originates at the scale of the single cell, acting as an automaton moving in the field of its neighbours. A large population must be considered, because particular patterns arise as collective cell behaviour.
The simulations reproduce various aspects of the morphology of enamel in the mouse incisor, correlating waviness with the speed of the commencement front, predicting the shape of the commencement front, the layering pattern and angle of decussation, and providing a mechanism for the transition from inner to outer enamel. Thus essentially all morphology at scales greater than the prism (cell) diameter can be predicted as resulting from strain stimuli, without reference to chemical stimuli.
While chemical stimuli do not appear explicitly in the theory, they are implicit in the response functions for the single cell. However, if chemical signals are involved in the process of pattern formation, it is enough that they aid in the expression of straincued motility. The capacity of a cell to exhibit the required motility change in response to strain can be assumed to be genetically encoded.
The two most influential parameters in the rules for cell response are v_{p}/v_{f}, which has been independently calibrated and is therefore unavailable for fitting data, and the shearstimulus strength, c_{1}. Thus the various aspects of morphology are predicted by manipulating a single adjustable scalar parameter.
The mechanism of straincued motility can trigger the transition in ordering from inner to outer enamel under smoothly changing rate factors, because of the existence of the transition boundary (figure 9). Once triggered, the transition itself could trigger gene expression and changes in the concentration of chemical species.
The fact that the rate of eruption of the mouse incisor and the rate of propagation of strain pulses are correlated introduces the possibility that the rate of eruption is governed by the propagating strain pulse associated with the commencement front.
In the formation of most other organs, cell proliferation, differentiation and death are all involved, whereas the theory of amelogenesis treats migration only. Nevertheless, kinetic effects can have a major role in pattern development in other organs [5,17,20]. Organs whose development begins with invagination of a mesenchymal–epithelial bilayer to create an expanding doubly curved process immediately possess the characteristics needed to trigger straincued ordering transitions.
The theory of patterning via straincued motility relates observed morphology to the response of a single cell to strain stimulus, which in turn implies information about the biochemical functionality of the cell. Thus the theory could provide a key for interpreting the patterns in fossilized enamel in terms of the response characteristics of cells at that epoch and may aid the identification of evolutionary steps in cell behaviour.
Acknowledgements
The author is indebted to Drs S. N. White, M. L. Snead, C. E. Smith, R. Lacruz, P. Gilbert and T. Bromage, for many remarks on the known characteristics of amelogenesis and enamel morphology, and correction of errors; and to several anonymous reviewers for very beneficial feedback.
 Received March 24, 2013.
 Accepted April 2, 2013.
 © 2013 The Author(s) Published by the Royal Society. All rights reserved.