Assessments of bilateral asymmetry with application in human skull analysis

Abstract

As a common feature, bilateral symmetry of biological forms is ubiquitous, but in fact rarely exact. In a setting of analytic geometry, bilateral symmetry is defined with respect to a point, line or plane, and the well-known notions of fluctuating asymmetry, directional asymmetry and antisymmetry are recast. A meticulous scheme for asymmetry assessments is proposed and explicit solutions to them are derived. An investigation into observational errors of points representing the geometric structure of an object offers a baseline reference for asymmetry assessment of the object. The proposed assessments are applicable to individual, part or all point pairs at both individual and collective levels. The exact relationship between the developed treatments and the widely used Procrustes method in asymmetry assessment is examined. An application of the proposed assessments to a large collection of human skull data in the form of 3D landmark coordinates finds: (a) asymmetry of most skulls is not fluctuating, but directional if measured about a plane fitted to shared landmarks or side landmarks for balancing; (b) asymmetry becomes completely fluctuating if one side of a skull could be slightly rotated and translated with respect to the other side; (c) female skulls are more asymmetric than male skulls. The methodology developed in this study is rigorous and transparent, and lays an analytical base for investigation of structural symmetries and asymmetries in a wide range of biological and medical applications.

About the Authors:

M. Hou

Roles Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Validation, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliation: Department of Engineering, Medical and Biological Engineering, Faculty of Science & Engineering, University of Hull, Hull, United Kingdom

ORCID logo https://orcid.org/0000-0003-3794-389X

M. J. Fagan

Roles Conceptualization, Formal analysis, Validation, Writing – original draft, Writing – review & editing

Affiliation: Department of Engineering, Medical and Biological Engineering, Faculty of Science & Engineering, University of Hull, Hull, United Kingdom

Introduction

Often the geometric structure of a biological object is represented by a collection of anatomical landmarks divided into three subsets: two of them contain the paired landmarks on two sides of an object respectively, while the third includes the landmarks shared by both sides. The two sides of the object are typically assumed to be geometrically symmetric about a plane in the closest proximity to the shared landmarks. Symmetry can take different forms and abnormal biological development may be reflected in noticeable structural asymmetries. As a fundamental notion, bilateral symmetry has been an important topic in morphometrics, e.g. [1], while other forms of symmetry have been found [2] traceable on the basis of bilateral symmetry. Investigation of bilateral asymmetry is important in biological studies for assessing developmental variability and instability, and also useful in clinical and surgical applications for diagnosing abnormalities and alignment of limbs and organs.

Due to bioecological variations such as genetic mutations and environmental changes together with observational errors caused by uncertainties in landmark placements and measurements, practically no biological object is exactly bilaterally symmetric. The sufficiency of one-side-only analysis of bilaterally symmetric objects was questioned and methods for use or recovery of missing information with two-side data were suggested [3]. In biological studies different ideas of asymmetry assessment have been used to characterise deviations from bilateral symmetry. The three most widely used notions are fluctuating asymmetry (small, random deviations), directional asymmetry (a consistent bias toward a given side) and antisymmetry (a consistent bias toward a random side) [1, 4]. Long standing research in seemingly simple fluctuating asymmetry has been fruitful but the topic remains controversial in particular when used to assess developmental instability [5–9]. The continuing debate demonstrates a lack of certainty about real connections between often subtle variations in the morphological and biological phenomena. Most studies have focused on fluctuating asymmetry as a measure of developmental instability, but many conflicting results may actually be due to inaccuracies in its assessment [4, 8]. Hence, the current study attempts to refine the three asymmetry notions by recasting them in the framework of analytic geometry to facilitate more reliable and accurate asymmetry assessments. The development of a rigorous analytical methodology for asymmetry assessments benefits the evaluation of developmental instability and the diagnosis of diseases causing structural deviations in biological objects.

Traditional approaches [1, 10] use a number of length measurements of the two-side structure of an object to investigate bilateral symmetry. While these measurements capture some morphometric traits in a population, they do not offer a direct assessment of the deviations of individual objects from bilateral symmetry. A more reliable but less effective method appears to be the Euclidean distance matrix analysis (EDMA) [11]. It was adopted in a study of human skull asymmetry [12], where all linear distance differences between two-side landmarks to middle landmarks of a skull were used in the analysis. The full set of linear distances has however a high level of redundancy [13] and only indirectly indicates 3D asymmetry. More systematic studies, e.g. [14, 15], consider morphometric traits collectively, and use a sophisticated combination of the two methods: Procrustes superimposition and statistical analysis in a software tool [16]. The former is about object alignment after reflection of anatomical landmarks on one object side about a prespecified reference, while the latter provides a statistical means of assessing degrees of asymmetry differences between the reflected landmarks and the original ones on the other side. Such a combination has been used in many applications to assess, for instance, human skulls [17] and femurs [18]. Using a mixed method of the Procrustes landmark registration and traditional left-right trait calculations, a negative correlation of bodily fluctuating asymmetry has been established with bodily attractiveness [19].

The current study examines bilateral asymmetry directly in an analytic geometry setting and derives explicit solutions for several types of asymmetry assessments. These types are naturally linked to assessment of the two known categories: object symmetry and matching symmetry [14, 20], meaning symmetry of two inseparable parts, and that of two separate parts of an object, respectively. The proposed asymmetry assessments are applicable to individual landmark pairs and objects, or all of them as a whole in a population. The direct treatment allows evaluation of bilateral asymmetry at both individual and collective levels and determination of its nature (fluctuating asymmetry, directional asymmetry or antisymmetry) in a population. The exact relationship of the proposed asymmetry assessments with the widely used Procrustes method is examined. An application of the developed method in this study to a large collection of human skulls provides some quantitative insights into the nature and extent to which human skulls are bilaterally asymmetric under influence of observational errors of skull landmarks.

