A model of developmental canalization, applied to human cranial form

Developmental mechanisms that canalize or compensate perturbations of organismal development (targeted or compensatory growth) are widely considered a prerequisite of individual health and the evolution of complex life, but little is known about the nature of these mechanisms. It is even unclear if and how a “target trajectory” of individual development is encoded in the organism’s genetic-developmental system or, instead, emerges as an epiphenomenon. Here we develop a statistical model of developmental canalization based on an extended autoregressive model. We show that under certain assumptions the strength of canalization and the amount of canalized variance in a population can be estimated, or at least approximated, from longitudinal phenotypic measurements, even if the target trajectories are unobserved. We extend this model to multivariate measures and discuss reifications of the ensuing parameter matrix. We apply these approaches to longitudinal geometric morphometric data on human postnatal craniofacial size and shape as well as to the size of the frontal sinuses. Craniofacial size showed strong developmental canalization during the first 5 years of life, leading to a 50% reduction of cross-sectional size variance, followed by a continual increase in variance during puberty. Frontal sinus size, by contrast, did not show any signs of canalization. Total variance of craniofacial shape decreased slightly until about 5 years of age and increased thereafter. However, different features of craniofacial shape showed very different developmental dynamics. Whereas the relative dimensions of the nasopharynx showed strong canalization and a reduction of variance throughout postnatal development, facial orientation continually increased in variance. Some of the signals of canalization may owe to independent variation in developmental timing of cranial components, but our results indicate evolved, partly mechanically induced mechanisms of canalization that ensure properly sized upper airways and facial dimensions.


