Analysis of 3D-models of artificially deformed crania, using geometric morphometry

Artificially deformed crania are much rarer employed in craniometric studies than undeformed skulls (Tiesler, 2012: 33; 2014: 4). One of the reasons for this is the difficulty of describing varying patterns of deformation of such crania, using traditional typological schemes (Natahi et al., 2019; Kazarnitsky, Kapinus, Grigoriev, 2021). Another reason is the influence of the deformation on the initial size and shape of the crania, which obstructs the use of their data for studying population history. But the advent of novel methods of morphometric analysis—geometric morphometrics (GM) (Bookstein, 1991, 1997; Slice, 2005; Vasiliev, Vasilieva, Shkurikhin, 2018; and others)—stimulates researchers to turn back to this “inconvenient” object of study (Mayall, Pilbrow, 2018; Natahi et al., 2019; Gromov, Kazarnitsky, 2020; and others).

Unlike “traditional” morphometrics (TM), GM does not employ linear measurements, but rather

coordinates of landmarks placed on objects of any shape, including those problematic for measuring. An additional advantage of GM is the possibility to exclude the influence of the absolute dimensions, which is possible at the stage of central adjustment, scaling, and rotation. Multidimensional statistical methods, including the principal component analysis (PCA), are applicable for analyzing GM data in just the same way they are for TM data (Vasiliev, Vasilieva, Shkurikhin, 2018: 33–124).

In this study, we turn to the collections of crania of the Okunev archaeological culture, known for their artificial deformation (Zhirov, 1940; Benevolenskaya, Gromov, 1997; Gromov, 1998). The history of the study of deformation in those collections reflects the development of methods for describing the shape of the cranial vault. The Okunev culture occupied the Minusinsk basin in Southern Siberia from the 26th to 18th centuries BC (Maksimenkov, 1965; Vadetskaya, 1986: 27–36; Polyakov, Svyatko, 2009). Sites of this culture can be divided into several chronological horizons, gradually changing one another during the continuous development of the culture over a long period of time (Lazaretov, 2019). Three of those horizons contained cranial collections of a satisfactory sample size: the early Uibat (26th–25th centuries BC) and Tas-Khazaa (24th century BC), and the relatively late Chernovaya (22nd–20th centuries BC).

The first mention of the deformed crania from the Okunev Ulus cemetery was made by E.V. Zhirov (1940) in the time when the term “Okunev culture” had not yet been coined. This conclusion by Zhirov was based on visual assessment of the shape of the skull—an easy and convenient though subjective method. Further, Zhirov developed a visual scheme of classification of the types of artificial deformation (Ibid.). However, the absence of objective criteria was likely the reason why later, when the Okunev archaeological culture was singled out by G.A. Maksimenkov (1965), some researchers expressed doubts regarding the cause that produced the specific shape of the Okunev skulls (Alekseev, Gokhman, Tumen, 1987). Only a morphometric study employing metric variables describing the occipital and parietal regions of the cranial vault (indexes of the shape and height of the occiput) was able to demonstrate the artificial nature of the changes in the shape of Okunev skulls (Benevolenskaya, Gromov, 1997). But those indexes have a disadvantage of being dependent on the reference points of the coordinate system put forward by Y.D. Benevolenskaya (1976: 25–28) for describing the occipital region of the skull. If the deformed area does not fit into the region containing the reference points (i.e. the pole of the occiput and the projection of the opistion point on the sagittal suture), the identification of the deformation becomes impossible (Gromov, 2004).

The next step in studying the artificially deformed crania of the Okunev people was related to the transition from measuring single linear distances and calculating their indexes to studying the shape of the sagittal outline in general. This became possible due to the use of GM methods. These are particularly relevant, because in the Okunev skulls, the deformed area can be rather small and then only subtly affects the shape of the vault. In such cases, the severity of the deformation cannot be described via linear measurements.