Materials and methods

Notations

Set {xi} collects a given number of scalars or (column) vectors. Symbol ′ stands for the transpose of a vector or matrix, and ‘det’ and ‘tr’ abbreviate the determinate and trace of a square matrix respectively. |a| is the absolute value of scalar a, and is the norm of vector or matrix A. The diagonal matrix with n scalars {ai} is denoted by diag(a1, …, an), and becomes the identity matrix, denoted by I, when ai = 1 for all i. In singular value decomposition A = UΛV′, diagonal matrix Λ consists of the singular values, and the columns of matrices U and V contain the left- and right-singular vectors of A respectively. Denoted by a × b is the cross product of vectors a and b. With sgn(0) = 0, sgn(x) = x/|x| for scalar x ≠ 0 is the sign function. Denoted by is a random variable following the normal distribution with mean μx and standard deviation σx, similarly indicates a normal random vector with mean and variance matrix Σx. A+ is the Moore-Penrose inverse of matrix A. Pr(A) is the probability of event A. All quantities under consideration are of real numbers.

Definitions of symmetry

Explicitly or implicitly symmetry always has a reference. Fig 1 illustrates reflection symmetries of a pair of points (p, q) about a point r0, a line (n0, r0) and a plane (n, d), called midpoint, midline and midplane respectively, in a 3D space. Here, n0 is the direction vector of the line passing r0, n the normal to the plane, and d the (signed perpendicular) distance to the plane from the origin of a reference frame which is a 3D Cartesian coordinate system. Both n0 and n are unit vectors.

[Figure omitted. See PDF.]

Fig 1. Point, line and plane reflections between 3D points p and q.

https://doi.org/10.1371/journal.pone.0258146.g001

Respectively on to references r0, (n0, r0) and (n, d), the (orthogonal) projections of point p are r0, p + I0(r0 − p) and p + (d − n′p)n with .

Definition 1 A pair of points (p, q) is said to be symmetric about a reference if q reflects p (or vice versa) through the reference, namely if (1)

Definition 2 Sets {pi} and {qi} with the same number of elements are said to be symmetric about a reference if pair (pi, qi) for all i is so.

Definition 3 Sets {pi} and {qi} with the same number of elements are said to be symmetric about a reference by transformation (a, R, t) (scaling, rotation and translation) on {qi} if pair (pi, aRqi + t) for all i is symmetric about the reference.

To be linked to symmetry notions used in biological applications, the point, line and plane symmetries are associated with spherical, cylindrical and mirror symmetries respectively, while symmetry by transformation generally covers scaling, radial and glide symmetries. The three actions of a transformation do not have to be applied altogether. For instance, (a, I, 0) means no translation or rotation, but scaling, while (1, R, 0) indicates rotation only.

Types of asymmetry assessment

Consider an object having {pi} and {qi} respectively on its two sides, and {ri} shared by both sides. Each of {pi} and {qi} consists of m points, and {ri} k points, with integers m > 0 and k ≥ 0. Hence, an object under consideration has at least one point on each side, and some or possibly no point shared by both sides.

Asymmetry index.

Due to structural deviations and observational errors in an object, inevitably none of the relationships in (1) can hold exactly in real applications. The deviation of (pi, qi) from exact symmetry, called the asymmetry index, is defined by, for i = 1, 2, …, m,(2)When transformation (a, R, t) is applied to qi to have aRqi + t, the corresponding update of ei in (2) is denoted by . Indices {ei} and can be used for asymmetry assessments at individual and collective levels.

Assessment types.

By least squares of asymmetry indices {ei} or , a midpoint r0, midline (n0, r0) or midplane (n, d) as a symmetry reference and a transformation (a, R, t) of points on one side of the reference, say {qi}, can be determined, leading to different types of asymmetry assessments.

1. Type 0 (middle fitting): A reference is determined by best fitting to {ri}.

2. Type 1 (side balancing): A reference is determined by best balancing {pi} and {qi} on its two sides.

3. Type α (mixed types 0 and 1): A reference is determined by simultaneously fitting to {ri} and best balancing {pi} and {qi} on its two sides, with weightings 1 − α and α respectively, for 0 ≤ α ≤ 1.

4. Type 2 (side transforming): A general transformation (a, R, t) of {qi} is determined by making {pi} and {aRqi + t} most symmetric about a given reference.

The least squares criteria for determination of these references and transformations with respect to the four assessment types are(3)where δi stands for the perpendicular distances of point ri to these references:(4)

Clearly, types 0 and 1 are special cases of type α, and in the case of no transformation, namely (a, R, t) = (1, I, 0), type 2 is reduced back to type 0, 1 or α depending on which of the three references is considered. It is possible to further combine types α and 2 to form just one general assessment type by considering the criterion Jβ = (1 − β)Jα + βJ2 with 0 ≤ β ≤ 1. Having such a universal type would be theoretically attractive and effective for algorithm implementation, but is not pursued for sake of clarity and accessibility of the key methodology developed in this study.

Specifications of references and transformations.

