**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–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–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,

**Citation:** Profico, A.; Zeppilli, C.; Micarelli, I.; Mondanaro, A.; Raia, P.; Marchi, D.; Manzi, G.; O'Higgins, P. Morphometric Maps of Bilateral Asymmetry in the Human Humerus: An Implementation in the R Package Morphomap. *Symmetry* **2021**, *13*, 1711. https://doi.org/10.3390/ sym13091711

Academic Editor: Antoine Balzeau

Received: 19 August 2021 Accepted: 10 September 2021 Published: 16 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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–24], as well as in paleontological samples [25–30]. Studies of professional athletes who play unimanual or bimanual sports, such as tennis [5,31–34], throwing and swimming [34–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).

**Figure 1. Top**: morphomap extracts shape information as equiangular semilandmarks from the periosteal (blue) and endosteal surface (orange). **Bottom**: the cross section at 1% of the biomechanical length.

## *2.2. Asymmetry and Cross-Sectional Geometry*

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

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

**Table 1.** *morphomapAsymmetry*: description of the main arguments.


The workflow implemented in *morphomap* is as follows:

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

**Figure 2.** Workflow of the function *morphomapAsymmetry*. In (**A**), one of the two long bones is mirrored, and the two matrices of cortical thicknesses (MM) are computed (**B**). The differences between the two MMs may be computed by calculating (i) the arithmetic difference (**C**), (ii) the percentage change of the target with respect to the reference side (**D**), (iii) the difference between the two MMs and their mean (**E**).

## *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:


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<sup>2</sup> and beta coefficients. R<sup>2</sup> 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.
