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

Abstract

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 profile, 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 significant 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.

1. 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.

2. Materials and Methods

2.1. 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).

2.2. Asymmetry and Cross-Sectional Geometry

On each individual, we calculated the polar moment of inertia (J mm⁴) 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.

2.3. 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. morphomapAsymmetry: description of the main arguments.

Argument	Definition
mshape1	First long bone processed with morphomapShape
mshape2	Second long bone processed with morphomapShape
standandize	If TRUE, the matrices of cortical thickness are standardized using the average biomechanical length between sides.
plot	If TRUE, the map of cortical thickness asymmetry is returned.
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	If set to 1, mshape1 is defined as reference; if set on 2, mshape2 is defined as reference.
rem.out	If TRUE, outliers are removed from the matrices of cortical thickness.
scale	If TRUE, the matrices of cortical thickness are scaled from 0 to 1.
gamMap	If TRUE, gam smoothing is applied.

). 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:

• Load the output of the first long bone processed with morphomapShape.
• Load the output of the second long bone processed with morphomapShape (Figure 1).
• Specify if one of the two input objects needs to be mirrored (Figure 2A).
• Calculate the cortical thickness map of the entire diaphysis in both long bones (Figure 2B).
• Standardize the cortical thicknesses by dividing the matrices of cortical thickness by the biomechanical length (optional).
• Choose the method of visualization by setting the argument type to:
  • type = “diff”
    Calculate the differences between the cortical thickness maps of the two long bones (Figure 2C).
  • type = “onMean”
    Calculate the differences between the two cortical maps and their mean (Figure 2E).
  • type = “relChange”
    Compute a cortical map as the percentage change of one side (target) with respect to the other one (reference) (Figure 2D).
• 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).

2.4. 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:

• PCA of the mean morphometric maps calculated by averaging left and mirrored right side (symmetric component).
• 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.

2.5. 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² and beta coefficients. R² 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.

3. 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. 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³ (Ruff 2000). D = working at desk job; B = working in building/mining companies; A = working in the armed forces.

ID	Sex	Age	Weight	Height	Occupation	JL	JR	JLAT%
100221	male	34	91	188	D	2.84	3.60	5.90
101358	male	21	70	168	B	1.50	1.57	1.07
101510	male	26	86	195	A	2.39	2.96	5.33
102253	male	46	100	183	D	1.62	1.92	4.26
102436	male	37	83	183	D	3.10	3.22	0.91
102602	male	32	86	193	B	2.26	2.40	1.56
103530	male	22	73	168	B	1.51	1.72	3.26
103862	male	31	91	191	A	2.11	2.26	1.76
104373	male	34	86	170	B	2.47	2.56	0.86
108039	male	38	82	188	A	1.65	1.92	3.77
114405	male	27	79	178	B	1.62	1.97	4.81
116546	male	25	86	184	A	2.82	2.43	3.72
116833	male	45	109	183	A	2.40	2.61	2.08
117662	female	24	68	173	A	1.29	1.50	3.85
118646	male	26	84	180	A	1.92	2.84	9.63
121289	male	25	102	184	B	2.74	2.98	2.10
123096	male	24	75	184	B	2.81	3.14	2.78
123240	male	31	82	180	D	1.64	2.05	5.61
125527	female	26	50	157	D	0.69	0.79	3.35
127137	female	54	59	158	D	1.38	1.32	1.20
129131	male	33	79	173	A	1.80	1.87	0.93
129352	male	30	82	180	A	1.56	1.86	4.29
130388	male	37	68	168	A	1.67	1.82	2.16
130964	male	25	68	175	D	2.87	3.42	4.42
132233	male	29	80	166	B	1.77	2.13	4.64
132433	male	27	70	185	D	2.06	2.30	2.72
139871	female	49	61	157	B	1.49	1.56	1.27
140368	male	34	80	185	B	3.17	3.83	4.74
141318	male	32	77	178	B	2.04	2.24	2.32
143365	female	47	72	165	A	1.16	1.10	1.32
143984	male	30	91	183	D	2.51	2.63	1.17
144071	female	35	57	160	A	1.04	1.08	0.94
144977	male	24	75	175	D	2.06	2.24	2.09
146626	male	27	89	163	A	0.98	1.16	4.24
147949	male	29	91	178	A	1.82	2.63	9.15
150608	male	34	63	158	B	1.81	2.01	2.64
152567	male	24	82	185	A	2.46	2.74	2.70
156886	male	29	77	178	B	1.78	1.97	2.51
158402	male	35	109	185	B	2.32	2.37	0.58
162065	male	28	64	170	B	1.56	1.88	4.76
166116	female	36	56	162	B	1.07	1.17	2.22
171170	male	32	106	182	A	2.96	3.75	5.90
171479	male	27	100	185	D	2.50	2.82	3.06
173218	male	22	84	163	D	1.95	2.06	1.36
175725	male	28	86	196	B	3.12	3.46	2.54
176660	female	24	68	173	A	0.79	0.84	1.47
177679	male	38	77	168	A	2.75	2.82	0.68
178078	male	31	73	168	B	2.10	2.31	2.32
180030	male	40	82	180	A	2.06	2.33	3.05
188902	female	26	77	165	D	1.25	1.12	2.82
190756	female	24	48	162	D	0.84	0.98	3.82

). 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.

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 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.

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).

In multivariate regressions, body weight (R² = 0.19, b = 0.22, p < 0.01), height (R² = 0.22, b = 0.31, p < 0.01) and biomechanical length (R² = 0.16, b = 0.14, p < 0.01) significantly predicted asymmetry, while age and occupation score are not statistically significant predictors (Figure 7).

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).

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²) 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² 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.

4. 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.