In terms of a few data vectors, centroids and matrices defined as(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)references and transformations for assessment types α and 2 are specified as follows with verifications given in S1 Appendix, and types 0 and 1 specified through type α.

The references for type α assessment with 0 ≤ α ≤ 1 are(17)where nb and ns are the eigenvectors corresponding to the biggest and smallest eigenvalues of Xα respectively. Transformations for type 2 assessment are(18)with(19)and (U, V) coming from the singular value decomposition A = UΛV′. The general specification in (18) and (19) remains valid for special (a, R, t), such as (a, I, 0), (1, R, 0) or (1, I, t), etc. This is because a, R and t are considered respectively as independent scalar, matrix and vector in the minimisation.

Assessments of asymmetry

Asymmetry assessments at a collective level are of interest in biological evolutionary studies, while at an individual level, they are useful in medical and clinical applications.

At individual level.

The average and root-mean-square deviation (RMSD) of asymmetry indices {ei} from an object are(20)Pair (pi, qi) of an object is asymmetric if its asymmetry index ei ≠ 0, and the object itself is so if at least one asymmetry index from {ei} is nonzero, namely er ≠ 0. In practical applications, er should be greater than a threshold which takes into account deviations from perfect symmetry caused by observational errors of asymmetry indices. If er exceeds the threshold, ‖e‖ indicates the extent to which the two sides of the object are unbalanced.

At collective level.

For object j in a population, the set of its point pairs is denoted by {(pij, qij)}, and, if any, {rij} the set of points shared by its both sides. The associated asymmetry index eij is calculated from (2) with (pij, qij) replacing (pi, qi). Accordingly, the previously defined terms ri, r0, (n0, r0) and (n, d) now become rij, r0j, (n0j, r0j), and (nj, dj) respectively. If symmetry by transformation is considered, qij should be replaced by aj Rj qij + tj for all i and each j.

Point pairs {(pi, qi)} can be considered random with {(pij, qij)} being their samples drawn from a population. Similarly, their asymmetry indices {ei} defined in (2) have samples {eij}. Moreover, object asymmetry index e defined in (20) has the averages of {eij} over i as its sample. Asymmetries of point pairs and objects in a population are defined below, followed by their assessments.

Asymmetry of point pairs: Pair (pi, qi) in a population is said to be asymmetric if its asymmetry index ei is statistically nonzero, namely Pr(‖ei‖ = 0) = 0, and the asymmetry is

1. a). fluctuating if with ;

2. b). directional if with ;

3. c). antisymmetric if with and a mixing index 0 < γ < 1;

where , and are the vector-valued means, Σi, Σai and Σbi the variance matrices.

Asymmetry of objects: Objects in a population are said to be asymmetric if object asymmetry index e is statistically nonzero, namely Pr(‖e‖ = 0) = 0, and the asymmetry is

1. a). fluctuating if with ;

2. b). directional if with ;

3. c). antisymmetric if with and a mixing index 0 < γ < 1;

where , and are the vector-valued means, Σ, Σa and Σb the variance matrices.

Assessments: Assessing asymmetries of point pairs and objects in a population is practically carried out by assessing the same in its samples. In respect to a sample of size l, denote the average and RMSD of asymmetry indices {eij} over j by(21)and those over i by(22)While asymmetry index ei of (pi, qi) in (2) has {eij} with varying j as its sample, indices e and er in (20) have {e⋅j} and {erj} in (22) as their samples respectively.

For assessing asymmetries of point pairs, index eir ≠ 0 approximates zero probability of ‖ei‖ = 0 in the population, sample mean ei. approximates population mean , and the sampling distribution of {eij} approximates the population distribution of ei.

For assessing asymmetry of objects, indices erj ≠ 0 for all j in a sample approximate zero probability of ‖e‖ = 0 in the population, sample mean e⋅j approximates population mean , and the sampling distribution of {e⋅j} approximates the population distribution of e.

In practice, distributions of ei and e can be graphically examined with 3D scatter plots of their samples {eij} and {e⋅j}, and/or histograms for individual components of the samples. Moreover, criteria eir ≠ 0 and erj ≠ 0 should be replaced by them exceeding some thresholds counting for observational errors of asymmetry indices.

Internal referencing for asymmetry index alignment.

Asymmetry indices {eij} are dependent on objects’ orientations. The dependence introduces artefacts to the sampling distribution of {eij} since unlikely these objects have been aligned with each other when measurements of them were being taken. Removal of the artefacts can be achieved by referring sample indices {eij} to an arbitrary internal frame uniformly specified for every object. This amounts to using instead of {eij} in asymmetry assessments, where is a rotation matrix determined by points of object j in the same way for all objects.

A choice is with its direction vectors(23)with and . In (23), nj for object j is the normal to the midplane determined in the same way as n in (17), the sign function makes point to the side of {pij}, and is arbitrary as long as it is not parallel to , namely . For example, as inspired by rα in (6),(24)could be such a choice, where scalar γ and integers kt and mt are fixed for all objects with 0 ≤ γ ≤ 1, 0 ≤ kt ≤ k and 0 ≤ mt ≤ m.

Removal of the artefacts caused by external referencing from asymmetry indices can also be achieved by aligning the objects to each other through pre-processing the coordinate data, for instance, using the Procrustes superimposition. Clearly, in (21) and (22), RMSDs {eir} and {erj} of asymmetry indices {eij} are coordinate independent. Hence, if only asymmetry but not its nature is concerned in an application, the removal or alignment is not needed.

Observational errors.

