Morphometric Maps of Bilateral Asymmetry in the Human Humerus: An Implementation in the R Package Morphomap

: In biological anthropology, parameters relating to cross-sectional geometry are calculated in paired long bones to evaluate the degree of lateralization of anatomy and, by inference, function. Here, we describe a novel approach, newly added to the morphomap R package, to assess the lateralization of the distribution of cortical bone along the entire diaphysis. The sample comprises paired long bones belonging to 51 individuals (10 females and 41 males) from The New Mexico Decedent Image Database with known biological proﬁle, occupational and loading histories. Both males and females show a pattern of right lateralization. In addition, males are more lateralized than females, whereas there is not a signiﬁcant association between lateralization with occupation and loading history. Body weight, height and long-bone length are the major factors driving the emergence of asymmetry in the humerus, while interestingly, the degree of lateralization decreases in the oldest individuals.


Introduction
In bioarchaeology and anthropology, it is of interest to infer the physical activities, occupations and behaviours of past populations from skeletal material [1,2]. During life, the distribution of cortical bone is influenced by loading history [3][4][5], and bone remodelling seems to be significantly associated with high-frequency daily action [6]. Asymmetry of loading, as is common in many physical activities, occupations and behaviours, can be expected to lead to asymmetry of bone form. Thus, to fulfil the goal of inferring past lifestyles often requires the assessment of differences in bone shape and cortical thickness distributions among antimeres [5,7].
Different models have been proposed to explain how bone is remodelled in relation to loading [8][9][10][11], although bone adaptation and remodelling has a sizeable physiological and environmental (i.e., nutritional) component. The comparison of antimeric bones from the same individual offers the opportunity to identify asymmetry of loading history [12] while ignoring the confounding, presumed bilaterally equal effects of genetics and nutrition. Yet, even the comparison of paired bone elements is not entirely without issues, since inflammatory processes [13] may trigger osteogenesis in distant regions [14,15], and differences in patterns of asymmetry in the upper limb have been found with ageing [16] and long-term disuse [17], in addition to loading history per se. Despite these caveats, traditional methods that rely on calculations of the percentage change of cross-sectional geometric parameters (total area, cortical area, area moments of inertia) on the humerus have provided useful insights into activity patterns in modern [12] and archaeological populations [18][19][20][21][22][23][24], as well as in paleontological samples [25][26][27][28][29][30]. Studies of professional athletes who play unimanual or bimanual sports, such as tennis [5,[31][32][33][34], throwing and swimming [34][35][36], provide an interesting natural experiment. Studies of their long bones allow assessment of the extent to which asymmetry of cortical thickness and whole-bone morphology exists between the dominant (e.g., the racket arm) and non-dominant arm. Younger starters show a higher index of "strength" than older, suggesting that intense activities during adolescence lead to greater subperiosteal expansion [37].
Current methods to evaluate asymmetry in long bones often involve comparison of shape and biomechanical parameters (cross sectional geometry) between antimeres based on a limited number of sections along the diaphysis [19,38,39]. More recently, Wei et al. [40] have extended such analyses to multiple closely placed sections along the entire diaphysis, calculating asymmetry of cortical thicknesses and polar moments of area (J) using the R software tool, morphomap [41].
Here, we assess how asymmetry in the distribution of cortical thickness in the human humerus is related to physical activity level, sex, body mass and weight based on data from recently deceased individuals, with known occupational, lifestyle and medical history curated in the New Mexico Decedent Image Database (NMDID) [42].
We tested the hypotheses that: (a) the degree of asymmetry does not differ among sexes or among three different occupation groups; (b) the difference in distribution of cortical bone and degree of asymmetry are not influenced by age, weight, height, humerus length and occupation. The hypotheses we tested have significant implications for the evaluation of asymmetry in archaeological populations and in extinct human species.

Data Preparation and Processing
From the New Mexico Decedent Image Database [42], we selected 51 individuals (41 males, 10 females) with known occupation, ranging in age at death from 21 to 54 years. We selected individuals who had been in the armed forces or worked in building/mining or at a desk job to test the methodology using groups with distinctly different occupational histories (likely high vs. low loading).
We extracted from NMDID metadata associated with occupational history for each individual. Then, we computed occupation scores relating activity to energy cost (see Appendix 4.1 from [43]) and duration in years for each occupation. Missing data for duration in years are estimated by calculating the expected working based on the formula, (age at death -18 years)/ the number of recorded occupations.
A total body CT scan is available for each individual at a slice resolution of 0.5 mm with 16 × 0.75 mm collimation, 120 KVp and 300 mAs. From these scans, we cropped the left and right humerus. In order to create 3D models defined by only bony material, the image stacks were automatically segmented using the Otsu algorithm available in the morphomap R package. The 102 resulting 3D models (51 left humeri and 51 right humeri) were aligned following the protocols proposed by Ruff [44]. From each 3D model, we extracted 61 cross sections from 20% to 80% of the biomechanical length along the bone shaft. At each cross section, we defined 24 paired equiangular semilandmarks on the external and internal outlines centred at the barycentre of the cross section. The production of the cross sections is automatically executed in morphomap by using the functions morphomapCore and morphomapShape ( Figure 1).

Asymmetry and Cross-Sectional Geometry
On each individual, we calculated the polar moment of inertia (J mm 4 ) at 40% of the biomechanical length on both sides using the function morphomapCSG. We avoided standardization of J (on body mass and bone length), because we analyzed the percentage of lateralization (JLAT%) using the following equation: JLAT% = (|(JR − JL)|/JM) × 100, where JM = (JR + JL)/2.

Description of the Function MorphomapAsymmetry
The new function, morphomapAsymmetry, embedded in morphomap facilitates the mapping and analysis of bilateral asymmetry in long bones (Table 1). We provide three strategies to map differences in the distribution of cortical thickness between the two sides: (i) the difference between sides (type="diff"); (ii) the difference from the mean (type="onMean"); (iii) the relative change in thickness (type="relChange") of one side with respect to the other.
The workflow implemented in morphomap is as follows: 1. Load the output of the first long bone processed with morphomapShape. 2. Load the output of the second long bone processed with morphomapShape ( Figure 1). 3. Specify if one of the two input objects needs to be mirrored ( Figure 2A). 4. Calculate the cortical thickness map of the entire diaphysis in both long bones ( Figure   2B). 5. Standardize the cortical thicknesses by dividing the matrices of cortical thickness by the biomechanical length (optional). 6. Choose the method of visualization by setting the argument type to: a. type="diff"

Asymmetry and Cross-Sectional Geometry
On each individual, we calculated the polar moment of inertia (J mm 4 ) at 40% of the biomechanical length on both sides using the function morphomapCSG. We avoided standardization of J (on body mass and bone length), because we analyzed the percentage of lateralization (J LAT %) using the following equation: J LAT % = (|(J R − J L )|/J M ) × 100, where J M = (J R + J L )/2.

Description of the Function MorphomapAsymmetry
The new function, morphomapAsymmetry, embedded in morphomap facilitates the mapping and analysis of bilateral asymmetry in long bones (Table 1). We provide three strategies to map differences in the distribution of cortical thickness between the two sides: (i) the difference between sides (type = "diff"); (ii) the difference from the mean (type = "onMean"); (iii) the relative change in thickness (type = "relChange") of one side with respect to the other. type Defines the method to calculate the differences in cortical thickness between the two long bones: "diff" a map of arithmetic difference between reference and target is computed; "onMean" the morphometric map of asymmetry is defined by computing the differences from the mean for each long bone; "relChange" the morphometric map is computed by calculating the relative change in cortical thickness, expressed as percentage difference between reference and target long bones reference The workflow implemented in morphomap is as follows: 1.
Load the output of the first long bone processed with morphomapShape.

2.
Load the output of the second long bone processed with morphomapShape ( Figure 1).

3.
Specify if one of the two input objects needs to be mirrored ( Figure 2A).

4.
Calculate the cortical thickness map of the entire diaphysis in both long bones ( Figure 2B).

5.
Standardize the cortical thicknesses by dividing the matrices of cortical thickness by the biomechanical length (optional). 6.
Choose the method of visualization by setting the argument type to: a. type = "diff" Calculate the differences between the cortical thickness maps of the two long bones ( Figure 2C). b. type = "onMean" Calculate the differences between the two cortical maps and their mean ( Figure 2E). c. type = "relChange" Compute a cortical map as the percentage change of one side (target) with respect to the other one (reference) ( Figure 2D).

7.
2D Plot the map of differences in cortical thickness between the selected specimens. The difference map is displayed after "unrolling" the long-bone shaft to produce a 2D plot, starting and ending at the anterior (A) border passing through the lateral (L), posterior (P), and medial (M) borders ( Figure 2C). Calculate the differences between the cortical thickness maps of the two long bones ( Figure 2C). b. type="onMean" Calculate the differences between the two cortical maps and their mean ( Figure  2E). c. type="relChange" Compute a cortical map as the percentage change of one side (target) with respect to the other one (reference) ( Figure 2D). 7. 2D Plot the map of differences in cortical thickness between the selected specimens.
The difference map is displayed after "unrolling" the long-bone shaft to produce a 2D plot, starting and ending at the anterior (A) border passing through the lateral (L), posterior (P), and medial (M) borders ( Figure 2C).  First long bone processed with morphomapShape

Description of the Function MorphomapMapPCA
We processed the right and left humerus in 51 individuals selected from the NMDID using the R package morphomap (Profico et al. 2021). The individuals belong to three different categories for occupation: "building-mining" (called "building" from now on), "army" and "desk". We extracted 61 cross sections from the humeri and built a multivariate dataset of cortical thicknesses along the entire diaphysis on both sides.
To decompose the total variance of the sample into symmetric and asymmetric components, we performed two different Principal Component Analyses (PCA) on each dataset: 1.
PCA of the mean morphometric maps calculated by averaging left and mirrored right side (symmetric component).

2.
PCA of the matrices obtained by in each individual subtracting the mean matrix of cortical thicknesses from the matrices of cortical thicknesses of the left and mirrored right sides (asymmetric component).
The function morphomapPCA requires two inputs, the left and right sets of long-bone semilandmarks, obtained from morphomapShape. The user can select if the calculation of the symmetric and asymmetric component is performed on semilandmark coordinates or on the values of cortical thickness computed from these along the diaphysis.

Relation between Cortical Thickness Asymmetry Humerus and Biological Variables
Commonly, some limitations apply in evaluating patterns of lateralization (i.e., asymmetry). Analyses are usually limited to a single (e.g., at 40% of the total biomechanical length) or a few cross sections. In addition, the investigation is restricted to the use of univariate and exploratory statistics. Here, we present two different strategies to evaluate the relative contribution to asymmetry of different predictors (i.e., weight, height, age and occupation).
To assess how asymmetry in the distribution of cortical thickness varies in relation to occupation, age, weight, height and biomechanical length, we performed a multiple regression with these variables as independent and maps of differences in cortical thickness between the left and right side as dependent variables. Specifically, the cortical maps of differences between sides are created by subtracting the mean matrix of cortical thicknesses from the matrices of the left and mirrored right sides from each linear regression we computed R 2 and beta coefficients. R 2 quantifies the strength of the relationship between the model and the dependent variable. The beta coefficient describes the rate of change of differences in cortical thickness between sides for every unit of change in the independent variables. In addition, we measured the proportion of variance in total asymmetry related to each independent variable using multivariate regression analysis. Lastly, we applied the variance partitioning method [45] to measure the portions of variance of total asymmetry shared by independent variables. The method calculates the explanatory power of different variables in relation to the same response variable (or matrix). We used redundancy analysis to determine the partial effect of each variable (i.e., weight, height, age and occupation) on the response variable (magnitude of asymmetry of cortical thickness between sides). We used alpha (significance) level = 0.05 for all statistical tests.

Results
The asymmetry in polar moment of inertia, J, calculated at 40% of the biomechanical length shows a general trend of lateralization ranging between 0.36% and 10.37%. In all but 4 individuals, J is larger on the right side (Table 2). There are no statistically significant differences between occupation and sex group means (Figure 3), as determined by two-way ANOVA of J among sexes or occupations.
The first two PCs of the symmetric component of cortical thicknesses account for 78% of the total variance (PC1 = 72.33%; PC2 = 5.79%) (Figure 4). On average, the two sexes are separated along PC1 with the females towards the positive limit, and males, the negative (t 14.48 = 4.66, p < 0.01). The three occupation groups largely overlap but with the "building", "army" and "desk" groups approximately distributed in this order from negative to positive limits of PC1. Table 2. Description of the sample and calculation of lateralization. For all the individuals, we report age, weight, height and occupation. We calculated the biomechanical length of the humerus and the polar moment of inertia, J, at 40% of the biomechanical length on both sides (JL and JR) and the degree of lateralization, J LAT %, between the two sides. Values of J are multiplied by 10 3 (Ruff 2000). D = working at desk job; B = working in building/mining companies; A = working in the armed forces.  The first two PCs of the symmetric component of cortical thicknesses account for 78% of the total variance (PC1 = 72.33%; PC2 = 5.79%) (Figure 4). On average, the two sexes are separated along PC1 with the females towards the positive limit, and males, the negative (t14.48 = 4.66, p < 0.01). The three occupation groups largely overlap but with the "building", "army" and "desk" groups approximately distributed in this order from negative to positive limits of PC1.

ID
The visualisations of the morphometric maps represented by the extremes of PCs 1 and 2 highlight a general increase in cortical thickness at negative values of PC1. PC2 represents a different pattern of thickening/thinning of the humeral diaphysis. With increasingly negative values of PC2, the cortical bone is thicker in the mid and distal portion of the medial and anterior margins and between the posterior and lateral margins. Conversely, with more positive values, the proximal portion of the diaphysis is thicker anteriorly ( Figure 4).  The PCA of the asymmetric component ( Figure 5) indicates how the cortical thickness of the entire diaphysis differs from symmetry. The distance of points from the origin indicates the degree of asymmetry represented by the first two PC scores. Following Mardia et al. [46], we found that directional (mean difference between sides) and fluctuating (average differences from the symmetric mean) components account for 16.00% and 84.00% of the total variance, respectively. PC1 explains 41.05% of the total variance. This axis describes generalised thickening or thinning of the cortex as seen in the morphometric maps representing the extremes of this PC. Scores on PC1 indicate that the right-sided cortex tends to be thicker than the left (with a few exceptions, plausibly explained by handedness). PC2 (7.47% of the total variance) shows a different pattern of asymmetry. The morphometric maps indicate that, from  Conversely, with more positive values, the proximal portion of the diaphysis is thicker anteriorly (Figure 4).
The PCA of the asymmetric component ( Figure 5) indicates how the cortical thickness of the entire diaphysis differs from symmetry. The distance of points from the origin indicates the degree of asymmetry represented by the first two PC scores. Following Mardia et al. [46], we found that directional (mean difference between sides) and fluctuating (average differences from the symmetric mean) components account for 16.00% and 84.00% of the total variance, respectively. A variance partitioning analysis was performed to evaluate the percentage of variance in asymmetry associated with weight, height, age and occupation score. The combination of all tested variables explained 24.72% of the variance in asymmetry calculated from the PC scores of asymmetric component. Weight (2.09%) and height (2.39%), and  PC1 explains 41.05% of the total variance. This axis describes generalised thickening or thinning of the cortex as seen in the morphometric maps representing the extremes of this PC. Scores on PC1 indicate that the right-sided cortex tends to be thicker than the left (with a few exceptions, plausibly explained by handedness). PC2 (7.47% of the total variance) shows a different pattern of asymmetry. The morphometric maps indicate that, from positive to negative limits, this PC represents posterior and anterior thinning of the diaphysis ( Figure 5).
While the differences between sexes (F = 18.40, Df = 1, p < 0.01) in asymmetry are significant, as indicated by two-way ANOVA (Figure 6), there are no statistically significant differences in asymmetry between the occupation groups. cortical thickness measured at each semilandmark). At each cell, the explained variance (R 2 ) and slope (Beta coefficient) can be mapped to visualise the relationship between cor tical map asymmetry and the independent variables. Maps can be drawn to represent the (exactly opposite) effects of these relationships on the right or left sides. Such maps are presented in Figure 8. The maps of R 2 indicate that each independent variable is associated with a different pattern of asymmetry, with different localised regions showing an asso ciation with each variable. Unsurprisingly, regressions of height and biomechanica length produce the most similar diagrams. On average, the slopes (beta coefficients) o asymmetry of cortical thickness on the independent variables indicate that cortical thick nesses tend to be greater in the right arm (likely this is a mostly right-handed sample) The maps of figure 8 show values only for those regions where the regression is signifi cant. Weight, height and mechanical length are associated with a larger rate of increase in cortical thickness in the right (usually the dominant) compared to the left humerus. In contrast, with increasing age, cortical thickness decreases more slowly in the dominan arm. Figure 6. Violin plots of the total length of the displacement vectors of asymmetry in relation to sex (left) and occupation (right, A = "army", B = "building", D = "desk"). Violin plots illustrate the density distribution of the data. Figure 6. Violin plots of the total length of the displacement vectors of asymmetry in relation to sex (left) and occupation (right, A = "army", B = "building", D = "desk"). Violin plots illustrate the density distribution of the data.
The lengths of the vectors in Figure 5 indicate the magnitude of asymmetry of cortical thickness, represented by the first two PC scores. The sum of the vectors from the entire matrix of PC scores represents the overall magnitude of asymmetry. In this sample, it is correlated with the index of lateralization (J LAT %), calculated using the polar moment of inertia, J, (correlation = 0.58, p < 0.01).
A variance partitioning analysis was performed to evaluate the percentage of variance in asymmetry associated with weight, height, age and occupation score. The combination of all tested variables explained 24.72% of the variance in asymmetry calculated from the PC scores of asymmetric component. Weight (2.09%) and height (2.39%), and their interaction (16.31%) explain the largest portion of asymmetry (Figure 8).