Our analysis of the sagittal outlines of the deformed crania has demonstrated systematic differences in the patterns of deformation in the early and late Okunev samples (Gromov, Kazarnitsky, 2022). The deformed area in the skulls from the earlier horizons is relatively small and is placed near lambda . It affects the upper part of the occipital squama and the posterior part of the parietal bones, while it likely does not affect the height of the vault. In the later samples, the deformation is stronger pronounced and covers almost the whole parietal region, decreasing the vertical diameter of the skull. However, the variation of the sagittal outline does not reflect all the consequences of the artificial deformation. In order to obtain information regarding the variation of the cranial vault in general and to develop a protocol for analyzing 3D-models, we performed 3D-scanning of the same collections.

Material and methods

Well-preserved skulls displaying visually discernible signs of deformation were selected from the samples of the Okunev culture. The high requirements to the state of preservation have led to a reduction of the sample size as compared to our previous research of the sagittal outlines (Gromov, Kazarnitsky, 2022): 26 male and 19 female skulls were studied instead of 35 and 28, respectively. Crania from the following sites were sampled: Verkhny Askiz I, Uibat III and V, Uibat-Charkov, Itkol I and II (excavated by S.V. Khavrin, A.A. Kovalev, I.P. Lazaretov, A.V. Polyakov, S.V. Morozov). All the collections are presently kept at the Museum of Anthropology and Ethnography RAS.

The skulls were scanned using the RangeVision Spectrum device and the RangeVision ScanCenter NG 2021.2 software, using a medium-sized calibration field. The precision of placing a 3D-landmark was 0.06 mm; projector resolution 1280 × 800 pixels; camera resolution 2048 × 1536 pixels. Each 3D-model was initially exported as a cloud of 1.5 to 2.0 points in the PLY format (3–4 millions of polygons). But later, lighter models (22,500 points, 45,000 polygons) proved to be more convenient.

As our main objective was to analyze variation of the cranial vault, the facial parts of the 3D-models were removed during obtaining the sample via the tools available in RangeVision ScanCenter NG. The need for this removal was related to our intention to use the virtual digitizer of the Geomorph package in R. But later, preference was given to the Landmark IDAV software, and the step of removing the facial parts turned out to be redundant.

Further, the PLY-files were converted to the ASCII format in order to present the raw data as a text rather than a binary document. This was necessary to process the data in various programs, i.e. Landmark IDAV or R-Studio. The type of formatting was changed in MeshLab via switching off the “Binary encoding” variant of extension when exporting a file.

The shape of an object can be described in GM using two types of landmarks. The first is landmarks in a narrow sense, which are placed on homologous elements of the shape defined by strict morphological criteria. Landmarks of this type, in turn, can be divided into three categories: points at the contact between homologous elements; points of maximum curvature (concavity or convexity); and marginally delimiting points. The traditional craniometric points (Alekseev, Debets, 1964: 41–48) all belong to this first type of landmarks. The second type (semilandmarks) was developed for describing smooth spherical or extremely complex shapes, and was, thus, more suitable for our study. Semilandmarks are the points placed at an equal distance from each other along a line connecting “true” landmarks (of the first type). The more semilandmarks are placed the better their array describes the outline of an element of shape; therefore, landmarks of this type are sometimes referred to as “outline semilandmarks” (see (Vasiliev, Vasilieva, Shkurikhin, 2018: 46–53; Pavlinov, Mikeshina, 2002)).

Both types of landmarks were employed in the present study. Six anatomical landmarks and 450 semilandmarks were placed on each of the 3D cranial models in Landmark IDAV (Wiley, 2006). The following craniometric points were employed (Alekseev, Debets, 1964: 45–47): 1) frontomalare orbitale, left; 2) nasion; 3) frontomalare orbitale, right; 4) porion, right; 5) opistion; 6) porion, left. The surface delimited by those points (excluding the cranial base) was evenly filled with semilandmarks divided into two symmetric networks (patches) covering the left and right sides of the cranial vault, 225 (15 × 15) semilandmarks each (Fig. 1). The coordinates of the landmarks and semilandmarks were exported in