Uncertainties in the placement and measurement of points, and inaccuracy in numerical calculation of asymmetry indices contribute to the observational errors in asymmetry analysis. To separate effects of structural deviation ϵ from those of observational error ε on asymmetry assessments, it is useful to express the points as(25)with(26)where, for simplicity, object index j is dropped. As the average of itself and the mirrored copy of its correspondence, and are perfectly symmetric about the reference on to which ri, if any, has projection .

It is reasonable to assume observational errors as noise vectors satisfying Gaussian distribution: with σn being the standard deviation for each of the three components in a noise vector. If structural deviations are set to be zero (ϵp,i = ϵq,i = ϵr,i = 0), simulations of the simple data model (25) and (26) can generate useful information for determination of thresholds used to distinguish between effects of structural deviations and those of observational errors on asymmetry assessments.

Flowcharts of asymmetry assessments

To facilitate use of the developed scheme, three explanatory flowcharts are shown in Figs 2 and 3. The procedures in Fig 2 are for plane asymmetries at individual levels. If point or line asymmetry is considered, according to (2), (17) to (19), few straightforward adjustments are needed in the fourth block of type α assessment procedure, and in the first and fifth blocks of type 2 assessment procedure.

[Figure omitted. See PDF.]

Fig 2. Flowcharts of plane asymmetry assessments of types α and 2 at individual levels.

https://doi.org/10.1371/journal.pone.0258146.g002

[Figure omitted. See PDF.]

Fig 3. Flowchart of asymmetry assessments at collective levels.

https://doi.org/10.1371/journal.pone.0258146.g003

Two testing datasets

Human skull dataset.

A subset consisting of 359 skulls was selected from the collection of 889 human skulls in [21]. It includes skulls of African Americans (69 females and 71 males), and European Americans (100 females and 119 males). The only criterion for inclusion of a particular skull in this current study was the completeness of 10 pairs of side landmarks and 8 middle landmarks, as summarised in Table 1 and indicated in Fig 4. There was however one exception, namely removal of sample ID876 (European American female) due to an obvious measurement error of its left asterion. 3D landmark coordinates of the skulls in the subset obtained from the collection in [21] are used in the asymmetry assessments. To determine midplanes, the eight middle landmarks are used in type 0 assessments, while the ten side landmark pairs are used in type 1 assessments, and all of them are used in type α assessments. Once a midplane has been determined for asymmetry referencing, only the ten side landmark pairs need consideration and hence are used in the subsequent analysis of skull asymmetries.

[Figure omitted. See PDF.]

Fig 4. Landmarks indicated in different views of a skull.

https://doi.org/10.1371/journal.pone.0258146.g004

[Figure omitted. See PDF.]

Table 1. 2 × 10 side landmarks and 8 middle landmarks.

https://doi.org/10.1371/journal.pone.0258146.t001

Simulated data.

This dataset contains simulated skulls whose landmarks are obtained by adding normally distributed independent random values of zero mean and a given standard deviation to each of the three landmark coordinates of the perfectly symmetric skulls. The number of the simulated skulls is the same as that of real skulls. The identical landmarks of the real skulls have been used to produce the perturbed symmetric skulls in the simulated data. A perfectly symmetric skull is represented by the landmarks described in (26). A standard deviation of 1mm for landmark displacement errors is chosen to count mainly for inaccuracy of several typical measurement methods, ranging from physical measurements by a digital calliper [22] or MicroScribe [23], to digital measurements by cone beam computed tomography [24]. By the three-sigma rule of thumb, this standard deviation corresponds to a maximal error of 3mm in landmark placement and measurement. Asymmetry indices calculated from the simulated dataset will be used as baselines for comparison with those generated from real skull dataset in asymmetry assessments.

Results

The theoretical development of asymmetry assessments is applied to the skull datasets as a case study. In view of the nature of skull asymmetry, only plane asymmetry is examined. In type α assessment, α = 0.5 is used, and in type 2 assessment, the midplane is the one determined in type 1 assessment. For alignment of asymmetry indices, parameters γ = 1/2, kt = k/2 and mt = m/2 are used in (24). The abbreviations for categorisation are AA: African American, EA: European American, F: Females, M: Males.

In some plots, where asymmetries of 359 individual skulls are concerned, the horizontal axes represent skull indices up to that total number. In other plots, when comparisons among the four skull categories are considered, the skull indices shown for the horizontal axes are only up to 119 which is for EA male skulls, the largest number of the skulls in the four categories.

Noise effects on symmetry

With respect to different assessment types, landmark-pair asymmetry of perturbed symmetric skulls caused by observational errors are produced. As expected, at the population level, asymmetry of every landmark pair of perturbed symmetric skulls is fluctuating. S1 Fig in S2 Appendix indicates the landmark-pair asymmetry of perturbed symmetric skulls caused by observational errors. As an example of type 1 assessment, scattered points and histograms of asymmetry index of the 6th landmark pair are shown in S2 Fig in S2 Appendix.

With respect to different assessment types, the first sub-figure in S3 Fig in S2 Appendix indicates indistinguishable variations in the RMSDs of asymmetry indices for the perturbed symmetric skulls. Scatter plots in S3 Fig in S2 Appendix show distribution ranges of skull asymmetry indices {e⋅j} with respect to different assessment types. These responses from the perturbed symmetric skulls provide comparing baselines for asymmetry assessments of the real skulls. As will be seen, in general, the baseline asymmetry indices are minor compared with large magnitudes of asymmetry indices from the real skulls.