Discussion and Conclusions
Studies of patterns of lateralization in archaeological populations often suffer from the lack of a reference sample with known loading history to contextualise the findings. Additionally, the assessment of lateralization is commonly limited to the analysis of one or just a few levels along the diaphysis of paired long bones. The publication of the NMDID [42] offers the prospect of directly relating skeletal form (total body CT-scan) with known biological profile and loading history (metadata with 60 variables). The time and effort in gathering skeletal data is much reduced by the functions available in the R package morphomap [41], a recently published toolkit providing functions to extract from CT data, the segmented long bone of interest and based on that. Here, we further extend morphomap to visualize and analyze asymmetry in paired long bones. Specifically, the new implementation: (i) performs a PCA on the symmetric and asymmetric component of form variation; (ii) creates morphometric maps of symmetric and asymmetric variation on single individuals or on entire samples from PC scores; (iii) calculates the total magnitude of asymmetry of cortical bone distribution, quantifying deviations from symmetry.
PCAs of both symmetric and asymmetric components indicate that cortical bone distribution differs between sexes, but not between the occupation groups considered in our analyses. On average, male individuals possess thicker and more asymmetric humeri than females. All our measurements of magnitude of asymmetry (cross-sectional geometry and vector lengths from PC scores) present a general pattern of right lateralization. This finding is consistent with previous studies indicating 90% preference for right-handedness in modern humans [47][48][49]. Further, the analyses of morphometric maps indicate that males are more asymmetric than females. However, since males are larger on average, their more asymmetric cortical thickness might be guided by allometric effects. In contrast, the index of lateralization based on cross-sectional geometry (JLAT%) does not significantly differ between sexes. This contrast may be due to the fact that JLAT% is calculated at a single Multivariate regression was used to assess the relationship between asymmetry in cortical thickness and the independent variable of interest (i.e., weight, age, height, biomechanical length and occupation scores) at each cell of the morphometric maps (i.e., the cortical thickness measured at each semilandmark). At each cell, the explained variance (R 2 ) and slope (Beta coefficient) can be mapped to visualise the relationship between cortical map asymmetry and the independent variables. Maps can be drawn to represent the (exactly opposite) effects of these relationships on the right or left sides. Such maps are presented in Figure 8. The maps of R 2 indicate that each independent variable is associated with a different pattern of asymmetry, with different localised regions showing an association with each variable. Unsurprisingly, regressions of height and biomechanical length produce the most similar diagrams. On average, the slopes (beta coefficients) of asymmetry of cortical thickness on the independent variables indicate that cortical thicknesses tend to be greater in the right arm (likely this is a mostly right-handed sample). The maps of Figure 8 show values only for those regions where the regression is significant. Weight, height and mechanical length are associated with a larger rate of increase in cortical thickness in the right (usually the dominant) compared to the left humerus. In contrast, with increasing age, cortical thickness decreases more slowly in the dominant arm.