Fig. 1. Landmarks ( a ) and semilandmarks ( b ) on a 3D-model of the cranial vault.

1 – left patch; 2 – right patch.

the PTS format for individual skulls, and in the DTA format for the whole sample.

The following procedures of the Generalized Procrustes Analysis (GPA), including calculation of the mean shape and Procrustes coordinates, PCA, and visualization of the main trends of variation as “heat maps” were carried out in R-Studio. The R language for statistical analyses targets a wide range of users, which makes it relatively easy to master and employ. A large number of open-source program packages have been created to date, using this language. All the functions of R are available for editing and addition by any researcher. By using various packages and editing existing functions, it is possible to develop new statistical algorithms for particular research aims. There are a number of GM packages for creating databases and statistical analyses in R. Three of these were employed in this study: Morpho/Rvcg, Geomorph, and Arothron (Schlager, 2017; Adams, Otarola-Castillo, 2013; Antonio et al., 2021).

Results and discussion

Most of the analyses were carried out in the Morpho and Rvcg packages, created by the same author (Schlager, 2017) and related to each other. The first of the packages contains functions for employing the GM techniques. The second, an immediate extension of the first, can be used for modeling and surface warping. The Geomorph package (Adams, Otarola-Castillo, 2013) largely duplicates the functions of Morpho, and could be used as the main tool of the analysis. However, in this study, we employed only one function of Geomorph, absent in Morpho—computation of a mean 3D-model.

The Arothron package (Antonio et al., 2021) was rather an accessory tool, containing some additional functions widening the possibility of presenting and visualizing the results. This package was specifically developed for visual presentation of anthropological objects. Some of its functions are related to working with 3D-models, virtual reconstructions, and restoration of skeletal elements, while others facilitate the visual representation of the results. One of the functions can be used for creating “heat maps”, which demonstrate local differences in shape via mapping some areas of the object as “warm” (mainly orange) or “cold” (mainly blue) colors.

The algorithm of working in R we followed in this study can be broadly divided into three stages: import of the data, statistical analysis, and visualization of the results in the form of plots and “heat maps”. Morpho and Geomorph were used for the first two stages, while some functions of Arothron and other graphic editors were employed at the third stage.

The read.ply and read.pts functions were utilized for importing the data to the workspace of R. These functions upload PLY-files of 3D-models and PTS-files containing landmark and semilandmark coordinates. These raw data were then transformed to arrays using the list2array function of Morpho in order to prepare the data for further operations.

The mean shape (i.e. the configuration of landmarks describing the average shape of a particular sample) was calculated at the second stage by the ProcGPA function of Morpho. It was extracted using the aggregate and vecx functions, and then a PCA was carried out employing the procSym function of the same package. The PCA values for individual skulls, as well as the proportion of the total variance described by the PCs, were exported for drawing scatterplots in an external graphics editor.

At the last stage, the “heat maps” visualizing the morphological changes associated with particular PCs were created. The areas displaying the highest level of morphological variation are depicted in the 3D-models via the gradations of color (Fig. 2–4), ranging from red and orange (maximum and intermediate expression of a trait—the degree of convexity of an area of the cranial vault in this particular case) to violet and blue (weak or the weakest convexity of an area). The areas exhibiting no shape change were colored in white.

Building a “heat map” requires an averaged 3D-model, and coordinates of the mean shape and the individuals displaying extreme values of the PCs. The averaged 3D-models were created in Geomorph, using the model of one of the skulls in the sample whose shape was the closest to the mean shape. Such skulls have been detected using the findMeanSpeс function, and then the averaged 3D-model was calculated by warpRefMesh . The extraction of the extreme morphological variants was beginning with the creation of the shapes representing the maximal and minimal values of the PCs, using the restoreShapes function in Morpho. These were further transformed employing the thin-plate spline interpolation (Morpho, tps3d ). When all the necessary data were extracted, the “heat maps” representing the “poles” of the morphological variation in the sample were created using the localmeshdiff function of Arothron.