Landmark-pair asymmetry

The RMSDs in Fig 5 suggest overall asymmetries in most of the 10 landmark pairs, induced by bioecological factors in the studied human skulls. In particular, landmark pairs of female skulls appear more asymmetric than those of male skulls. Also, as shown in S4 Fig in S2 Appendix, asymmetry indices of landmark pairs of female skulls have more variations than those of male skulls do.

[Figure omitted. See PDF.]

Fig 5. Asymmetry of landmark pairs.

The plots are for the perturbed symmetric skulls (model) and original skulls with respect to different assessment types (0, 0.5, 1 and 2) and in different categories (AAF, AAM, EAF and EAM), where in brackets are the means and standard deviations of the RMSDs for asymmetry indices.

https://doi.org/10.1371/journal.pone.0258146.g005

Fig 6 indicates that asymmetry of landmark pairs is not fluctuating. As an example, for type 0 assessment, scattered points and component histograms of asymmetry indices of the 6th landmark pair are shown in Fig 7, where the components of the asymmetry indices appear to follow slightly skewed normal distributions with clearly non-zero means of x and z components. Along with similar plots of scattered points and histograms for types 1 and 2 assessments shown in S5 and S6 Figs in S2 Appendix, these plots indicate that the 6th landmark pair is directionally asymmetric.

[Figure omitted. See PDF.]

Fig 6. Non-fluctuating asymmetry of landmark pairs.

The plots are for the perturbed symmetric skulls (model) and original skulls with respect to different assessment types (0, 0.5, 1 and 2) and in different categories (AAF, AAM, EAF and EAM), where in brackets are the means and standard deviations of the sample-mean norms of asymmetry indices.

https://doi.org/10.1371/journal.pone.0258146.g006

[Figure omitted. See PDF.]

Fig 7. Type 0 assessment: Directional asymmetry of landmark pairs.

Scattered points and component histograms of asymmetry indices {e6j} in type 0 assessment from the 359 perturbed symmetric skull and all real skulls are shown, where in brackets are the means and standard deviations of components of {e6j}.

https://doi.org/10.1371/journal.pone.0258146.g007

Skull asymmetry

Fig 8 shows that the majority of the skulls are asymmetric and on average female skulls are more asymmetric than male skulls. In type 0 assessment, Fig 9 shows that the skull asymmetry is not fluctuating, but directional due to clear deviations of the mean of {e⋅j} from that of perturbed symmetric skulls. In type 1 assessment, Fig 10 shows that the asymmetry remains directional since in y-direction the mean of {e⋅j} is much greater than that of perturbed symmetric skulls although the x and z components of the mean are close to those of perturbed symmetric skulls. Interestingly, in type 2 assessment shown in Fig 11, the skull asymmetry becomes completely fluctuating since the scale of the scatter plot (first plot) remains at 10−13 mm as in the third scatter plot of S3 Fig in S2 Appendix, and the amounts of the associated rotations and translations are very small in general, with the means at 2° and 0.7 mm respectively. There are however 3 or 4 outliers (about 1% of the total skulls) requiring much larger rotational and/or translational adjustments. The angular norm in Fig 11 is calculated from the three rotation angles as a vector associated with each rotation matrix. These angles are converted from the quaternions which are first calculated from each rotation matrix.

[Figure omitted. See PDF.]

Fig 8. Skull asymmetry in different categories.

RMSDs of asymmetry indices of the first 119 perturbed symmetric skulls (model) and original skulls with respect to different assessment types (0, 0.5, 1 and 2) and in different categories (AAF, AAM, EAF and EAM) are shown, where in brackets are the means and standard deviations of the RMSDs.

https://doi.org/10.1371/journal.pone.0258146.g008

[Figure omitted. See PDF.]

Fig 9. Type 0 assessment: Directional asymmetry of skulls.

Scattered points and histograms of asymmetry indices of the 359 perturbed symmetric skull and all real skulls are shown, where in brackets are the the means and the standard deviations of components of {e⋅j}.

https://doi.org/10.1371/journal.pone.0258146.g009

[Figure omitted. See PDF.]

Fig 10. Type 1 assessment: Directional asymmetry of skulls.

Scattered points and histograms of asymmetry indices of the 359 perturbed symmetric skull and all real skulls are shown, where in brackets are the the means and standard deviations of components of {e⋅j}.

https://doi.org/10.1371/journal.pone.0258146.g010

[Figure omitted. See PDF.]

Fig 11. Type 2 assessment: Fluctuating asymmetry of skulls.

Scattered points of averaged asymmetry indices and norms of rotation angles and translation displacements of the 359 perturbed symmetric skull and all real skulls are shown, where in brackets are the means and standard deviations of the norms.

https://doi.org/10.1371/journal.pone.0258146.g011

Procrustes analysis

As detailed in the Discussion section, few straightforward adjustments of the specification of type 2 assessment in (18) and (19) give rise to explicit solutions needed in the Procrustes analysis. Use of these explicit and direct transformations easily avoids undesirable rotations and reflections in the orthogonal Procrustes transformation and superimposition used for object reflection and alignment respectively.

For object symmetry assessment, midplane (nj, dj) for referencing object j is the one used for type 0 assessment, and the superimposition is applied to and {pij}, where represents the reflection of qij through the jth midplane. For matching symmetry assessment, reflection Rj and translation tj are determined for object j first and the superimposition is then applied to {Rjqij + tj} and {pij}. Accordingly, after the reflections and superimpositions, assessment indices and are used for object and matching asymmetries respectively. In the superimposition, all objects are aligned to the first skull, namely the ordinary Procrustes analysis is used.

