How the tooth got its stripes: patterning via strain-cued motility

Brian N. Cox


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. Discrete-cell 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 secretory-stage ameloblasts from those that are not yet migrating, the initiation and sustainment of layered, fracture-resistant decussation patterns (cross-plied microstructure) and the transition from decussating inner enamel to non-decussating 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 symmetry-breaking mechanism underlying pattern formation. Many experiments have shown that strain is a factor in cell migration, proliferation and differentiation [311]. 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 cell-adhesion energies differ [1315]. 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 branch-like 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 pattern-forming 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. Strain-cued 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 105 secretory-stage 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’ [2228], 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 secretory-stage ameloblasts

When secretory-stage 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.

Figure 1.

The migration of an ameloblast population as a continuous curved sheet of cells that expands away from the DEJ as enamel forms underneath is typified by this image of the early bell stage of odontogenesis in the human molar (adapted by Cox [29] from Nanci [30]). Odontoblasts migrate down to form dentine and ameloblasts migrate up, with a thin boomerang-shaped wedge of enamel (black) already formed in their wake. Axes (x, z) are global coordinates.

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 secretory-stage migration, at least over short-time frames, and oriented with their long axis normal to the population sheet and approximately parallel to their neighbours (figure 2a). Secretory-stage ameloblasts do not possess lamellipodia or filipodia.

Figure 2.

Schematic of ameloblasts arrayed on the growth front: (a) view showing length of cells perpendicular to growth front and (b) view from above showing approximately equi-axed cross-sections of cells. Advance normal to the growth front and lateral motion are exemplified for two cells. (c) The angles θ and γ define the orientation of a trajectory relative to the growth front. Coordinates (X1, X2 and X3) are aligned locally with growth front.

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 xi. 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 vf. 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 one-dimensional locus in three-dimensional space, xi(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 xi is used to represent an individual cell, the set of all points {xi, i = 1, … , N} contains information about cell shape: the cell cross-sections can be estimated by the Voronoi construction (see the electronic supplementary material, appendix SA2). To pursue the kinetic theory, only the shapes of cell cross-sections 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 secretory-stage 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’ cross-sectional 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 secretory-stage 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 non-secretory. As a consequence of the progressive commencement of migration, the first enamel formed is wedge-shaped, 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 non-secretory 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 steady-state configuration, translating in a self-similar 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.

Figure 3.

(a) Enamel-bearing surface of mouse incisor. The curves' ‘commencement front’ and ‘transition boundary’ between inner and outer enamel recreate the dash-dot curve of fig. 22 and dotted curve of fig. 23 in experimental data of Warshawsky & Smith [32]. (b) Shape of enamel layer on transverse section (perpendicular to long axis of tooth).

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).

Figure 4.

Nature: enamel morphology in the mouse incisor (wild-type), viewed on a section transverse to the axis of the tooth and inclined to coincide with prism layers. Prisms appear as whitish lines. Part of one prism layer is visible behind another (area marked ‘decussation’). Beyond the curve marked ‘transition,’ prisms switch to a pattern of invariant ordering. Image courtesy Dr Shane White.

Towards the end of amelogenesis in mammals, a transition occurs from so-called ‘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 discrete-cell 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, convex-outward two-dimensional surface in three-dimensional 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 close-packed array, consistent with observations (C. E. Smith 2012, personal communication).

Figure 5.

Typical shape of the DEJ for simulating the mouse incisor, showing the growth front at one time. The growth front meets the DEJ at the current location of the commencement front. Successive growth fronts are increasingly distant from the DEJ. All dimensions in millimetres. For ρx = 0.3 mm.

Cells begin to migrate when passed by the commencement front, which moves at velocity vp. The initial velocity assigned to each cell has magnitude vf, 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 vf, 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, vf(tt0), where t0 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 = vfδt.

The problem consists of a system of coupled nonlinear differential equations to determine the set of migration trajectories {xi(t), i = 1, … , N}, with one equation for each cell. A cell's velocity vi = dxi/dt varies in proportion to a vector stimulus function Si:Embedded Image 4.1where s is the distance advanced by the growth front, and Si is a function of the relative positions of the contacting neighbours of the ith cell, xjxi, j = 1, … , ni, with ni the current number of neighbours, and possibly the relative velocities vjvi.

