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.

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 .
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 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).
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  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 ij 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 ij given the observed network N O (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 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, Q) is completely determined by a partition P of bones into groups and the group-to-group interaction probability matrix Q. 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 σ i and j to group σ j we have that 24 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 σi σj = n 1 σi σj + n 0 σi σj is the total number of possible sutures between groups σ i and σ j , and H(P) is a function that depends on the partition only ∑ = + + .
In practice, it is useful to "thin" the sample {P i }, 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 −H(P) /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: 1. articulations affected by craniosynostosis have lower reliability scores than articulations unaffected (one-sided test); 2. articulations that close during normal development have lower reliability than those that persist in the adult (one-sided test); 3. 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.
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.