Although the procedure described above may not be exactly the same as those adopted by other investigators in applications of the Procrustes method to asymmetry assessments, it effectively captures the essence of the method and links it to the scheme developed in the current study.

Fig 12 shows a brief comparison between type 0 and 2 assessments and the Procrustes analysis of asymmetries of landmark pairs and objects for the perturbed symmetric skulls and original skulls. Two observations and some comments are made below:

1. a). In asymmetry assessment of landmark pairs, both methods work well for producing the RMSDs of asymmetry indices for real skulls distinguishable from those for perturbed symmetric skulls. The difference between real and model skull RMSDs is greater in the current method than in the Procrustes method. This also true after scaling (not shown).

2. b). In asymmetry assessment of skulls, both methods produced the identical RMSD of asymmetry indices for real skulls, and comparable ones for perturbed symmetric skulls.

[Figure omitted. See PDF.]

Fig 12. A comparison with the Procrustes analysis.

The plots are for the perturbed symmetric skulls (model) and original skulls (real) with respect to type 0 assessment (top) and type 2 assessment (bottom) in this work, and object symmetry (top) and matching symmetry (bottom) of Procrustes method; RMSDs of asymmetry indices of landmark pairs (left), and objects (right).

https://doi.org/10.1371/journal.pone.0258146.g012

The adjustments to the Procrustes analysis as explained above avoid making it unfairly disadvantaged in such a brief comparison. This leads to equivalence of the two methods when used for assessing asymmetries of objects, as confirmed by observation b). Although observation a) indicates validity of both methods for assessing asymmetries of point pairs of these objects, it does suggest an advantage of the current method over the other. The reason behind this lies in different alignments of asymmetry indices by the two methods. Since object alignment plays no role in asymmetry assessment of individual skulls, it is logic to have observation b). Recall that in Figs 9–11, examination of distributions of skull asymmetry indices {e⋅j} requires removal of the artefacts induced by unknown orientations of the objects, alignment of these indices for all skulls in the sample has been used as well according to (23) and (24).

As a final remark, the Procrustes method could be improved by using the generalised Procrustes superimposition which is an iterative procedure. A detailed study of that is beyond the scope of the current work.

Discussion

General consideration

In an analytic geometry setting, this study has rigorously defined several asymmetry assessments and derived explicit solutions to them. Recast in the framework of analytic geometry, notions of fluctuating asymmetry, directional asymmetry and antisymmetry have also been refined. The main computations used for the asymmetry assessments are calculations of the eigenvalue/eigenvector and singular value decomposition of data matrices. Most popular software packages for numerical calculations can be used to complete these computations. The methods have been generally developed for 3D objects, but are applicable with minor adjustments to 2D objects, where point and/or line asymmetry is concerned. An extension of the treatments to nD objects for n > 3 is also possible, which is of theoretical interest at least.

The illustration of the assessments in a study of human skull bilateral asymmetry shows the efficiency of the methods. Effects of uncertainties in the landmark placement and measurement on asymmetry assessment have been examined, which helps differentiate the effects of observational errors and bioecological factors. The current study assumes that the correspondences between points in the two halves of an object are known. If this is not the case, permutations of points in one half of the object need to be combined with asymmetry assessments. However, the developed treatments do not apply to an object with different numbers of landmarks on its two sides and in particular when no clear correspondence between any subsets of the points on its two has been established.

Assessment types

Different assessment types evaluate bilateral asymmetry of an object from different perspectives. Type 0 assessment for evaluation of deviations from point, line or plane symmetry is useful when there are middle points centred around a midpoint, midline or midplane. Type 1 assessment applies to the situations where no middle points are recorded, or they are not well centred and better interpreted as being scattered around a midpoint, midline or midplane. As a trade-off between types 0 and 1, type α could be more convenient and useful for asymmetry assessment in practical applications since types 0 and 1 are its special cases. Type 2 assessments are normally used for examining asymmetries of objects with two disjoint sides.

If an object is examined sequentially by all three assessment types 0, 1 and 2, asymmetry indices will be largest for type 0 assessments and smallest for type 2 assessments, and the asymmetry index associated with type α assessment is normally between those of types 0 and 1. While the solution of fitting a line or plane to a set of points used in type 0 assessment is known and particularly implied in a general solution to linear manifold fitting [25], a simple direct derivation of type α and 2 assessments for point, line and plane asymmetries can be found in S1 Appendix. With respect to the categories of object symmetry and matching symmetry [14, 20], type α belongs to the former, while type 2 to the latter.

Relationship with Procrustes method

In [14], notions of the object symmetry and matching symmetry in terms of reflection and permutational transformations on a data matrix were formulated, and asymmetry assessments were based on statistical and geometrical analyses of the Procrustes shape manifold. The Procrustes method is normally used to assess object and matching symmetries [20]. Originally developed for studying similarity of objects, the well-known Procrustes method, see e.g. [26, 27], allows rotation or reflection along with scaling and translation to be applied to points on one side to be compared with those on the other. Conceptually symmetry is secondary to similarity, which makes a general method developed for one not best suit the other.