Discussion and Conclusions
Studies of patterns of lateralization in archaeological populations often suffer from the lack of a reference sample with known loading history to contextualise the findings. Additionally, the assessment of lateralization is commonly limited to the analysis of one or just a few levels along the diaphysis of paired long bones. The publication of the NMDID [42] offers the prospect of directly relating skeletal form (total body CT-scan) with known biological profile and loading history (metadata with 60 variables). The time and effort in gathering skeletal data is much reduced by the functions available in the R package morphomap [41], a recently published toolkit providing functions to extract from CT data, the segmented long bone of interest and based on that. Here, we further extend morphomap to visualize and analyze asymmetry in paired long bones. Specifically, the new implementation: (i) performs a PCA on the symmetric and asymmetric component of form variation; (ii) creates morphometric maps of symmetric and asymmetric variation on single individuals or on entire samples from PC scores; (iii) calculates the total magnitude of asymmetry of cortical bone distribution, quantifying deviations from symmetry.
PCAs of both symmetric and asymmetric components indicate that cortical bone distribution differs between sexes, but not between the occupation groups considered in our analyses. On average, male individuals possess thicker and more asymmetric humeri than females. All our measurements of magnitude of asymmetry (cross-sectional geometry and vector lengths from PC scores) present a general pattern of right lateralization. This finding is consistent with previous studies indicating 90% preference for right-handedness in modern humans [47][48][49]. Further, the analyses of morphometric maps indicate that males are more asymmetric than females. However, since males are larger on average, their more asymmetric cortical thickness might be guided by allometric effects. In contrast, the index of lateralization based on cross-sectional geometry (J LAT% ) does not significantly differ between sexes. This contrast may be due to the fact that J LAT% is calculated at a single level along the diaphysis (40%v of bone biomechanical length), whereas the calculation of absolute lateralization from PC scores takes into account the entire diaphysis. In fact, J LAT% calculated at 70%, 75%, 76%, 78% and 79% of bone length is statistically different between sexes.
Regression analyses on morphometric maps show that, as body weight, height and longbone biomechanical length increase, so does asymmetry. In contrast, with increasing of age, asymmetry decreases (i.e., the oldest individuals are less asymmetric than the youngest individuals). Interestingly, loading history (occupation scores) does not affect the pattern of asymmetry. The main effect of body proportions (weight and height) as a key factor in determining the degree of asymmetry is also confirmed by the partitioning of variance analysis performed on values of total asymmetry calculated from the PCA.
This study was able to test the hypothesised relationships between loading history and cortical bone distribution in the humeral shaft using the unique and extensive collection curated in the NMDID [42]. Despite the quality of these data, our analyses do not show a significant association between occupation and asymmetry. This analysis was confined to occupations that would be expected to lead to extremely different occupational loading histories (army, building vs. desk) in order to emphasise any effect of occupation. This suggests either that there is no strong difference in the effect of these occupations on loading history, or that occupational history does not reflect the full loading history, and that our categorisation by occupation is inadequate to describe individual loading history. Further studies are required to clarify this finding, which is potentially of great importance in archaeological and forensic contexts.