Bone Fusion in Normal and Pathological Development is Constrained by the Network Architecture of the Human Skull.

Craniosynostosis, the premature fusion of cranial bones, affects the correct development of the skull producing morphological malformations in newborns. To assess the susceptibility of each craniofacial articulation to close prematurely, we used a network model of the skull to quantify the link reliability (an index based on stochastic block models and Bayesian inference) of each articulation. We show that, of the 93 human skull articulations at birth, the few articulations that are associated with non-syndromic craniosynostosis conditions have statistically significant lower reliability scores than the others. In a similar way, articulations that close during the normal postnatal development of the skull have also lower reliability scores than those articulations that persist through adult life. These results indicate a relationship between the architecture of the skull and the specific articulations that close during normal development as well as in pathological conditions. Our findings suggest that the topological arrangement of skull bones might act as a structural constraint, predisposing some articulations to closure, both in normal and pathological development, also affecting the long-term evolution of the skull.

Introduction

Craniofacial articulations are sites of primary bone growth and remodeling; adequate formation and maintenance of these articulations is therefore important for a healthy development of the head and brain. The timely closure of bone articulations is a normal process that takes place during skull development. Craniosynostosis is a pathological condition with an estimated prevalence of about 5 in 10,000 live births[1], in which one or more articulations between cranial bones (frontal, parietal, temporal, and occipital) close prematurely, leading to the fusion of these bones. This premature fusion of bones, if not treated surgically, can cause head malformations due to compensatory growth of other joints[2], sometimes provoking severe brain damage due to an increase of intracranial pressure[3]. Craniosynostosis can occur in isolation, as non-syndromic craniosynostosis[4, 5], or as part of a variety of congenital disorders, such as Apert and Crouzon syndromes[6].

In general, it is not well understood which factors predispose some articulations but not others to close prematurely. It is known that both genetic and non-genetic factors participate in the formation and maintenance of craniofacial articulations through life. The number of genes identified to be carrying mutations associated with craniosynostosis has grown in the last two decades[7]; for example, more than 60 genes have been shown to carry mutations associated with craniosynostosis[7]: some of them show specificity for a suture in the context of a syndrome (e.g., ASXL1 and metopic suture in the Bohring-Opitz syndrome), others predispose to more than one type of craniosynostosis (e.g., FGFR2 in coronal, sagittal, and multi-suture synostoses), while most of them are not specifically associated with suture development, but to osteogenesis in general (e.g., ALX4, EFNA4, and TGFBR2). Non-genetic factors are even less specific than genetic ones and include, among many others, bio-mechanical stress, hypoxia, and use of drugs or smoking during pregnancy[5, 8–11].

Here, we address the susceptibility of articulations to close from a theoretical standpoint, by modeling the skull as a network in which nodes and links formalize bones and their articulations at birth (Fig. 1). This network model is thus a mathematical representation of the entire pattern of structural relations (i.e., physical contacts or articulations) among skull bones. Anatomical network models have been used before, for example, to identify developmental constraints in skull evolution[12, 13], analyze the evolution of tetrapod disparity in morphospace across phylogeny[14], and model the growth of human skull bones[15]. A recent comparison of network models of craniosynostosis conditions showed that, despite the associated abnormal shape variation, skulls with different types of craniosynostosis share a same general pattern of network modules[16].

Figure 1

The arrangement of bones in the human skull at birth modeled as a network (top). Nodes represent bones and links represent articulations among bones (cartilaginous and fibrous joints). Red links are articulations associated with craniosynostosis conditions; dashed links are articulations that close during the normal development of the skull. Note that the metopic suture between the left and right frontal bones closes in both pathological and normal development. Drawings illustrate the shape of the head in some of the conditions studied (bottom).

We infer the susceptibility of craniofacial articulations to close prematurely using the reliability formalism developed for network models[17]. A common feature of the topology of complex networks such as the skull is that one can identify groups of nodes (bones) that have well-defined patterns of connections (i.e., craniofacial articulations or synarthrosis) with other groups of nodes[17]. Such formalism allows us to identify connections that are not ‘expected’ to occur in the context of the entire topology of the network. Since the network represents the actual anatomy of the skull at the newborn stage, the biological processes behind the ‘topological unexpectedness’ of some articulations can also be interpreted as the result of the developmental processes that shape the anatomy of the skull; for example, position of ossification centers, growth patterns, and/or presence of functional matrices[15, 18]. If the architecture of the skull is driving (or influencing) the closure of articulations, we surmise that there is a relationship between the susceptibility of a pair of bones to fuse and the ‘topological unexpectedness’ of their articulation. To quantify such ‘unexpectedness’, we use the link reliability score, that is the probability that a connection exists in the network given the observed (neonatal) topology of the skull[17]. A low score means that the presence of this articulation is rare, that is, not commonly expected in the given arrangement of bones (see Methods for details on how this is estimated). Importantly, the link reliability formalism has been used in other complex systems to accurately predict missing and spurious interactions in social, neural, and molecular networks[17], to predict harmful interactions between pairs of drugs[19], and to predict the appearance of conflicts in teams[20]. Here we use the reliability formalism to investigate whether the topological arrangement of bones predicts which articulations are more susceptible to close in development; in other words, we want to assess if the architecture of the skull acts as an agent that constrains the fusion of bones.