In the current study, through direct formulations in 3D analytic geometry, the object symmetry is assessed by type α with 0 ≤ α ≤ 1, and the matching symmetry by type 2 assessments of deviations from point, line or plane symmetry. An advantage of the current treatment is its direct formulation with explicit solutions. As shown in S1 Appendix, type 2 assessment is readily reducible to a minor variation of the standard orthogonal Procrustes and Wahba’s problems [28–30]. Also, as explained below, with few straightforward adjustments, the transformation specified for type 2 assessment becomes transformations that best meet the needs of asymmetry analysis using the Procrustes method.

In general, with or without translation, rotation of one of two identical halves does not make a symmetric pair, and conversely, reflection of half of a symmetric pair breaks the symmetry. Hence, in asymmetry analysis, only reflections with translations in the Procrustes transformations should be kept, which can be ensured by simply setting a = 1, gi = pi and using slightly adjusted R = Udiag(1, 1, −detUV′)V′ in the solution (a, R, t) of type 2 assessment given in (18) and (19). The adjustments amount to re-defining type 2 assessment as to find a reflection along with a translation of points on one object side so that the transformed points become closest to those on the other. For the purpose of asymmetry assessment, this re-definition is mathematically equivalent to the original one. Similarly, for aligning two objects, only rotations with translations in the Procrustes superimposition should be retained, which can be ensured by the transformation (a, R, t) of type 2 assessment if all points of one object are re-denoted by {gi}, and those of another by {qi} in (18) and (19).

Asymmetry of human skulls

Asymmetry of mammalian skulls is well known, see, e.g. [31, 32]. Studies of asymmetry of human skulls have initially used length and angle measurements and lately the Procrustes analysis of landmarks, see, e.g. [12, 17, 33–35]. The current study shows that with respect to a midplane fitted to the shared landmarks (type 0 assessment) or to side landmarks for balancing the two sides (type 1 assessment), asymmetry of the studied human skulls is generally directional, which is consistent with the previous discoveries in the literature (e.g. [33, 34]). In general, in the studied population, female skulls appear clearly more asymmetric than male skulls. Nevertheless, no significant cranium directional asymmetry differences were found in [35] between the sexes in a modern Greek population. Finally, the current study shows that, after a gentle coordinate transformation (type 2 assessment) is applied to the landmarks of one side, the skull asymmetry becomes extraordinarily fluctuating, and the associated rotational angles and translational displacement are minor. This interesting feature does not seem to have been reported elsewhere.

Other applications

Apart from analysis of rotational and translational symmetries for discovery of the nature of bilateral symmetry and exploration of scaling, radial and glide symmetries, type 2 assessment of matching symmetry is useful in other applications. For instance, forensic and archaeological investigations are often interested in matching one part of an object to another from a number of candidates. In medical and clinical applications, a large deviation from symmetry of two sides of an organ can be a diagnostic factor for certain diseases affecting that organ. Alignment of one object with another to make a symmetric pair is likely required in various surgical operations. Investigations on asymmetry of species in evolutionary biological studies at macroscopic and microscopic levels benefit research on developmental instability and its causes in a population. The number of publications on asymmetry assessments in these diverse disciplines is huge, and specific referencing is spared here. The current study has attempted to provide investigators in wide application areas of asymmetry assessments with a solid and transparent analytical base.

Supporting information

S1 Appendix. Derivation of Eqs (17) to (19) [36].

https://doi.org/10.1371/journal.pone.0258146.s001

(PDF)

S2 Appendix. Additional simulation results.

https://doi.org/10.1371/journal.pone.0258146.s002

(PDF)

Citation: Hou M, Fagan MJ (2021) Assessments of bilateral asymmetry with application in human skull analysis. PLoS ONE 16(10): e0258146. https://doi.org/10.1371/journal.pone.0258146

References

1. Palmer A.R. and Strobeck C., ‘Fluctuating asymmetry: Measurement, analysis, patterns’, Annu Rev Ecol Syst, 1986 Jan;17(1):391–421.

2. Savriama Y. and Klingenberg C.P., ‘Beyond bilateral symmetry: geometric morphometric methods for any type of symmetry’, BMC Evol Biol, 2011 Sep;11:280. pmid:21958045

3. Cardini A., ‘Lost in the other half: Improving accuracy in geometric morphometric analyses of one side of bilaterally symmetric structures’, Syst Biol, 2016 Nov;65(6):1096–1106. pmid:27288476

4. Loehr J., Leinonen T., Herczeg G., O’Hara R. and Merilä J, ‘Heritability of asymmetry and lateral plate number in the three spine stickleback’, PLoS ONE, 2012 Jul;7(7):e39843. pmid:22792190

5. Brakefield P.M. and Breuker C.J., ‘The genetical basis of fluctuating asymmetry for developmentally integrated traits in a butterfly eyespot pattern’, Proc R Soc B, 1996 Nov;263(1376):1557–1563.

6. Lens L., Van Dongen S., Kark S. and Matthysen E., ‘Fluctuating asymmetry as an indicator of fitness: can we bridge the gap between studies?’, Biol Rev, 2002 Feb;77(1):27–38. pmid:11911372

7. Markow T.A., ‘Evolutionary ecology and developmental instability’, Annu Rev Entomol, 1995 Jan;40(1):105–120.

8. van Dongen S., ‘Fluctuating asymmetry and developmental instability in evolutionary biology: past, present and future’, J Evol Biol, 2006 Nov;19(6):1727–1743.

9. van Valen L., A study of fluctuating asymmetry. Evolution, 1962 Jun;16(2):125–142.

10. Graham J.H. and Özener B., ‘Fluctuating asymmetry of human populations: A review’, Symmetry, 2016 Dec;8(12):154.

