A comparison of metrics for quantifying cranial suture complexity

Cranial sutures play critical roles in facilitating postnatal skull development and function. The diversity of function is reflected in the highly variable suture morphology and complexity. Suture complexity has seldom been studied, resulting in little consensus on the most appropriate approach for comparative, quantitative analyses. Here, we provide the first comprehensive comparison of current approaches for quantifying suture morphology, using a wide range of two-dimensional suture outlines across extinct and extant mammals (n = 79). Five complexity metrics (sinuosity index (SI), suture complexity index (SCI), fractal dimension (FD) box counting, FD madogram and a windowed short-time Fourier transform with power spectrum density (PSD) calculation) were compared with each other and with the shape variation in the dataset. Analyses of suture shape demonstrate that the primary axis of variation captured attributes other than complexity, supporting the use of a complexity metric over raw shape data for sutural complexity analyses. Each approach captured different aspects of complexity. PSD successfully discriminates different sutural features, such as looping patterns and interdigitation amplitude and number, while SCI best-captured variation in interdigitation number alone. Therefore, future studies should consider the relevant attributes for their question when selecting a metric for comparative analysis of suture variation, function and evolution.

Traditionally, methods applied to the suture complexity problem were based on length ratio measurements. Sinuosity Index (SI) (Equation 1) (Westerman, 1971) was the earliest method, followed soon after by Suture Complexity Index (SCI) (Equation 2) (Saunders, 1995), which utilised the foundation provided by SI but created an additional complexity factor (CF) to place a greater emphasis on interdigitations. Both SI and SCI utilise Euclidean geometry methods to generate a length ratio measure of complexity. A revolution in the field of shape occurred following the pivotal work of Mandelbrot who developed fractal geometry (Mandelbrot, 1982) whereby shapes of expanding symmetry are produced from a series of iterations that appear remarkably like patterns in nature. Sutures form natural fractals, which led to the implementation of fractal dimension methods to quantify suture complexity as a statistical index that is widely applicable across diverse morphologies (Long, 1985;Boyajian & Lutz, 1992;Long & Long, 1992;Monteiro & Lessa, 2000;Skratz & Walocha, 2003). Specifically, fractal dimension measures the recurving of a suture. Sutures can be reflected as transect line data for the implementation of fractal dimension methods. One-dimensional methods of fractal dimension are reported to be suitable for transect line data, these include: box counting, Hall-Wood estimator, madogram estimator, variogram estimator, and spectral and wavelet estimator (Gneiting et al. 2012). Previous studies evaluating cranial suture complexity (Skrzat & Walocha, 2003) implemented the box count method (Equation 3). However, the madogram method (Equation 4) has since been recommended for use due to its higher efficiency and robustness (Gneiting et al. 2012). where, p is order of p of a stochastic process with stationary increments for the madogram approach p = 1; L ≥ 2; = log( / ); ̅ is the mean of , … , ; ̂( / ) is the classical method of moments estimator Equation 4 Given that sutures possess a huge degree of morphological complexity, it has been argued that previous methods utilising a single value of complexity are overly simplistic, produce ambiguous results, and are unable to distinguish between disparate morphologies (Lutz & Boyajian, 1995;Gildner, 2003). Quantifying sutures with periodic signals has been suggested to provide a more comprehensive description of suture complexity and morphology (Allen, 2006). Consequently, Fourier approaches, which summarise suture morphology as a series of sine and cosine functions (Gildner, 2003), have been suggested to form a stronger basis for studies considering morphological variation such as those utilising Elliptical Fourier Analysis (EFA) (Crampton, 1995;Tort, 2003;Carlo et al. 2011;Bardua et al. 2018, Emmons et al. 2018. EFA methods, however, deal with closed outline shapes such as an entire bone circumference. As sutures largely present as open outlines they violate the basic assumptions of EFA (leaking and aliasing) (Allen, 2006). Instead, for open outlines, it is possible to apply a windowed short-time Fourier transform (STFT) approach, which utilises coefficients calculated from a discrete Fourier transform (Equation 5) (Gildner, 2003;Allen, 2006;Wu et al. 2007). STFT obtains a rigorous mathematical decomposition of the suture morphology which can then be summarised to a single point estimate with a power spectrum density analysis (PSD) (Equation 6) (Allen, 2006). However, Fourier methods have seldom been applied to suture morphological variation thus far.
where, S is the skip frame, , is the Fourier coefficient of the k th harmonic in the m th frame (window) along the length of the original discretely sampled ( ).

Specimens
A sample of 79 specimens was generated to assess and test methods available for quantifying suture morphology. This test dataset was selected to sample a diverse range of suture morphologies and therefore encompassed a broad range of mammalian taxa, including extant eutherians (n=41), extant marsupials (n=4), extinct eutherians (n=31) and extinct marsupials (n=3) (electronic supplementary material: table S1). Multiple sutures were sampled across this sample dataset rather than a single homologous suture, as the goal was to determine the most appropriate metric for a sample with a broad range of morphologies and not to assess the morphological variation within the sample dataset.
Three-dimensional isosurfaces (n=79) were loaded in the rgl R package (Adler & Murdoch, 2019), to capture 2D images of the sutures. A single homologous suture was not used across all specimens; instead, many different sutures exhibiting a range of morphologies were captured by the 2D images (figure 1), as the purpose of the study was to determine the most appropriate metric over a range of suture morphologies and complexities.

Morphometric data collection
Open suture outlines were obtained from the 2D suture images (n=79), by the manual positioning of semi-landmarks along the individual suture outline from start to end, using the StereoMorph R package (Olsen & Haber, 2018) (electronic supplementary material: figure S1). Two-dimensional semi-landmarks were subsequently resampled at 500 per suture to ensure that suture complexity was accurately captured.

Shape analysis
To remove all non-shape elements, specifically rotation, translation, and isometric size (Rohlf & Slice, 1990) from the resampled 2D semi-landmarks, generalised Procrustes' analysis (Gower, 1975) was performed. This was implemented using the 'gpagen' function in the geomorph R package (Adams & Otárola-Castillo, 2013;Adams et al. 2019) to centre, scale and rotate the 2D semi-landmark coordinate data. The Procrustes' superimposed 2D semi-landmarks were used to produce a principal components analysis using the 'plotTangentSpace' function in the geomorph R package (Adams & Otárola-Castillo, 2013;Adams et al. 2019). This PCA was used to analyse the distribution of specimens, based on suture shape, in a multivariate morphospace for comparison to suture complexity metrics. Moreover, the PCA was used to identify the major axes of shape variation captured by the sample dataset and whether these reflected variations in complexity.

Complexity analysis
The five methods to determine suture complexity ( (n=79). The methods applied used a range of approaches to capture complexity, with SI and SCI focussing on a linear length approach, fractal dimension employing a statistical index of complexity, and PSD being calculated on Fourier-based method to capture periodic signals using the sum of trigonometric functions. All methods generated a single value indicative of complexity, with low complexity reflected by a low value and higher levels of complexity reflected by a higher value.
Sinuosity index calculations were implemented in base R (v.3.6.0; R Core Team, 2019) using a standard distance between two points equation (Equation 7) in order complete the SI equation (Equation 1, figure 2) (Westerman, 1971). Suture complexity index utilised the SI value calculated in base R (v.3.6.0; R Core Team, 2019) and a complexity factor multiplier which was calculated from the interdigitation lobes (major and minor), as outlined by Saunders (1995) (figure 2). Fractal dimension was calculated in the fractaldim R package (Sevcikova et al. 2014), using both the box counting and madogram methods by implementing the 'fd.estim.boxcount' and 'fd.estim.madogram' functions respectively. STFT was computed using the 'stft' function in the e1071 R package (Meyer et al. 2019). The PSD of each suture was calculated using the STFT results, by averaging the squared STFT coefficients over each frequency across the local transforms and summing the averages at each harmonic, as described by Allen (2006).
where dx is the distance between the x-coordinates and dy is the distance between the ycoordinates

Comparative analysis of the complexity metrics
Subsequent to the implementation of the five complexity methods (SI, SCI, FD box, FD madogram, and PSD), the appropriateness and effectivity of each metric was assessed using multiple comparative approaches. In order to identify which of the metrics captured the functionally relevant suture morphologies and the variation in complexity, comparisons were made between the shape data and complexity scores. From the PCA produced using the 2D semi-landmark data, PC scores for each PC axis capturing >5% variation were extracted for every specimen. Correspondence between complexity score and PC score were quantified with Pearson's correlation coefficient (r) and visualised with a heatmap using the 'corrplot' function from the corrplot R package (Wei et al. 2017).
The significance (p<0.05) of the correlations were tested using the 'cor.mtest' function (Wei et al. 2017). To further determine whether the complexity metrics captured the functionally relevant morphological variation of the sample dataset, a heatmap reflecting complexity scores for each specimen was mapped onto the PCA of Procrustes superimposed 2D semi-landmarks.
To assess the similarity among the five different complexity metrics, Pearson's correlation coefficient (r) was quantified and visualised as a heatmap using the 'corrplot' function in the corrplot R package (Wei et al. 2017). The significance (p<0.05) of the correlations were tested using the 'cor.mtest' function (Wei et al. 2017).
In order to assess variation and alignment among the methods across our sample, we conducted a further principal components analysis with the specimen complexity scores, rather than shape, obtained by the five complexity methods. The contribution of each variable on the PC axes (i.e. the scaled PC loadings or eigenvectors) were represented on a correlation circle, where each variable equated to one of the five complexity metrics. As each variable trajectory approaches the circle limit, the greater the contribution the method has on the complexity variation explained by the PC axis. The correlation circle was plotted using the factoextra R package (Kassambara & Mundt, 2017). A second PCA of complexity was plotted using the 'fviz_pca_biplot' function in the factoextra R package (Kassambara & Mundt, 2017), to plot both specimens (n=79) and variables (n=5) using the complexity scores, in order to summarise the relationships between the methods and specimens and therefore characterise the aspects of complexity captured by the various metric.   Figure S2. Extreme suture morphologies for the PC axes explaining >5% of the overall variation from PCA of 2D semi-landmarks, with the negative extreme indicated on the left and positive extreme indicated on the right. Figure S3. Least and most complex sutures identified by each complexity method (SI, SCI, FD box counting, FD madogram, PSD). Figure S4. PC loadings for each of the five methods (SI, SCI, FD box counting, FD madogram, PSD) on the PCA of complexity scores, for the PC axes contributing to >5% of the overall variation: (a) PC1; (b) PC2; (c) PC3.