Methods

A network model of the skull

We built a network model of the human skull at birth based on anatomical descriptions[21] and information of ossification timing and fusion events[22]. The nodes and links of the network model formalize the bones and articulations of the skull, respectively (Fig. 1). For simplicity, we use bone in a broad sense to refer both to bonny elements (e.g., a parietal bone) and well-formed cartilaginous templates of the future bones (e.g., the ethmoidal bone). Likewise, we use the term articulation to refer to the cartilaginous (synchondroses) as well as fibrous joints (sutures) of the skull. We are aware that each type of skeletal element and articulation has different biological properties, which might be hard to compare in some contexts. However, our theoretical analysis focuses on a higher level of abstraction, that of topology (i.e., the arrangement of constituent parts), aiming to extract relevant information from the sole topological structure of the skull. Thus, specific properties of nodes (e.g., cellular origins, ossification mechanisms) and of articulations (e.g., contact areas, tensile properties) have not been included in the present model (see ref. 23 for a review of examples of how anatomical network analysis abstractions have successfully been applied in different anatomical contexts).

Topological Organization of the Neonatal Skull

The topological organization of the skull varies during pre- and postnatal development. We have chosen to work with the skull configuration at birth because it allows a broader comparison between closed and persistent articulations, both in normal and pathological conditions. What follows is a summary of the bones present at birth that we used to build the neonatal skull network model (for details, see refs 20 and 21).

The occipital bone at birth consists of four units: a ventral basilar part, a more dorsal occipital plate, and two lateral parts. Around the fourth year the occipital plate and the lateral parts fuse into one unit. Around the sixth year the basilar part is also fused together. During adulthood (about 18–25 years) the occipital bone and the sphenoid bone fuse into a single unit. The frontal bone at birth consists of two halves separated by the metopic suture. Around the eighth year the metopic suture obliterates and the two halves of the frontal fuse into one single bone (although in some individuals the suture endures and left and right frontals are present through life). The premature fusion of the metopic suture is one of the craniosynostosis conditions included in the present study (see Fig. 1). Each temporal bone at birth consists of two parts: the petromastoid and the squama (to which the tympanic ring has united shortly before birth). Around the first year the petromastoid and squama fuse into a single unit. The temporal bone has a tight relationship with two small structures: the ear ossicles (maellus, incus, and stapes) and the styloid process (tympanohyal part and stylohyal part). The former structures develop partially embedded within the temporal bone, while the latter structures fuse with it during the first years of development. For simplicity, we have decided not to consider these structures as separate nodes in the network model; instead, we include them within the temporal bone in order to focus on the main skeletal units of the skull. The sphenoid bone at birth consists of three parts: a central body (including the small wings) and two lateral parts or alisphenoids (comprising the great wings and the pterygoid processes). Around the first year the sphenoid body and the alisphenoids fuse together. As we already mentioned, the sphenoid and the occipital fuse into a single unit during adulthood. The ethmoid bone is still a cartilaginous template at birth, which will later ossify endochondrally to form the ehtmoid bone. The maxilla and premaxilla (one of each per side) at birth are still separated by a suture that can persist until well into adulthood. Each zygomatic bone consists of one single skeletal structure at birth, although sometimes can be divided horizontally in an upper and a lower part. The vomer at birth consist of two lamellae, which fuse together at puberty (although sometimes there are traces of their paired laminar origin). Finally, the lacrimals, nasals, inferior nasal conchae, palatines, and parietals are well-formed skeletal units at birth (although the parietal and palatines still will continue growing some time after birth). At times, the parietal bone can be divided by a longitudinal suture in an upper and a lower part (as this is a deviation of the more common pattern found in humans, we did not include this phenotype in our network model).

Estimation of Link Type Probability Using Stochastic Block Models

Stochastic block models are good models to describe the patterns of connections in complex networks. In such models, nodes are assigned to groups and the probability of a link existing between two pairs of nodes is given by a matrix that specifies the connectivity rate between nodes belonging to pairs of blocks. For a given network, good stochastic block models are those that group nodes that have a similar pattern of connections; for instance, in our case we could group together nodes vomer and palatine since both tend to connect to similar nodes (sphenoid, ethmoid, maxilla) along with a disconnection to similar nodes (e.g., parietal, zygomatic, frontal). Within this description, links between pairs of nodes that belong to groups that are densely interconnected are more likely than those links between pairs of nodes belonging to groups that are sparsely connected. For instance, in the previous example an articulation existing between nodes palatine and maxilla is much more likely than a suture between nodes palatine and parietal. Biologically, the probability that a pair of bones connect depends on the developmental processes that determine the spatial location and growth patterns (direction and speed) of each ossification center, as well as the presence of functional matrices[15].