Possible changes to the set of contacting neighbours of cell i are critical to the value of Si. Therefore, as migration proceeds, the map of contacting neighbours is re-computed 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 Si, a cell is influenced only by its contacting neighbours.

While Si is expressed for computational convenience in terms of cell positions and velocities, this is equivalent to stating a dependence of Si on the local strain and strain-rate 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 strain-cued 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 two-dimensional curved sheet, rather than being able to spread into a three-dimensional assembly. This mimics the universal observation for secretory-stage 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 non-physical 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 history-dependent 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 strain-cued 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 Si could be prescribed a priori and used in a ‘bottom-up’ approach to define the parameters of a multi-cell 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 single-cell 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 strain-cued motility is a feasible mechanism for pattern formation.

5. Correlation of wavy microstructure and the ‘commencement front’

In prior work [29], a stimulus Si was introduced that depended linearly on fluctuations in the separations dij = |xjxi| of cell i from each of its contacting neighbours:Embedded Image 5.1aandEmbedded Image 5.1bwhere (dijd)/d is the average strain between the two cells' centres; Embedded Image is a unit vector from cell i to neighbour j; and wij 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 wij appears in the electronic supplementary material, appendix SE.) The linearity of equation (5.1) was justified by re-analysing experimental observations of how changes in the cytoskeleton of aortic cells correlate with the magnitude of applied strain, as well as by model-based interpretation of enamel waviness [29,34].

The tendency of the rule of equation (5.1) to restore a preferred cell centre-to-centre separation combined with rate-based ‘inertia’ generates wave propagation in the population (see the electronic supplementary material movie, ‘pulse 9.7.10’; [29]). The factor Embedded Image in equation (4.1) governs the rate of propagation of strain disturbances, vp.