Introduction
The precision and stability with which animal development brings about the adult phenotype has puzzled generations of scientists and is still not well understood.In the 1940s and 50s, the early developmental and theoretical biologist Conrad H. Waddington postulated developmental mechanisms that "canalize" the phenotype in the presence of genetic and environmental perturbations [1][2][3].A similar concept of auto-regulatory mechanisms in development was proposed independently and at about the same time by the Russian biologist Ivan Schmalhausen.Today, we know that a disruption of developmental canalization due to environmental or genetic perturbations underlies several of the diseases that are typically considered of multifactorial or unclear aetiology and that can show a spectrum of symptoms (e.g., metabolic syndrome, autism spectrum disorder, alcohol spectrum disorder; [4][5][6][7][8][9][10]).A related concept in evolutionary biology is "phenotypic robustness," which refers to the ability of biological systems to produce functional phenotypes in the face of genetic or environmental variation [11][12][13].
In these contexts, the term canalization can refer to phenomena at two different levels, a population pattern and a mechanism in individual development.In population genetics, canalization denotes the reduced effect of an allele, which can lead to a decrease in genetic variation within a population of individuals [13][14][15].The second meaning of canalization, to which we refer as "developmental canalization" and which is the focus of this paper, denotes mechanisms that reduce or compensate deviations of individual development from its target trajectory.Such deviations can arise from both environmental and genetic perturbations (e.g., [16]).Developmental canalization can reduce phenotypic variation in a population and can (but not necessarily does) underly the population genetic concept of canalization if the canalizing mechanisms are subject to allelic variation.Genetic and phenotypic variances can also reduce or increase due to multiple other mechanisms and phenomena (see below and, e.g., [13,[17][18][19]).
Developmental canalization, also referred to as targeted growth, has been studied extensively in human and mammalian postnatal growth [17,18,[20][21][22][23][24].For instance, relatively short individuals tend to undergo accelerated growth, whereas unusually tall individuals show a decreased growth rate to finally achieve an adult body height within the normal range of variation.In human medicine and animal science, these phenomena are also called "catch-up growth" and "catch-down growth," which are often observed after periods of malnutrition or pathological conditions (e.g., [20,24,25]).The term compensatory growth, by contrast, usually denotes the developmental compensation of one perturbed trait by another trait, e.g., the hypertrophy of organs or parts of organs after other parts have been removed (some authors, however, use the terms catch-up growth and compensatory growth synonymously).In the same sense, Waddington envisioned development as a "homeorhetic" process: perturbed by environmental or internal fluctuations, individual development varies around a steady trajectory or target trajectory-which he termed the chreod (a combination of the Greek words for "determined" or "necessary" and for "pathway").When the actual trajectory deviates from its target, it tends to return to the chreod due to canalizing processes, as illustrated by his influential metaphor of the epigenetic landscape (Fig 1; see [26] for a review of this concept and it's impact on modern life science).Individual development is represented by a rolling marble in this landscape.When it deviates from the valley floor (the chreod or target trajectory), e.g., due to environmental influences, the marble will eventually roll down the wall and return to the floor, until it is pushed away again.The steepness of the wall represents the strength of canalization, the tendency to return to the chreod.
Waddington and many later researchers emphasized the role of canalization in evolution.Canalization can facilitate the accumulation of hidden genetic variation that is "blind" to selection; a breakdown of developmental canalization (or gene-gene and gene-environment interactions of canalized alleles) can release this genetic variation and give rise to episodes of rapid evolution and divergence [11-13, 15, 28, 29].Robustness and canalization may also be subject to evolution themselves [30][31][32].Although it remains unclear how selection can enhance robustness, mechanisms that canalize development against environmental and genetic perturbations are considered a prerequisite for evolvable complex life [33][34][35].The nature of these mechanisms, however, remains poorly understood.Basically, two lines of explanations have been proposed [19,35].One focuses on the role of particular genes, such as Hsp90 and other chaperones, in buffering the effects of internal and external perturbations of development [10,[36][37][38].Quantitative trait loci (QTLs) that affect the variance of a trait but not the trait mean are interpreted as evidence for genetic loci that specifically stabilize the phenotype [39,40].However, alleles that destabilize development do not necessarily indicate an evolved basis of genetic canalization [13,19].The second line of research explains canalization as an inherent property of the developmental-genetic architecture, such as where the branching points relate to alternative differentiations.The steepness of the walls represent the strength of canalization, the tendency of individual development to return to the chreod (modified from [27], after [2]).https://doi.org/10.1371/journal.pcbi.1008381.g001nonlinear developmental processes, redundancies among pathways, epistatic interactions, averaging effects of independent sources of variation, and epigenetic effects of tissue interactions [19,[41][42][43][44][45][46][47][48][49][50][51].
Most empirical studies on canalization have compared the variance of phenotypic traits across different genotypes or environments.But reduced cross-sectional variance (i.e., variance across individuals of the same age or developmental stage) can result both from increased canalization of developmental perturbations (stronger or faster reversion to the target trajectory) as well as from reduced effects of genetic or environmental variation on development in the first place (i.e., phenotypic development was never perturbed by these factors).Also the variance of newly arising perturbations of development can change with age or differ between groups.As the target trajectory is usually unobservable and likely to vary across individuals, these phenomena cannot be disentangled by simple comparisons of phenotypic variance.It is even unclear if and how a target trajectory (Waddington's chreod) is directly encoded in the organism's genetic-developmental system or, instead, emerges as an epiphenomenon (e.g., [18]).Variation in developmental timing and episodic growth can also lead to statistical signals similar to that of developmental canalization (see below).
Here, we formalize Waddington's notion of developmental canalization or targeted growth by an extended autoregressive model.We show that under certain assumptions the strength of canalization and the amount of canalized variance in a population can be estimated, or at least approximated, from longitudinal data (measurements of the same specimens at different age stages), even if the target trajectories are unobserved.We extend this model to multivariate measures and discuss reifications of the ensuing parameter matrix.We apply the univariate and multivariate approaches to longitudinal data on human cranial size and shape.

Modeling canalized development
Waddington's view of developmental canalization, as represented by the metaphor of the epigenetic landscape, is closely matched by an Ornstein-Uhlenbeck (OU) process, a mean-reverting continuous-time stochastic process.The OU process, x t , is defined by the following stochastic differential equation: where θ > 0 and σ > 0. As a model of individual development, x t represents a quantitative phenotypic trait at time t, dx t its developmental change, and θ the "strength" by which x is attracted by μ, which corresponds to the target trajectory or the valley floor of the epigenetic landscape (Waddington's chreod).W t represents random and independent perturbations of development.The larger the deviation from the chreod, the faster the reversion (because the reversion is a constant fraction of the perturbation).The slope of the epigenetic landscape thus increases with the distance from the valley floor, giving rise to a concave cross-section of the landscape.
An OU process is a Markov process, that is, the change dx t depends only on its current value x t , not on any earlier values of x.An OU process is also a stationary process, i.e., even though x fluctuates due to random perturbations, in the long run the mean of x approaches the mean reversion level, μ, and the variance of x stabilizes at s 2 2y .Hence, the variance of the developing trait is determined both by the strength of canalization and the variance of developmental perturbations.
Even though organismal development is a continuous process, it is empirically observed at discrete time points.The discrete-time analog of an OU process is an autoregressive model of order 1, AR(1): As for the OU process, the parameters μ and θ correspond to the target trajectory and the strength of canalization, respectively.Independent random perturbations arising at time t are represented by � t .For 0 < θ < 2, x t+1 approaches μ and reduces the deviation from the target trajectory.Complete canalization occurs for θ = 1, whereas 1 < θ < 2 corresponds to an overshooting canalization, leading to a damped oscillation of the developmental trajectory around the chreod.For θ = 0, the absence of any canalization, developmental perturbations accumulate and x shows a random walk.
In the above AR model, the mean reversion level μ is constant.For a realistic account of organismal development, however, the target trajectory must change over time and is also likely to vary across individuals unless they are genetically identical (isogenic) and develop in a homogenous environment.Hence, μ t is not a constant but a variable.Furthermore, empirical studies showed that cross-sectional phenotypic variance can change considerably during mammalian development [18,21,22,52,53]; the strength of canalization may thus differ across developmental stages.In principle, θ t may also differ among individuals, but if canalization is indeed an inherent property of the genetic-developmental system it is likely to be of similar magnitude in a population.For this reason, and for the sake of a tractable model, we will consider θ t constant across individuals but free to change over time, giving rise to an extended model: where c t = μ t+1 − μ t are the developmental changes of the target trajectories (the "deterministic" component of development).While the developmental perturbations (the "stochastic" component), � t , may be assumed to be uncorrelated across time points (see below for a discussion of this assumption), the developmental changes of the target trajectories can be correlated across time, e.g., when the trajectories continually diverge over a longer period.In this case, the deterministic part can itself be decomposed into a function of its lagged values plus new, independent variation: where g t are deterministic effects on development that newly arise at time t, e.g., the effects of developmental genes that are first expressed at this developmental stage.Except for variation in developmental timing (see below), the parameter z t is positive if allelic variation in developmental genes or regulatory factors induces a divergence of growth trajectories.Together, Eqs 3 and 4 represent a first-order autoregressive moving average process.In this model, the magnitude of the mean reversion is the same linear function for positive and negative deviations from the chreod, corresponding to a symmetric cross-section of the epigenetic landscape.In sufficiently large samples, the reversion may be estimated as a more complex, non-linear function of the perturbations.Here we will continue modelling the developmental reversion as the fraction θ t of the perturbation, which can be considered a linear approximation to a nonlinear first-order stochastic process [54].
In the AR(1), canalization is triggered just by the current deviations from the chreod, not by earlier ones.But some developmental perturbations may cause a delayed response, leading to a mean reversion with a longer time lag.This can be represented by an autogressive model of higher order, AR(p): for p � t.In words, the current mean reversion is a function θ t of the deviation at time t plus a function θ t−1 of the deviation at time t − 1, etc.If not known a priori, the maximal time lag can be estimated from the data, e.g., by the partial autocorrelation function, the partial correlations of x t with its own lagged values [55,56].

Developmental plasticity versus elasticity
The developmental target trajectory may be considered genetically determined, but many allele effects and developmental processes are context-dependent.Non-reversible effects of environmental factors on development, usually referred to as developmental plasticity, would then count as a deterministic component of development.By contrast, stochastic short-term perturbations that are subsequently canalized may be considered a developmental elasticity (see also [57, section 2.4]).But the same developmental mechanisms that canalize stochastic short-term perturbations (elasticity) may also influence the long-term effect of an environmental factor (plasticity).Similarly, a developmental mechanism that canalizes environmental perturbations may also buffer genetic perturbations of development (e.g., [58][59][60]) so that, again, the same developmental mechanisms can contribute to both the stochastic and deterministic components of development.Hence, even though stochastic and deterministic aspects of development are distinguished in Eq 3, this distinction is neither conceptually nor biologically obvious.As shown below, also the separate statistical estimation of these components can be difficult.

Estimating the strength of canalization
In Eq 3, the developmental change, x t+1 − x t , is a linear function, θ t , of the perturbation from the chreod plus newly arising variation: the deterministic changes, c t , and the stochastic perturbations, � t .Write d t = x t − μ t for the developmental perturbations at time t, and assume c t , � t , and d t to be mutually uncorrelated as well as � t and d t to be normally distributed with a mean of zero.Then the maximum-likelihood estimate of θ t is given by the least-squares regression of x t+1 − x t on d t (e.g., [56]).In the limit of large sample size, In practice, however, μ t is unobservable and one can regress the developmental change only on x t = μ t + d t instead of d t .Only if μ t is constant among individuals (as, perhaps, in an isogenic sample), the ensuing regression slope is an unbiased estimate of θ t .If μ t varies across individuals, the regression of x t+1 − x t on x t yields a biased estimate of θ t , the empirical estimate ỹt : assuming that cov(μ t , d t ) = 0.This shows that the variance of μ t biases ỹt towards zero, whereas long-term trends of the target trajectories (the effect z t μ t on c t ) negatively bias ỹt .Hence, with increasing (epi)genetic heterogeneity among individuals, the regression of developmental changes on the past trait values increasingly underestimates the strength of canalization.
What is the amount of cross-sectional variance reduced by canalization?By combining Eqs 3 and 4, the trait values at time t + 1 can be written as Under the abovementioned independence assumptions, the change in variance from time t to t + 1 can be decomposed into The variance of x increases by the variance of new stochastic perturbations, var(� t ), and by the variance of the target trajectories, which can be decomposed in newly occurring "deterministic" variation, var(g t ), and the continuation of earlier divergencies of the target trajectories, The variance reduction is largest for θ t = 1, when all deviations from the target trajectories are completely reverted (S1 Fig) .Note that even in the presence of canalization (i.e., 0 < θ t < 2), the cross-sectional variance of x can increase if the newly introduced variance exceeds the canalized variance.
Under the same independence assumptions, and because ðy 2 t À 2y t Þ � À 1:5 y for −0.1 < θ < 0.6, In words, for a small to intermediate canalization coefficient, the canalized variance can be approximated by 1.5 times the covariance between x t and the developmental change x t+1 − x t , as long as newly arising variation among the target trajectories is only weakly associated with the current differences among target trajectories.
In contrast to the estimated strength of canalization, ỹt , the estimation of the canalized variance is not biased by the variance of μ t .It is also unbiased if the target trajectories all change in the same way (parallel trajectories) or if new differences in target trajectories are unrelated to current differences.However, any long-term trend among target trajectories leads to an underestimation of the canalized variance.For example, if developmental trajectories of body size diverge over a longer time period, large individuals grow more than smaller individuals.This appears like an "inverse" canalization that amplifies existing individual differences, and ỹt as well as the estimated canalized variance are biased.

Measurement error
In practice, x is measured with error.It is well known that measurement error of the independent variable biases the least squares regression slope towards zero (e.g., [61]), which also affects the estimation of θ t .But independent error for consecutive measurements of a times series also induces a "regression to the mean" [62,63].Write for the measurement of x t with error η t , and ỹ� t for the empirically estimated θ t based on x � t .Assuming that measurement error is uncorrelated between time points and with both x t and x t+1 , the estimate equals because, intuitively, measurement error at time t "disappears" at time t + 1.Independent error thus leads to the same statistical signal as canalization.If estimated from repeated measurements, the variance of measurement error can be added to covðx � tþ1 À x � t ; x � t Þ to achieve a more realistic estimate of the canalized variance.Because Eqs 3, 4 and 7 together constitute a state space model, the parameters can also be estimated by a series of iterations of regression and prediction steps, known as Kalman filter, as long as the target trajectories are relatively stable and close to linear in the observed time interval and the variance of measurement error is known [56,64,65].

Variation in developmental timing
Variation in the onset of rapid growth episodes leads to an increase and subsequent decrease of cross-sectional phenotypic variation (S2 Fig) .For instance, variation in the onset of the pubertal growth spurt leads to an increased interindividual variance of body height in early adolescence, followed by a decrease of variance after the cessation of growth spurts (e.g., [18]).A very early or late pubertal growth spurt thereby resembles targeted growth.However, such a reduction of variance does not necessarily reflect mechanisms of targeted growth, just variation in developmental timing, which can have a genetic basis or be environmentally triggered.The resulting divergence and subsequent convergence of trajectories can thus be part of the deterministic or the stochastic components of development.However, trait variation due to perturbed developmental timing could also be subject to developmental canalization.In order to avoid a potentially spurious signal of canalization, individual trajectories can be sampled, if possible, at homologous developmental stages rather than at the same ages.

Application to cranial size
We studied developmental canalization in human cranial size using a longitudinal sample of 26 untreated individuals (13 girls, 13 boys) from the Denver Growth Study, a longitudinal Xray study carried out in the US between 1931 and 1966.On a total of 500 lateral radiographs, covering the age range from birth to early adulthood, 18 landmarks were digitized by Ekaterina Stansfield (Fig 2 ; for more details see [66]).Six semilandmarks were located on the external outline of the frontal bone, which were allowed to slide along the bone outline so as to minimize the bending energy of the landmark configurations around their Procrustes average [67,68].Further information on the Denver Growth Study and the original X-rays are available at https://www.aaoflegacycollection.org/aaof_collection.html?id=UOKDenver.As these are historic, anonymized, and publicly available data, no approval of an ethics committee was necessary.
We used the centroid size of these landmarks as a measure of cranial size [69,70].As the centroid size is the square root of the summed squared distances between the landmarks and their centroid, it is measured in units of pixel length, which is 0.1625 mm.Additionally, the area of the lateral projection of the frontal sinus (hereafter, frontal sinus size) was measured on every radiograph.Because the radiographs were not taken at exactly the same age for all individuals, we estimated individual cranial size and frontal sinus size at yearly intervals (2-17 years of age) by unweighted local linear regressions within moving age windows of 2.5 years for each individual.Changing the width of the age window or using a moving average algorithm instead of local linear regression affected the smoothness of the trajectories but lead to qualitatively similar results.
We estimated measurement error by measuring eight age stages for four individuals four times each.Because centroid size is a composite variable of all the landmark coordinates, its measurement error tends to be low.The variance of repeated measurements was in the order of 1-2% of between-individual variation for the different ages.The magnitude of measurement error did not correlate with age.
Average cranial size increased fastest until approximately five years of age.Sexual dimorphism in cranial size was relatively constant from 2 to 12 years of age and increased with the beginning of puberty (Fig 3A).The partial autocorrelation function indicated that only a time lag of 1 is relevant (Fig 3B ), i.e., if conditioning on x t , the earlier observations, x t−1 , x t−2 , . .., did not show a relevant correlation with x t+ 1 .We thus applied an autoregressive model of order 1 to these data.Iterative autoregressive-moving-average (ARMA) estimation and Kalman filter did not perform well, presumably because of relatively unstable, nonlinear growth trajectories.
The cross-sectional variance in cranial size was approximately halved from 2 to 5 years of age and started to increase again at about 10 years.The magnitude of canalized variance-estimated as 1:5 covðx � tþ1 À x � t ; x � t Þ-was high in the first few years and decreased to about zero at six years of age (Fig 4A).From 13 to 16 years, the "canalized" variance was slightly positive, reflecting the continual divergence of male and female growth trajectories.After removing sexual dimorphism from the data by subtracting the age-and sex-specific average from all individuals, this effect went away (Fig 4B).The weak signal of canalization at age 16 coincides with the pubertal growth spurt.Because measurement error was negligible (pooled over all age stages, the variance of repeated measurements was only 81.8), we did not correct for it.
The estimated strength of canalization, ỹ� t , was in the order of 0.5 at age 2-3 and approached zero at about 6 years (Fig 5).The estimates ỹ� t are biased by interindividual variation in the target trajectories.However, correcting for sexual dimorphism, which removes a considerable fraction of variance in target trajectories, only weakly affected the estimated strength of canalization (solid versus dashed lines in Fig 5).
We performed the same analyses also for frontal sinus size.Average fontal sinus size as well as the variance in sinus size increased throughout the studied age range (Fig 6).In contrast to cranial size, frontal sinus size did not show any signs of developmental canalization.The estimated canalized variance was even slightly positive, indicating that individuals have a tendency to develop small or large frontal sinuses early on.Removing sexual dimorphism slightly lowered all three variance components.

A multivariate model of canalization
The univariate model in Eq 3 can be extended to multivariate traits:  arising perturbations at time t, respectively.Note that the effect of x t on developmental change is represented by the p × p matrix Θ t , instead of a scalar as in the univariate case.The diagonal elements, Θ t,ii , quantify how the i-th variable at time t affects its developmental change from t to t + 1, conditioned on the other p − 1 variables.The off-diagonal elements, Θ t,ij , quantify how the i-th variable at time t affects the developmental change of the j-th variable.In this sense, they reflect the strength by which one variable compensates the perturbation of another variable.These elements don't have a correspondence in Waddington's epigenetic landscape.
In the absence of measurement error and variation between target trajectories and if ϵ t is multivariate normal with the zero vector as mean, the multivariate multiple regression of x t+1 − x t on x t yields an unbiased estimate of Θ t .In the limit of large sample size, assuming that c t , d t , and ϵ t+1 are mutually uncorrelated (which is increasingly unrealistic as p increases) and that the variance-covariance matrix cov(d t ) is invertible.As for the univariate case, variation between target trajectories and measurement error bias the empirical estimate Θ� t , and the regression of x � tþ1 À x � t on x � t generally underestimates the strength of canalization.Extending Eq 6 to the multivariate model shows that the developmental change in variance-covariance structure due to canalization is A common multivariate measure of total variance is the sum of the variances of all variables, Tr(cov(x t )).The developmental change in total cross-sectional variance due to canalization is If the elements of Θ t are small, the quadratic form Θ 0 t covðd t ÞΘ t is negligible, and the developmental change in variance-covariance structure is mainly determined by cov(d t )Θ t and its transpose.Because the effects of canalization and compensation on the cross-sectional variance-covariance structure can be approximated by the cross-covariance matrix cov(x t , x t+1 − x t ) if the strengths of canalization and compensation are small and the changes of the target trajectories only weakly related to their current states, and if the compensatory effects are approximately symmetric (i.e., Θ t,ij � Θ t,ji ).
In many morphometric contexts (e.g., for landmark shape coordinates, collections of interlandmark distances, or voxel gray values), the measured variables cannot be separately biologically interpreted, and biometric analyses are usually based on linear combinations of the variables.For such data, the elements of Θ t bear no biological meaning by themselves and are typically too numerous for inspection.Hence, Θ t may be expressed in its canonical form: where U and V are p × p matrices of left and right singular vectors, and Λ is a diagonal matrix of singular values.The singular vectors can be interpreted as multivariate traits (axes in phenotype space) with maximal canalization (these traits, however, may not necessarily be interpretable separately; e.g.[71,72]).For instance, the regression of the linear combination (x t+1 − x t ) 0 v 1 on x 0 t u 1 has the slope Λ 1,1 , where u 1 and v 1 is the first pair of singular vectors of Θ t , and Λ 1,1 is the largest singular value.This decomposition of Θ t allows one to reduce the measured variables to a small set of linear combinations, which can be interpreted as latent variables underlying the observed patterns of canalization and compensation.In the multivariate statistical literature, this is referred to as reduced rank regression [73,74].
As the change in variance-covariance structure attributable to canalization is approximately proportional to the cross-covariance matrix between trait values at time t and their developmental changes from t to t + 1, one can also compute linear combinations that maximize canalized variance (instead of the strength of canalization) by the singular value decomposition of the cross-covariance matrix: Here, the linear combinations (x t+1 − x t ) 0 v 1 and x 0 t u 1 have the maximal covariance Λ 1,1 (as a proxy of canalized variance) among all unit vectors u and v.The second pair of linear combinations, (x t+1 − x t ) 0 v 2 and x 0 t u 2 , is orthogonal to the first pair and has the second highest covariance, and similarly for subsequent pairs.This approach is equivalent to two-block partial least squares analysis (PLS), which is commonly used in morphometrics (e.g., [75][76][77]).As the singular values are always positive, interpretations in terms of canalization and compensation must be based on the similarity of the singular vectors, u i , v i , as illustrated by the example below.E.g., pure canalization is indicated by u 0 i v i � À 1, whereas compensation is indicated by u 0 i v i � 0. Note that a stable estimation of Θt requires a sample size that greatly exceeds the number of variables (n � p) because of the matrix inversion of cov(x t ).In practical applications, the SVD of the cross-covariance matrix cov(x t , x t+1 − x t ) is more reliable, and the associated singular values less affected by sample size than the singular values of Θt .

Application to craniofacial shape
The 500 configurations of 18 landmarks (Fig 2 ) were superimposed by generalized Procrustes analysis in order to standardize for variation in position, orientation, and scale.The resulting 36 shape coordinates were used to quantify craniofacial shape [69,70,78].Like for centroid size, we estimated craniofacial shape at yearly intervals by unweighted local linear regressions within a moving age window of 2.5 years for each shape coordinate and each individual separately.Sexual dimorphism was removed for each age group.We computed the cross-covariance matrix between each age group and the subsequent developmental change over one year, covðx � t ; x � tþ1 À x � t Þ, in order to explore the effect of canalization on the cross-sectional variancecovariance pattern.Using different age windows, longer time lags, or the first few principal components of the shape coordinates instead of all 36 coordinates led to very similar results as those presented here.
When estimated by 1:5 Trðcovðx � t ; x � tþ1 À x � t ÞÞ, canalization reduced total shape variance until about 6-7 years of age.Cross-sectional shape variance, however, decreased only slightly in this age period, indicating that newly added variance and canalized variance were similar in magnitude ( Fig 7).During puberty, craniofacial shape variance increased rapidly and decreased thereafter, reflecting variation in the timing of pubertal growth episodes.
However, pooling the variances of all shape coordinates conceals the heterogeneity of the different aspects on cranial shape [53,79].Fig 8A shows the first dimension of the PLS analysis between craniofacial shape at age 2 and the shape changes from 2 to 3 years of age (i.e., the first two pairs of singular vectors of covðx � 2 ; x � 3 À x � 2 Þ), accounting for 50.3% of total squared covariance between craniofacial shape and the subsequent shape change.The first left singular vector, u 1 , mainly represented the flexion of the cranial base and the relative size of the maxilla at age 2, which together determine the size of the pharynx.The right singular vector, v 1 , was very similar to the left singular vector but in opposite direction, which indicates a pattern of canalization: individuals with a relatively small pharynx at age 2 experience an increased flexion of the cranial base and reduced maxillary growth from age 2 to 3, thus leading to an increase of pharyngeal space.By contrast, individuals with a relatively large pharynx at age 2 show decreased pharyngeal growth.The individual PLS scores, i.e., the linear combinations The second PLS dimension, accounting for 14.9% of total squared covariance, mainly reflected the canalization of relative facial height.Individuals with a relatively high face at age 2 tend to grow less in facial height than individuals with a lower face, and vice versa (Fig 8B).The correlation between the corresponding scores was 0.70 and the regression slope ỹ� ¼ 0:35.Variance in facial height was canalized during the first 7 years of age and was amplified again during puberty.The third PLS dimension reflected the angulation and relative size of the posterior cranial base, including its effect on relative pharyngeal size (not shown), which was subject to canalization until 6 years of age (Fig 7C and 7D).Dimension 4, by contrast, which corresponded to the relative orientation of the face and the pharynx (not its size), did not show any signs of canalization; it even increased considerably in variance during puberty (Figs 9A, 7C and 7D).For the subsequent dimensions, the left and right singular vectors differed in direction and, hence, did not clearly reflect patterns of canalization; also their crosssectional variances did not decrease.
From about 5 to 10 years, the corresponding PLS analyses revealed another pattern of canalization, mainly affecting the thickness of the cranial base and, indirectly, also the dimensions of the upper pharynx (Fig 9B).Several craniofacial shape features showed strong canalization from 15 to 16 years of age, which likely owes to variation in the timing of pubertal growth periods (see also Fig 7).

Discussion
Target trajectories of individual growth usually are unobservable.But we showed that under certain assumptions the strength of developmental canalization as well as the canalized analysis between the shape variables at age 2 (first block) and the shape changes from age 2 to 3 years (second block).(A) The first PLS dimension represents the shape pattern with the largest canalized variance.Individuals with a relatively large pharynx due to a retroflexed cranial base and a relatively short maxilla at age 2 (low scores along u 1 ) tend to experience increased maxillary growth and basicranial flexion from age 2 to 3 years (low scores along v 1 ).By contrast, individuals with a highly flexed cranial base and relatively long maxilla, both of which constrain the pharynx, tend to undergo a more than average increase in pharyngeal size (high scores along u 1 and v 1 ).Note that the left deformation grids represent deviations from the mean shape at age 2, whereas the right grids represent deviations from the mean shape change form 2 to 3 years.Hence, the maxilla does not shrink in individuals with a large maxilla, but it grows less than average.(B) The second PLS dimension represents canalization of facial height: individuals with a relatively high face (inferiorly placed maxilla) undergo reduced facial growth, whereas individuals with a relatively low face (including a small pharynx and nasal cavity) show a more than average increase in facial height.Together, this indicates strong canalization of the relative size of the upper airways in early postnatal development.
https://doi.org/10.1371/journal.pcbi.1008381.g008variance can still be approximated by the regression slope and the covariance, respectively, between a trait value and its subsequent developmental change.The strength of canalization corresponds to the slope in Waddington's epigenetic landscape and reflects the ability or propensity of perturbed individual growth to return to the individual's target trajectory.The canalized variance, by contrast, is the actual average squared reduction of developmental perturbations in the observed sample.Hence, a strongly canalized trait can show little canalized variance if its growth is widely unperturbed.
These estimations rest on the assumption that the developmental perturbations (d t , � t ) are statistically independent of the target trajectories (μ t , c t ).This assumption seems plausible unless individuals with particular trait values or alleles are considerably more susceptible to environmental or genetic perturbations than other individuals.A second crucial assumption is that the developmental perturbations are uncorrelated across time points.For observations close in time this assumption does not necessarily hold because environmental influences or the expression of particular genes may persist.Positive correlations between perturbations would lead to an underestimation of canalization, but this bias can be reduced by keeping the time lag between observations sufficiently long.The estimated strength of canalization is biased by inter-individual variation across target trajectories (Eq 5).For isogenic samples, e.g., inbred lab strains, this bias may be small because the individuals can be assumed to have very similar target trajectories.In natural populations, including humans, the strength of canalization is underestimated.However, in our analysis we found that correcting for sexual dimorphism, which considerable reduces variation among target trajectories, only weakly affected the estimates of canalization.Hence, strong signals of canalization are likely to be identifiable also in natural populations.
The estimation of the canalized variance is more robust than that of the strength of canalization.It is not biased by variation among target trajectories, only by the covariance between newly arising differences among target trajectories and their previous states (e.g., long-term divergences of trajectories).In isogenic samples this covariance is close to zero, and in natural populations it decays with the length of the time lag between observations.In summary, even if the target trajectories cannot be observed or modelled, the strength of canalization can nonetheless be estimated if genetic and environmental heterogeneity between individuals is small.The canalized variance can even be estimated in natural populations if the time lags between the observations are sufficiently long.As measurement error biases both estimates, it is important to keep the error small and constant across time points and groups, or to estimate it by repeated measurements.
The study of developmental canalization for multivariate traits goes beyond the epigenetic landscape metaphor.For a single variable the strength of canalization is represented by the parameter θ t , whereas its representation for p variables requires a p × p matrix Θ t , comprising p canalization parameters and p(p − 1) compensation parameters.For a small number of biologically meaningful measurements, the estimates of these parameters can be directly interpreted.But for larger numbers of variables, or for measurements that cannot be separately biologically interpreted (shape coordinates, voxel gray values), the matrix Θ t must be decomposed into interpretable linear combinations.Multivariate studies of canalization, therefore, may mostly focus on the exploratory identification of features with strong or weak canalization rather than on a precise estimate of the strength of canalization.Linear combinations with maximal strength of canalization can be computed by a reduced-rank regression between the trait values and their subsequent change.Computationally more stable, however, is partial least squares analysis, which yields linear combinations with maximal covariance between the trait values and their subsequent change, which is approximately proportional to the canalized variance.Isotropic measurement error biases the covariances but has little effect on the extracted linear combinations.

Developmental canalization of human craniofacial form
Despite a relatively small sample, we found clear signals of developmental canalization that varied both across cranial regions and age stages.Craniofacial size was strongly canalized during the first five years of postnatal development.Canalization decayed thereafter and showed a small peak again during puberty, which presumably owes to variation in the timing of pubertal growth spurts.As a result, cross-sectional variance in craniofacial size almost halved between 2 and 5 years of age and continually increased thereafter.Decreasing phenotypic variance of cranial size and body size during postnatal development as well as catch-up growth after periods of retardation have also been reported for many other vertebrates and even insects [18,22,24,[80][81][82].Quantitative genetic studies found that most of the reduced variance is attributable to environmental and maternal effects, whereas additive genetic variance reduces less or even increases during postnatal development [21,22,52,[83][84][85], suggesting that variation in target trajectories has a strong genetic component.
We contrasted the developmental dynamics of craniofacial size with that of frontal sinus size (as inferred from the lateral X-rays), which is known to be highly variable in adults [86,87].For the oldest age group in our data, the coefficient of variation (CV; standard deviation/ mean) was 0.51 for frontal sinus size but only 0.03 for craniofacial size (S3 Fig) .The variance of frontal sinus size increased continually in the observed age range, but the CV decreased from 4.0 at age 2 to 1.0 at age 6.However, this rapid decrease of the CV does not imply developmental canalization: we found no signs of canalization for frontal sinus size.The estimated strengths of canalization were even slightly negative, indicating that early individual differences were amplified (rather than canalized) during postnatal development.Presumably, the sinuses lack any canalization because they "grow" only by osteoclast activity and their exact size is of little functional relevance [88].The CV for sinus size just decreased because the mean size increased faster during development than the standard deviation [51,89].Several earlier studies on targeted growth relied on the comparison of CVs across age groups, but these results should be interpreted with care.
Total variance of craniofacial shape slightly decreased from 3 to 5 years of age and then increased by about 30% until age 17.Postnatal reduction of cranial shape variance was also found in other mammals [18,[90][91][92][93].However, such an overall assessment of shape variance partly depends on the multivariate measure of variance and conceals the fact that different shape features (linear combinations of shape variables) can deviate considerably in their variational dynamics [53,79,94].We used partial least squares (PLS) analyses between the shape variables and their subsequent developmental change over one year to identify canalized features of craniofacial shape.We found multiple aspects of cranial shape that show strong canalization during postnatal development, particularly the flexion and thickness of the cranial base as well as the relative size and position of the upper jaw, all of which determine the dimensions of the nasopharynx.Other features, such as the orientation of the entire face and posterior cranial base, showed no canalization and continually increased in cross-sectional variance.
In the first few years of human postnatal development, the cranial base flexes and the upper jaw increases in length, both of which constrain the pharynx [66,95,96].Accordingly, obstruction of the upper airways is relatively common in early childhood [97].The precise and coordinated growth of the cranial components forming the pharynx thus is crucial to respiratory function.However, in contrast to the pharynx, most other functional components of the cranium, such as the braincase, the orbits, and the jaws, develop in tight interaction with the encapsulated tissues (brain, eyes, tongue and teeth), which tightly coordinates bone growth [95,[98][99][100].This lack of a "functional matrix" in the pharynx may make pharyngeal growth particularly prone to perturbations despite its functional importance.Indeed, we found that relative pharyngeal size (PLS dimension 1) had by far the highest variance at age 2 among all the extracted shape features, but this variance was continually canalized throughout postnatal development (Fig 7C).By contrast, facial orientation (PLS dimension 4), which is of little functional importance, continually increased in variance until it became the most variable component of adult facial shape.
Signals of developmental canalization can also result from variation in developmental timing and from episodic growth of the entire structure or of its components (e.g., [18], S2 Fig) .But our findings of canalization cannot be completely explained by these phenomena.The canalization of cranial size during early postnatal development was much stronger than that observed during and after puberty, and removing average sexual dimorphism in cranial size had little effect on these signals (Figs 4 and 5).Likewise, the shape features related to relative pharyngeal size showed considerably more canalization in early postnatal life than during puberty.Furthermore, developmental dynamics differed strongly between the different shape features, and the age periods of strong canalization did not coincide with the periods of fastest developmental change.Increasing the time lag for the estimations of θ t and the canalized variance smoothed the curves but did not lead to qualitatively different results.
In summary, craniofacial size as well as multiple aspects of craniofacial shape show a behaviour that resembles developmental canalization in the sense of an auto-regulated or homeorhetic process.Whether or not this behaviour involves a predefined, genetically encoded target trajectory (Waddinton's chreod) is still unclear.At least for the cranium, it seems most plausible that the targeted behaviour of development emerges from the epigenetic interaction of the developing tissues together with controlled gene expression.Developing soft tissue, such as the brain, eyes, and tongue, provide a matrix that affects and integrates the growth of the surrounding cranial elements.Likewise, sutural growth mechanically induced by stunted or overshooting growth of the adjacent bones canalizes overall cranial shape despite ample variation in the relative dimensions of the constituting bones [51].The lack of epigenetic growth induction in the upper airways may account for the hyper-variability of nasopharyngeal dimensions in infancy and early childhood.It remains to be studied which mechanisms (e.g., induction by the air flow) account for the pronounced canalization of the nasopharynx during later postnatal development.Understanding the pattern and timing of developmental canalization in the cranium can also aid the effective timing of orthodontic treatment and the reduction of relapse [101,102].The methods proposed here will allow for the comparison of developmental canalization and targeted growth across populations, genotypes, treatments, and age groups in future research.

Fig 1 .
Fig 1. Epigenetic landscape.Conrad Waddington's epigenetic landscape is a landscape over a phenotype space and a time dimension, illustrating the canalization of individual development (the marble) around a target trajectory (the valley floor).Branching valleys represent cell fates, or more generally, predetermined ontogenetic pathways ("chreod"),where the branching points relate to alternative differentiations.The steepness of the walls represent the strength of canalization, the tendency of individual development to return to the chreod (modified from[27], after[2]).

Fig 2 .
Fig 2. Landmark scheme.(A) A lateral cranial radiograph of an adult human man with 13 anatomical landmarks: basion (Ba), bregma (Br), clival point (CP), foramen ceacum (FC), nasion (Na), nasospinale (Nsp), greater wings of the sphenoid (Pmp), prosthion (Pr), posterior nasal spine (Pns), planum sphenoideum (PS), sella (S), sphenoidale (Sph), and glabella (Gl).Glabella as well as five points along the frontal bone were treated as semilandmarks; their exact locations were estimated by the sliding landmark algorithm.(B) Mean shape configuration of these 18 landmarks in the sample of 26 individuals, measured from birth to adulthood.The midsagittal outline of the frontal bone is represented by the black curve, and the shape of the cranial base and the maxilla by the two polygons.https://doi.org/10.1371/journal.pcbi.1008381.g002

Fig 7 .
Fig 7. Canalization of craniofacial shape.(A) Total craniofacial shape variance from 2 to 17 years of age.(B) Total canalized shape variance, approximated by 1:5 Trðcovðx � t ; x � tþ1 À x � t ÞÞ: (C) Variance of the first four dimensions of the two-block partial least squares (PLS) analysis between craniofacial shape at age 2 and the shape changes from age 2 to 3 (see Fig 8).In other words, these are the cross-sectional variances of the shape features with maximal canalized variance from age 2 to 3. (D) Approximated canalized variance of these four PLS dimensions.The first dimension (blue curve) shows canalization until 6 years and again from 9 to 13, leading to a reduction of cross-sectional variance throughout postnatal development.Dimensions 2 and 3 only show canalization until 6 years of age, whereas dimension 4 shows no signs of canalization (individual differences were even amplified rather than reduced).https://doi.org/10.1371/journal.pcbi.1008381.g007

Fig 8 .
Fig 8. Canalized shape features.Canalization of craniofacial shape from 2 to 3 years of age was explored by two-block partial least squares (PLS) analysis between the shape variables at age 2 (first block) and the shape changes from age 2 to 3 years (second block).(A) The first PLS dimension represents the shape pattern with the largest canalized variance.Individuals with a relatively large pharynx due to a retroflexed cranial base and a relatively short maxilla at age 2 (low scores along u 1 ) tend to experience increased maxillary growth and basicranial flexion from age 2 to 3 years (low scores along v 1 ).By contrast, individuals with a highly flexed cranial base and relatively long maxilla, both of which constrain the pharynx, tend to undergo a more than average increase in pharyngeal size (high scores along u 1 and v 1 ).Note that the left deformation grids represent deviations from the mean shape at age 2, whereas the right grids represent deviations from the mean shape change form 2 to 3 years.Hence, the maxilla does not shrink in individuals with a large maxilla, but it grows less than average.(B) The second PLS dimension represents canalization of facial height: individuals with a relatively high face (inferiorly placed maxilla) undergo reduced facial growth, whereas individuals with a relatively low face (including a small pharynx and nasal cavity) show a more than average increase in facial height.Together, this indicates strong canalization of the relative size of the upper airways in early postnatal development.

Fig 9 .
Fig 9. Further features of craniofacial shape.(A) The fourth dimension of the PLS between craniofacial shape at age 2 and the shape change from age 2 to 3 years.This dimension reflects the relative orientation of the face and posterior cranial base, but unlike the previous dimensions these shape features are not canalized; existing individual differences are even amplified.(B) Canalization of craniofacial shape from 6 to 7 years of age (first PLS dimension).Individuals with a relatively thick cranial base at age 6 undergo reduced thickening of the cranial base, whereas individuals with a thin cranial based show increased growth in this area.Again, this indicates canalization of the relative dimensions of the upper airways.https://doi.org/10.1371/journal.pcbi.1008381.g009 , for cranial size before and after removing sexual dimorphism (solid and dashed lines, respectively).
ϵ t ; ð8Þ where x t , μ t , c t , d t , and ϵ t are p-dimensional random vectors of trait values, chreod values, changes of the chreod values, developmental perturbations from the chreods, and newly Fig 5. Strength of canalization in craniofacial size.The estimated strength of canalization, ỹ� t