11. Lele S. and Richtsmeier J.T., ‘Euclidean distance matrix analysis: a coordinate-free approach for comparing biological shapes using landmark data’, Am J Phys Anthropol, 1991 Nov;86(3):415–427. pmid:1746646

12. DeLeon V.B., ‘Fluctuating asymmetry and stress in a medieval Nubian population’, Am J Phys Anthropol, 2007 Apr;132(4):520–534. pmid:17243154

13. Dokmanić I., Parhizkar R., Ranieri J. and Vetterli M., ‘Euclidean Distance Matrices: Essential theory, algorithms, and applications’, IEEE Signal Process Mag, 2015 Nov;32(12):12–30.

14. Mardia K.V., Bookstein F.L. and Moreton I.J., ‘Statistical assessment of bilateral symmetry of shapes’, Biometrika, 2000 Jun;87(2):285–300.

15. Klingenberg C.P., Barluenga M. and Meyer A., ‘Shape analysis of symmetric structures: quantifying variation among individuals and asymmetry’, Evolution, 2002 Oct;56(10):1909–1920. pmid:12449478

16. Adams D.C. and Otárola-Castillo E., ‘geomorph: an R package for the collection and analysis of geometric morphometric shape data’, Methods Ecol Evol, 2013 Apr;4(4):393–399.

17. Bigoni L., Krajíček V., Sládek V., Velemínský P. and Velemínská J., ‘Skull shape asymmetry and the socioeconomic structure of an early medieval central European society’, Am J Phys Anthropol, 2013 Mar;150(3):349–364. pmid:23283725

18. Mopin C., Chaumoître1 C., Signoli M. and Adalian P., ‘Developmental stability and environmental stress: A geometric morphometrics analysis of asymmetry in the human femur’, Am J Phys Anthropol, 2018 Sep;167(1):144–160. pmid:30175505

19. Brown W.M, Price M.E., Kang J., Pound N., Zhao Y. and Yu H., ‘Fluctuating asymmetry and preferences for sex-typical bodily characteristics’, PNAS, 2008 Set;105(35):12938–12943. pmid:18711125

20. Klingenberg C.P., ‘Analyzing fluctuating asymmetry with geometric morphometrics: Concepts, methods, and applications’, Symmetry, 2015 Jun;7(2):843–934.

21. Slice D.E. and Ross A.,’Geometric morphometric classification of crania for forensic scientists’, 3D-ID manual, version 1.0. 2010. https://www.3d-id.org/forensic-3d-coordinates.

22. Hassan B., van der Steltand P. and Sanderink G., ‘Accuracy of three-dimensional measurements obtained from cone beam computed tomography surface-rendered images for cephalometric analysis: influence of patient scanning position’, Eur J Orthod, 2009 Apr;31(2):129–134. pmid:19106265

23. Pandya A., Siadat M.-R. and Auner G, ‘Design, implementation and accuracy of a prototype for medical augmented reality’, Comput Aided Surg, 2005 Jan;10(1):23–35. pmid:16199379

24. Periago D.R., Scarfe W.C., Moshiri M., Scheetz J.P., Silveira A.M. and Farman A.G., ‘Linear accuracy and reliability of cone beam CT derived 3-dimensional images constructed using an orthodontic volumetric rendering program’, Angle Orthod, 2008 May;78(3):387–395. pmid:18416632

25. Späth H., ‘Orthogonal least squares fitting with linear manifolds’, Numer Math, 1986 Jul;48:441–445.

26. Goodall C., ‘Procrustes methods in the statistical analysis of shape’, J R Statist Soc B, 1991;53(2):285–339.

27. Schönemann P.H. and Carroll R.M., ‘Fitting one matrix to another under choice of a central dilation and a rigid motion’, Psychometrika, 1970 Jun;35(2):245–255.

28. Farrell J.L., Stuelpnagel J.C., Wessner R.H., Velman J.R. and Brook J.E., ‘A least squares estimate of spacecraft attitude (Grace Wahba)’, SIAM Rev, 1966 Jul;8(3):384–386.

29. Schönemann P.H., ‘A generalized solution of the orthogonal Procrustes problem’, Psychometrika, 1966 Mar;31(1), 1–10.

30. Eckart C. and Young G., ‘The approximation of one matrix by another of lower rank’, Psychometrica, 1936 Sep;1(3):211–218.

31. del Castillo D.L., Segura V., Flores D.A. and Cappozzo H.L., ‘Cranial development and directional asymmetry in Commerson’s dolphin, Cephalorhynchus commersonii commersonii: 3D geometric morphometric approach,’ J Mammal, 2016 Sep;7(5):1345–1354.

32. Howell A.B., ‘Asymmetry in the skulls of mammals’, Proc US Natl Mus, 1925;67(2599):1–18.

33. Chovalopoulou M., Papageorgopoulou C. and Bertsatos A., ‘Cranium asymmetry in a modern Greek population sample of known age and sex’, Int J Legal Med, 2017 May;131(3):803–812. pmid:27485096

34. LeMay M., ‘Asymmetries of the skull and handedness: Phrenology revisited’, J Neurol Sci, 1977 Jun;32(2):243–253.

35. Woo T.L., ‘On the asymmetry of the human skull’, Biometrika, 1931 May;22(3-4):324–352.

36. Wahba G, ‘A least squares estimate of spacecraft attitude’, SIAM Rev, 1965 Jul;7(3):409.