The algorithm outlined above was employed for exploring the trends in variation of the 3D-models of the deformed crania from the Okunev burials.

Fig. 2. Principal component analysis of the male sample (PC1 and PC2).

a – Uibat chronological horizon; b – Tas-Khazaa chronological horizon; c – Chernovaya chronological horizon.

Fig. 3. Principal component analysis of the male sample (PC2 and PC3). Legend same as on Fig. 2.

Fig. 4. Principal component analysis of the female sample (PC1 and PC2).

Legend same as on Fig. 2.

The first three PCs of the analysis of the Procrustes coordinates in the male sample account for 38 % of the total variance: PC1 – 15.4 %; PC2 – 12.6 %; PC3 – 10.1 %. The “heat maps” of the cranial models are depicted near the respective axes of the plot (Fig. 2–4) in four projections: left and right lateral, frontal, and occipital. The maps show the shapes associated with minimal negative and maximal positive values of the PCs.

The “heat maps” show that PC1 (see Fig. 2) describes the change from relatively tall vaults (with a more prominent parietal area) to relatively low crania (more bossing in the frontal and occipital parts). PC2 (see Fig. 2, 3) is associated with the variation of the sagittal outline and the width of the vault: a more convex sagittal outline (according to the position of the orange zones along the sagittal axis) is correlated with flatter lateral walls of the vault, and vice versa . These two PCs (see Fig. 2) together account for almost a third of the total variance, and demonstrate the difference between the individuals belonging to the early and late Okunev burials. The former (Uibat and Tas-Khazaa) are usually relatively tall, while the latter (Chernovaya) exhibit a relatively low cranial vault and more convex temporo-parietal areas.

The degree of occipito-parietal flattening is described by PC3 (see Fig. 3): this feature ranges from weak flattening, accompanied by a relatively high parietal curvature, to strong flattening, with more bossing frontal bone and the lower part of the occipital squama. The skulls of the three chronological horizons display the clearest separation in the plot (see Fig. 3): not only are the early groups different from the late groups, but also the two early samples are distinct. The chronologically earliest individuals from Uibat exhibit a combination of the least pronounced occipital-parietal flattening and the least convexity of the lateral walls with a tall cranial vault.

The analysis of the female sample has shown a clear separation of the skulls from different chronological horizons along the first PC, which accounts for 21 % of the total variance (see Fig. 4). This axis describes the variation trend from tall cranial vaults with the least flattening of the parietal bones and the upper part of the occipital squama (with the lateral walls of the walls being vertically oriented) towards relatively low crania displaying a weaker height of the curvature of the vertex (with the temporo-parieatal areas being more convex or bossing). The individuals from Uibat occupy the former pole of PC1, while the skulls from

Chernovaya occupy the latter. The two chronologically intermediate individuals from Tas-Khazaa appeared to be morphologically intermediate as well.

Thus, our analysis of 3D-models of the deformed crania of the Okunev people has demonstrated a gradual change in their morphology at the transition from the early to late periods of the development of the culture, which was likely associated with a modification of the construction of the cradle and/or some accessories placed therein (Benevolenskaya, Gromov, 1997: 293). The pattern of morphological changes with time is the same for males and females.

Conclusions

The analytical algorithm put forward in the present study fosters the application of GM methods for the analysis of 3D-models with free-ware R packages and the Landmark software. The advantages of the algorithm were demonstrated using an example of a sample of deformed crania of bearers of the Okunev archaeological culture. We can recommend using this algorithm in morphological studies of any objects, both biological and archaeological.

Our analysis has shown the differences in the localization of the deformed areas in the Okunev skulls from different chronological horizons. The study of the 3D-models has demonstrated that the increase in the severity of the parietal deformation (which is typical of the individuals from later burials) led to a decrease in the height of the cranial vault due to a flattening of its parietal region, and to an increased convexity of the lateral walls of the skull, frontal bone, and lower part of the occipital squama. Nevertheless, most of the features of shape differentiating the skulls of the early and late horizons reside in the mid-sagittal plane of the vault.