To mathematically formalize this intuition, we compute the reliability score, that is the probability that a link exists given the network of connections we observe (the newborn skull in our case) using stochastic block models as the basis for our inference algorithm. In practice, our algorithm samples the space of partitions of nodes into groups taking into account how good a given partition manages to classify nodes with similar patterns of connections into the same group. For each of these partitions, each link between a pair of nodes (i, j) has a specific probability. The reliability score of link N is then a weighted average of the probabilities of that link for each sampled partition. Mathematically, we formalize the previous arguments in a Bayesian framework as follows. Given a family of models ε, the probability that N given the observed network N (that is the matrix of connections) is[17] where the integral is over all the models M in ensemble ε. We can rewrite this equation using Bayes theorem and obtain[17, 24] Here, p(N |M) is the probability of the observed interactions given model M and p(M) is the a priori probability of a model, which we assume to be model-independent p(M) = const. In our approach, we assume that the family of stochastic block models is a good ensemble to describe the connectivity in a complex network (in our case that of the human skull). Therefore, each model M = (P, ) is completely determined by a partition P of bones into groups and the group-to-group interaction probability matrix . For a given partition P, the matrix element Q is the probability of an articulation joining a bone in group α with a bone in group β. Thus, if i belongs to group σ and j to group σ we have that[24] andwhere n is the number of articulations between bones in groups α and β and n is the number of disconnections between bones in groups α and β.

The integral over all models in ε can be separated into a sum over all possible partitions of the bones into groups, and an integral over all possible values of each Q . Using this together with Equation 2 to 4, and under the assumption of no prior knowledge about the models (p(M) = const), we havewhere and the integral is over all values of Q , and Z is the normalizing constant (or partition function). Since the dependence on the Q s factorizes, one can carry out analytically the integral over the Q s. Specifically, using that we can express Eq. 5 as,where the sum is over all partitions of bones into groups, n  = n  + n is the total number of possible sutures between groups σ and σ , and H(P) is a function that depends on the partition only

This sum can be estimated using the Metropolis algorithm[17, 25] as detailed next.

Implementation Details

The sum in Equation 6 cannot be computed exactly because the number of possible partitions is combinatorially large, but can be estimated using the Metropolis algorithm[17, 25]. This amounts to generating a sequence of partitions in the following way. From the current partition P , select a random bone and move it to a random new group giving a new partition P . If H(P ) < H(P ), always accept the move; otherwise, accept the move only with probability . By doing this, one gets a sequence of partitions {P } such that one can approximate the sum in Equation 6 as ref. 25 where S is the number of sampled partitions in {P }.

In practice, it is useful to “thin” the sample {P }, that is, to consider only a small fraction of evenly spaced partitions so as to avoid the computational cost of sampling very similar partitions which provide very little additional information. Moreover, one needs to make sure that sampling starts only when the sampler is “thermalized”, that is, when sampled partitions are drawn from the desired probability distribution (which in our case is given by e −/Z). Our implementation automatically determines a reasonable thinning of the sample, and only starts sampling when certain thermalization conditions are met. Therefore, the whole process is completely unsupervised. The source code of our implementation of the algorithm is publicly available from http://seeslab.info/downloads/network-c-libraries-rgraph/ and http://github.com/seeslab/rgraph.

Statistical Analysis

We performed independent Mann-Whitney U tests for the following comparisons: (1) articulations affected by non-syndromic craniosynostosis vs. articulations unaffected; and (2) articulations normally closed in development vs. articulations that persist in the adult; and (3) articulations that close in craniosynostosis vs. articulations that close during normal development. The effect size of the difference of means between groups in standard deviations was estimated using Cohen’s d. The statistical analysis was performed using JASP version 0.7.5.6.

We tested the null hypothesis of equal distribution between groups against the corresponding alternative hypotheses that:

articulations affected by craniosynostosis have lower reliability scores than articulations unaffected (one-sided test);

articulations that close during normal development have lower reliability than those that persist in the adult (one-sided test);

articulations affected by craniosynostosis have different reliability scores than those that close during normal development (two-sided test).

Results

The human skull at birth comprises 32 bones and 93 articulations, of which only a small fraction are associated with non-syndromic craniosynostosis conditions. We investigated the relationship between the link reliability score and the susceptibility of an articulation to close during normal development or due to craniosynostosis.