All constants in the linear-stimulus theory reduce by scaling invariances to a single adjustable parameter, the ratio vp/vf (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, vp/vf, and the width, W, of the tooth measured along the DEJ. Numerical solutions yieldEmbedded Image 5.2for the linear-stimulus 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 secretory-stage 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 vp. Using eruption data for a mature wild-type mouse [35], one infers vp = 250 ± 40 μm d−1.

The observed migration velocity for ameloblasts in the mouse incisor is vf ≈ 15 ± 2 μm d−1 (C. E. Smith 2010, personal communication). Thus vp/vf ≈ 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 sub-scale, 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. Shear-stimulus conjecture

Studies based on the linear-stimulus rule of equation (5.1) establish the importance of the rate ratio, vp/vf, 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 linear-stimulus 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, vij = vjvi, was decomposed into components acting along the vector Embedded Image between the two cells’ centres and perpendicular to Embedded Image. Denote the latter Embedded Image. Non-zero Embedded Image implies that either the interiors of cells i and j experience non-zero 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 Embedded Image; 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 Sij exerted by neighbour j on cell i so that the pair can move relatively freely past one another when Embedded Image 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 Embedded Image that actively promotes further sliding motion of cell i past a neighbour j when Embedded Image is already non-zero:Embedded Image 6.1awhere wij is the weight factor of equation (5.1) andEmbedded Image 6.1bwith c1 and c2 dimensionless scalar parameters. The shear stimulus owing to the jth neighbour acts in the direction of Embedded Image (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 c1 sets the magnitude of the shear stimulus. The exponential factor (parameter c2) limits the stimulus owing to large sliding velocities. Solutions that successfully mimic nature are presented below for c1 > 0, but solutions for c1 < 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 vp/vf and the parameters c1 and c2 of equation (6.1). From data for the mouse incisor (§5), vf = 13 μm d−1 and vp = 250 ± 40 μm d−1 [35]. The unknown parameters related to strain stimulus are c1 and c2. Simulations show that only c1 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 shear-stimulus 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 shear-stimulus effect on cell motility, governed by c1, 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 vp/vf, 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 c1 and vp/vf, decussation (or precursors of decussation) either occurs within a few microns of the DEJ or not at all.

Figure 6.

Simulation: discrete-cell section images (see the electronic supplementary material, appendix SA8) show the computed shapes of cells on the current growth front forming a decussating pattern, with cells colour-coded according to (a) angle γ (polar angle or ‘decussation angle’) and (b) angle θ (azimuthal angle) (figure 2c). Cells rise with the growth front out of the page, while moving laterally within the page.

Figure 7.

Simulation: A trajectory slice (see the electronic supplementary material, appendix SA8) from a simulation of 56 967 cells, showing the trajectories of cells that originated within a narrow strip on the DEJ. Cells in two rows migrate to the left, those in a third row to the right. When vp/vf and c1 are decreased, the pattern transitions to invariant ordering. Compare with figure 4.

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 steady-state solution in the wake of the commencement front as it leaves the transient domain, and, once established, continues indefinitely along the tooth. In the steady-state region, the initial cell motions and the self-organization of cells into decussation develop a system that translates invariantly along the incisor.

Typical steady-state decussation is illustrated in figure 6. These discrete-cell section images show the Voronoi polygons representing the cross-sections 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 colour-coded 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 two-domain 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 defect-free 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 steady-state 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 single-cell 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 single-layer 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 build-up 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.

Figure 8.

(a) Continuous field plot (see the electronic supplementary material, appendix SA8) shows predicted equivalent strain rate rising behind an assumed straight commencement front (arb. units). (b) Commencement front observed in a rat incisor [32]. (Ignore other boundaries marked in (b).)

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 c1 and vp/vf: transitions in and out of decussation

Whether decussation or invariant ordering arises depends on the values of c1 and vp/vf. In a population of 105 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 c1 > 0.5 (figure 9). On the other hand, if the shear stimulus is too strong (typically c1 > 1.0), orderly decussation gives way to noisy patterns. Therefore, relevant values of c1 lie in the fairly narrow interval 0.5 ≤ c1 ≤ 1.0.

Figure 9.

Domains of decussation and invariant ordering. The critical condition is marked as a band because of the influence of shape and other factors (see the electronic supplementary material, appendix SB4). Points A and B mark independent estimates of vp/vf for the beginning and end of amelogenesis in the mouse incisor, along with the calibration value c1 = 0.7. Point C marks a value of c1 for which transition to invariant ordering will occur when a decussating pattern has already been established. (Online version in colour.)

Recall that vp/vf = 16 (§5), with uncertainty of approximately 50 per cent. Given 0.5 ≤ c1 ≤ 1 and 10 ≤ vp/vf ≤ 20, figure 9 predicts that any combination of vp/vf and c1 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 vp/vf is substantial on the transition map of figure 9.

However, decreasing vp/vf alone proves insufficient to change a pattern once decussation is already established in a simulation; and the pattern is also maintained if the shear-stimulus parameter c1 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 c1 is decreased to values of approximately 0.2 in mid-simulation, 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 vp/vf 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 c1 is in determining the decussation angle, as follows.

6.6. Predictions of the decussation angle γ

The decussation angle γ depends mainly on the shear-stimulus parameter c1. When c1 is varied mildly with decussation already established, decussation is preserved but with changing γ (figure 10).

Figure 10.

When c1 is restricted to values needed to nucleate decussation, predicted values of γ are confined to 35° < γ < 55°.

Values of γ range from approximately 45° down to zero as c1 decreases from approximately 1 to 0.2 (the threshold for reversion to invariant ordering). However, if c1 is restricted to the values necessary to nucleate decussation, i.e. c1 ≥ 0.5, predicted values of γ are in the range 35°–55° (figure 10). In the mouse incisor, γ ≈ 35° (figure 4).

With c1 > 0.5 as required to nucleate decussation, γ should lie in the range 35°–55° early in migration and remain in this range unless c1 decreases. A decrease in c1 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 c1.

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, noise-limiting 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 strain-cued 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 strain-cued motility is a viable mechanism for creating enamel morphology.

The rules for cell behaviour must contain a linear-stimulus response to reproduce the wavy character seen in the enamel of the mouse and other animals. They must also contain a shear-stimulus 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 discrete-cell 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 strain-cued 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 vp/vf, which has been independently calibrated and is therefore unavailable for fitting data, and the shear-stimulus strength, c1. Thus the various aspects of morphology are predicted by manipulating a single adjustable scalar parameter.

The mechanism of strain-cued 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 strain-cued ordering transitions.

The theory of patterning via strain-cued 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.


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.


View Abstract