First, we compared the reliability score of those articulations that close during the normal development of the skull to those that persist in the adult. We find that sutures that normally close have significantly slightly lower reliability scores than those that do not (Mann-Whitney-Wilcoxon: one sided, W = 368, p-value = 0.047; Cohen’s d = −0.52) (Fig. 2); which is in agreement with our hypothesis that during normal development there is a tendency to close articulations that are topologically rare in the newborn skull.

Figure 2

Box plot comparing link reliability scores. Articulations associated with craniosynostosis have lower reliability than those that are not associated (left, white panel). Articulations that close during normal development also have lower reliability than those that will persist in the adult live (right, gray panel).

Next, we compared the reliability score of articulations that close prematurely in craniosynostosis to that of those articulations unaffected by this pathological condition (Fig. 2). We found that articulations associated with craniosynostosis have significantly lower reliability scores than unaffected articulations (Mann-Whitney-Wilcoxon: one-sided, W = 98, p-value = 0.006; Cohen’s d = −1.066) (Fig. 2); which shows that articulations associated to craniosynostosis are also unexpected from a topological point of view.

Interestingly, we find that while the reliability scores of articulations that close in craniosynostosis conditions tend to have lower scores than those that close during normal development, this difference is not statistically significant at a 5% significance level (Mann-Whitney-Wilcoxon: one sided, W = 15.5, p-value = 0.087; Cohen’s d = −0.964). While this marginal significance (p-value < 0.1) might be due to the reduced statistical power of two sample tests in small samples, our finding suggests that despite skull architecture being an important factor in the loss of sutures during both pathological and normal development, there are non-topological factors that further discriminate between normal and pathological loss of sutures. However, this result must be interpreted with caution due to the small sample size of both groups (N = 6 and N = 11, respectively); notice that Cohen’s d is in fact indicating a difference of means of a similar magnitude to that observed in the previous comparison (see also Fig. 2). Further details of the statistical analysis and score values are available in the Supplementary Information.

Discussion

Our results suggest that the whole arrangement of craniofacial articulations of the skull might act itself as a structural constraint, making some articulations more susceptible to closure than others. The presence of processes acting at the level of the entire skull (e.g., via bio-mechanical signaling) and predisposing bones to a premature fusion have been suggested before in the context of the functional matrix hypothesis[26]. In addition, we show that reliability scores can pinpoint articulations that are more susceptible to close during both normal and pathological development. A low reliability score identifies articulations that are ‘unexpected’ in the context of the network topology of the skull. Thus, we propose that the very arrangement of bones in the skull predisposes some articulations as targets of pathological conditions.

We are not yet in a position to offer a mechanistic explanation for the relationship reported here, which we believe may be related to the same developmental mechanism that regulate the compensatory growth of bones after premature synostoses[2, 27, 28]. However, because articulations that close during normal development also show low reliability scores (i.e., they are unexpected to occur or persist) compared to those articulations that persist in the adult skull, our findings also suggest that such mechanisms might not be different between normal development and pathological conditions. In fact, the signaling pathways in both cases are the same, notably, the FGF, TGF-β/BMP, and Wnt pathways, as well as their upstream and downstream targets. Moreover, it is known that polycistins act as mechanosensors, transducing tensile forces on the mesenchymal cells to promote osteogenesis at the cranial sutures via those pathways[29]. Thus, a possible explanation is that there is a link between the structural constraint caused by the network topology and the signaling pathways promoting osteogenesis via mechanosensors. For example, since tensile forces result from the mutual interaction among the growth fronts of each bone (i.e., the connection), the resulting connectivity pattern of the network must distribute these tensile forces in a very specific way, playing a critical role in the likelihood that some sutures, and not others, will close as the action of the mechanosensors will respond differently by activating or suppressing the osteogenic signaling pathways.

Pathological conditions of the human skull such as craniosynostosis are a medical and social problem that needs special attention from the research community. In addition, they represent medical examples of more general developmental and evolutionary processes found in all tetrapods[16, 30]. Both aspects, the medical and the biological, need and can be integrated in order to reach a better understanding that could lead to improve treatments as well as to further our knowledge about fundamental evolutionary questions. If, as our results suggest, the system of articulations of skull bones is able to self-regulate or to constrain the formation and maintenance of individual bone articulations, this might have consequences also at an evolutionary scale. In craniosynostosis conditions, the number of bones is reduced due to the early fusion of bones, much in the same way as the net reduction in the number of bones during vertebrate evolution[12, 31, 32]; as a consequence, it has been postulated that craniosynostosis could be used as an informative model for skull evolution[33]. Our results suggest that this is not a mere analogy, but that similar constraints would regulate the pattern of bone contacts in the skull, both in development and in evolution.

Supplementary Material

Supplementary Dataset 1

Supplementary Dataset 2