Analysis of In Vivo Skin Anisotropy Using Elastic Wave Measurements and Bayesian Modelling

In vivo skin exhibits viscoelastic, hyper-elastic and non-linear characteristics. It is under a constant state of non-equibiaxial tension in its natural configuration and is reinforced with oriented collagen fibers, which gives rise to anisotropic behaviour. Understanding the complex mechanical behaviour of skin has relevance across many sectors including pharmaceuticals, cosmetics and surgery. However, there is a dearth of quality data characterizing the anisotropy of human skin in vivo. The data available in the literature is usually confined to limited population groups and/or limited angular resolution. Here, we used the speed of elastic waves travelling through the skin to obtain measurements from 78 volunteers ranging in age from 3 to 93 years old. Using a Bayesian framework allowed us to analyse the effect that age, gender and level of skin tension have on the skin anisotropy and stiffness. First, we propose a new measurement of anisotropy based on the eccentricity of angular data and conclude that it is a more robust measurement when compared to the classic “anisotropic ratio”. Our analysis then concluded that in vivo skin anisotropy increases logarithmically with age, while the skin stiffness increases linearly along the direction of Langer Lines. We also concluded that the gender does not significantly affect the level of skin anisotropy, but it does affect the overall stiffness, with males having stiffer skin on average. Finally, we found that the level of skin tension significantly affects both the anisotropy and stiffness measurements employed here. This indicates that elastic wave measurements may have promising applications in the determination of in vivo skin tension. In contrast to earlier studies, these results represent a comprehensive assessment of the variation of skin anisotropy with age and gender using a sizeable dataset and robust modern statistical analysis. This data has implications for the planning of surgical procedures and questions the adoption of universal cosmetic surgery practices for very young or elderly patients.


Introduction
The skin is a vital organ for a range of bodily functions including protection from the environment and temperature regulation [1]. It is constantly under varying amounts of tension and must be able to withstand significant flexion and deformation for daily tasks like locomotion. Understanding the mechanical properties of the skin is important for many different applications and industries: in the cosmetic industry, products must be assessed in terms of emolliency and hydration of the skin; in the design of anthropomorphic devices like crash test dummies and surgical simulators [1], the mechanical behaviour of the skin must be replicated as closely as possible; and in a surgical setting, a thorough understanding of the skin's mechanical properties is essential. For example, understanding skin growth through tissue expansion is necessary for breast reconstruction and burn victims [2].
Previous publications have examined many different mechanical properties of skin including viscoelasticity [3], the nonlinear stress-strain relationship [4,5], failure properties [6,7] and anisotropy [8,9]. The experimental methods employed have included extension [8], suction [10], torsion [11], indentation [12] and expansion [2], amongst others. In this paper, we focus mainly on the anisotropic nature of the in vivo skin tension which was first noted in the 19th century by anatomist Karl Langer [13]. Preferred lines of tension have become known as "Langer Lines" or skin tension lines and are used by surgeons to select the optimum orientation of skin incisions so as to reduce scarring [14]. Identification of these patient-specific lines is non-trivial and surgeons must rely on generic maps or an imprecise pinch test that requires significant experience to interpret [15][16][17]. Recent research has shown that minimising the skin tension across wounds is the single most important factor in scar prevention that is within a surgeon's control [18,19]. To that end, quantitative knowledge of both the direction and the anisotropic nature of skin tension lines is an essential component of wound closure. A deeper understanding of how they vary across a population may provide further optimisation of closure techniques, particularly for elderly or very young patient cohorts.
More recently, researchers have sought to develop techniques to determine the in vivo orientation of skin tension lines objectively: these include those using suction based devices [20], in vivo extensometry [21,22] and elastic wave propagation [15,3]. These papers have shown that these techniques can successfully identify skin tension lines, and that their orientation is patient-specific, depending on many different factors including location, age, health, BMI, ethnicity and hydration [15]. However, with the exception of Ruvolo et al. [3], none of these papers have comprehensively considered the level of anisotropy of these skin tension lines in vivo, i.e. how the tension levels in two orthogonal directions differ, and how that aspect varies by age and gender. In Ruvolo et al., 239 volunteers ranging in age from newborn to 75 years old underwent testing using elastic wave propagation. Ruvolo et al. noted that previous studies had found only a weak dependence on elastic wave speed with age [23][24][25]; however, these studies all employed poor resolution angular data (measurements taken every 45°) and, having undersampled, they may have missed important information. While great care was taken by Ruvolo et al. to overcome this issue by sampling every 3°, their anisotropy results are reported in terms of the classic "anisotropic ratio", which is a simple ratio between the fastest and slowest wave speed. In the current study we propose, instead, to report the eccentricity of an ellipse fit to the circular data, which may offer a more representative and robust measure of the in vivo anisotropy. Additionally, the current study includes a sizeable dataset with a large range of ages (78 individuals, age 3-93) in contrast to Laiacona (19 subjects, age [18][19][20][21][22][23][24][25][26][27][28][29][30] [20], Boyer (20 subjects, age 20-65) [21] and finally, Hermanns (110 subjects, age 19-93) [25], who did not include infants.
Finally, previous papers have all employed hypothesis testing to support their conclusions. It is now commonly accepted that there are significant issues with the use of p-values in scientific research [26]. Bayesian methods for data analysis provide a principled framework for inference, uncertainty assessment and inclusion of prior information [27]. These methods are flexible, capable of handling complex correlation structures and can eliminate the need for p-values and Null Hypothesis Significance Testing (NHST) [28][29][30].
The objective of this paper is to determine the level of in vivo skin anisotropy and determine how it varies with age and gender. Here we used elastic wave propagation to determine the speed of surface waves traveling through the skin of 78 subjects (age 3-93). Bayesian statistical methods were then employed to examine the significance and effects of age and gender. Furthermore, we examine how skin anisotropy is affected by skin tension and discuss its implications for surgical practice.

Data Collection
The Reviscometer ® (Model RVM 600, Courage & Khazaka Electronic GmbH) is used to examine the mechanical properties of the skin. The device consists of a handheld probe connected to a central controller and a laptop (see Fig. 1a).
The tip of the probe contains two piezoelectric transducers that are 2 mm apart. One transducer emits a Rayleigh surface wave in the form of an acoustic pulse on the skin surface, the other detects the resulting wave and records the time taken for the wave to propagate across the surface of the skin, in one orientation. A hollow plastic fixture also facilitates precise measurements at 10° increments (see Fig. 1b) allowing us to see how the surface wave speed varies for different angles.
By default, the measurement is in arbitrary units called "Resonance Running Time" (RRT). The device was calibrated by assuming the wave speed follows that of a Rayleigh wave travelling on an unstressed, incompressible, linear elastic, isotropic material [31]. The wave speed is then related to the stiffness through: where E is the Young modulus, is the material density and v is the wave speed [32]. Specifically, the Young moduli of 3 elastomers (Techsil 25 Silicone, Polyurethane and MVQ Elastomer) were determined with tensile tests and the average RRT (3 tests) determined for each sample. The conversion for each material is detailed in Table 1. The average RRT was found to be 0.284 s.
Ethical approval for the study was granted by the Human Research Ethics Commitee at University College Dublin (25-18-75). A total of 78 subjects were tested with 37 female and 41 male volunteers, aged between 3 and 92 years of age, see Table 2.
Measurements were obtained on either the left or right volar forearm approximately 5 cm proximal to the wrist, see Fig. 1b. This site was chosen as a convenient, flat surface with minimal body hair, veins or bone.
For each volunteer, two configurations were explored: the "natural configuration", where measurements were carried out on the skin with no interference and the "stretched configuration", where an additional stretch was applied to the skin in the direction of the fastest traveling surface wave. The purpose of this protocol was to demonstrate that an increase in skin tension corresponds to an increase in the wave speed (or equivalently, a decrease in the arrival time) which can be measured in vivo by the Reviscometer ® . Measurements were taken from 0° to 360° in 10° increments giving a total of 36 observations. This method was repeated three times per volunteer and an average was calculated. 1 Using the "natural configuration" data, the direction of the Langer line was identified to the nearest 10° by finding   the direction in which the shortest arrival time was recorded [15]. Then, using surgical tape, the skin was stretched in the direction of the identified Langer line. The "stretched configuration" test was then repeated a further three times and an average was calculated. 1

Anisotropy Measurement
As discussed in the "Introduction", there is a need to quantify skin anisotropy and understand its relationship to skin tension. A number of previous studies have been performed where measures of skin anisotropy are calculated. The most commonly used measure is the ratio of the maximum and the minimum measured value (arrival time or wave speed) [24,3,33,15]. In our study, using the Reviscometer ® , this Anisotropic Ratio (AR) is the ratio of the maximum and minimum RRT values: While this measure can be indicative of the degree of anisotropy, it is also very sensitive to outliers in the data. Furthermore, if measurements are taken from 0° to 360° (as is often the case), this measurement ignores much of the available data and uses only the maximum and minimum values. In this paper we suggest an alternative measure of anisotropy that is less susceptible to individual outliers and considers all measurements from 0° to 360°. We consider all RRT observations and fit an ellipse to them. The eccentricity of this fit ellipse is indicative of the material anisotropy. Assuming an ellipse is an appropriate model for our data, we plot the raw data by allowing the arrival time to be the distance from the origin and the angle to be the angle of inclination from the positive side of the x axis, see Fig. 2a.
Using this representation in Cartesian coordinates, an ellipse can be fit to the raw data using the least squares approach detailed in Ref. [34], which is implemented in the function "EllipseDirectFit" from the R package "conicfit" [35,36], see Fig. 2b. We can then extract the geometric parameters from the ellipse and use them to infer real-world attributes of the skin. All code used can be found in the public GitHub repository accompanying this paper:https:// github.com/matt-nagle/Analysis-of-in-vivo-skin-anisotropyusing-elastic-wave-measurements-and-Bayesian-modelling.
The geometric parameters extracted from each ellipse were: the lengths of the semi-major and semi-minor axes and the angle between the semi-major axis and the positive x axis (tilt angle). The tilt angle of the ellipse provides the direction of the slowest traveling wave and 90° from this is the fastest traveling wave which corresponds to the direction of Langer lines [15].
The lengths of the semi-major and semi-minor axes allow us to calculate both the area, A, and the eccentricity, e, of the fit ellipse using Eqs. 3 and 4 respectively: where a is the length of the semi-major axis and b is the length of the semi-minor axis.
The area relates to an average measure of arrival time in all directions. The smaller the area, the faster all waves are traveling on average. Following Eq. 1, we can relate this wave speed directly to skin stiffness. This parameter is independent of the anisotropic nature of the measurements.
Eccentricity relates to the material anisotropy; an eccentricity of 0 indicates a circle, i.e. the speed of the wave does not vary depending on the angle and the material is perfectly isotropic. As the eccentricity increases from 0, the difference between the slowest wave and the fastest wave also increases, i.e. the skin demonstrates more and more anisotropy. An eccentricity of 1 indicates a straight line which is perfectly transversely isotropic. In practice, eccentricity values in our study mostly fell between 0.5 and 0.9.

Simulation Study
To evaluate the performance of the two different measures of anisotropy (AR vs e) a simulation study was performed. Simulation studies are computer-based experiments that use artificially generated data to examine the performance of different methods. Knowledge of the underlying data generation mechanism enables a thorough evaluation and comparison [37]. In short, simulated data was generated following a known regular shape, random noise was added to the data, then the two measures of anisotropy were compared to the known true values. Ellipses with a fixed value of 160 RRT for the semi-minor axis with eccentricities e = [0.5, 0.7, 0.9] were selected as being representative of our dataset. Noise was added to the ellipses using a random draw from a normal distribution with mean 0 and standard deviation . Four different values of were used, ranging from very low to high amounts of noise, = [1, 10, 20, 30] , see Fig. 3.
For each set of points, both measures of anisotropy (eccentricity of the fit ellipse and the AR) were calculated and stored. This procedure was performed 10,000 times for each value of .

Bayesian Data Analysis
As discussed in the "Introduction", Bayesian methods for data analysis are often appealing as they avoid some of the potential issues with p-values and NHST [26]. In general, in a frequentist approach, a model coefficient is a single deterministic fixed value. On the other hand, in a Bayesian framework, model coefficients are random quantities which are assumed to have probability distributions that convey prior beliefs and the uncertainty around their value. The aim is to perform inference on the "posterior distribution" of the parameters, taking into account both the prior knowledge and the evidence provided by the observed data. Inference in this setting is typically performed using Markov Chain Monte Carlo (MCMC) methods [27], which are employed to derive a large chain of estimates for each model coefficient. Each coefficient estimate in the chain is a draw from the posterior distribution and considering a large number of estimates gives us the shape of this distribution. This allows us to quantify not only the mean value of the coefficient ("posterior mean") but also the uncertainty we have in each coefficient by considering the spread of the distribution.
In contrast to the confidence interval in the frequentist approach, Highest Posterior Density Intervals (HPDI) can be defined which directly relate to a probability, i.e. there is a probability of 0.95 that the true value of the coefficient lies within the 0.95 HPDI. Thus, if we are interested in covariate significance, we can simply examine the posterior distribution of the corresponding coefficient. For example, if 0 is within the 0.95 HPDI for a certain coefficient then we do not have enough evidence to suggest that the coefficient is significantly different from zero and must conclude the covariate does not have a significant effect on our outcome variable. Bayesian approaches are often much more flexible and allow for greater model complexity; this is especially important for the current application as we need to account for correlated target variables, as well as circular target variables which would be very difficult to do with a frequentist approach [27].
To examine the influence that the subject age, subject gender and the applied additional stretch have on the skin properties anisotropy, average stiffness and stiffness in the direction of the Langer line (measured by the eccentricity, area and length of the semi-minor axis of the fit ellipse respectively), a Bayesian multivariate outcome regression model was built. The multivariate approach was selected to account for the fact that e and A are both calculated using the semi-major and semi-minor axes and are therefore correlated. The model is of the form: where is a three-dimensional vector representing the intercept, B is a 3 × 3 matrix of coefficients for age, gender and configuration, and E ∼ N(0, Σ) is the error term. Note that as the eccentricity is a parameter bounded by 0 and 1, a regular linear model would not be suitable as we cannot obtain values of e > 1 or e < 0 regardless of the input values. Therefore, a log of the eccentricity was taken as the outcome variable. For numerical stability the Area was normalised by a factor of 1000.
Inference for a Bayesian model of the form Eq. 5 can be performed via Markov Chain Monte Carlo (MCMC) methods using the function "stan_mvmer" implemented in the R package "rstanarm" [38][39][40]. We use default non-informative Area/1000 Semi-minor Axis priors and obtain draws from the posterior distribution of the regression parameters.

Simulation Study Results
In the section "Anisotropy Measurement" we suggested that the eccentricity of a fit ellipse would be a more robust measure of anisotropy than the ratio of the max and min values. An illustrative example in Fig. 4 shows how outliers could drastically affect the anisotropy ratio, providing misrepresentative results. However, a more systematic evaluation is required to directly compare the two methods, see the section "Simulation Study".
The results of the simulation study can be seen in Fig. 5. Unsurprisingly, when the noise is very low ( = 1 ), both anisotropy measures perform extremely well, they are tightly distributed and centered over the true values. We can see that, as the level of noise increases, the distributions of both measures widen and stray away from the true value. However, we can see that for low to high amounts of noise the eccentricity is far more robust than the AR: the distribution tends to widen but it is still distributed around (or very close to) the true known value. In contrast, as the noise increases, the AR distribution strays further and further away from the true value and the tails of the distribution become more pronounced with extreme outliers appearing. It should also be noted that the AR consistently overestimates the degree Fig. 4 Schematic demonstrating how outliers (in grey) might change the AR dramatically (in this example by increasing the Max RRT and decreasing the Min RRT beyond a range that represents the data) but only slightly elongate the fit ellipse. The red dotted line represents the underlying ellipse with eccentricity 0.7 before noise or outliers were added. The fit ellipse can be seen in blue and has an eccentricity of 0.756. The true AR for the underlying ellipse would be 1.4; however, due to the outliers, the AR in this case would be 2.7.
of anisotropy, introducing a bias into our measurement, whereas the eccentricity provides a more consistent estimate on average (albeit underestimating the actual eccentricity in some scenarios with high noise).
Analysis of particular simulations which resulted in outliers for the eccentricity/AR (where the measure was far from the true value) can be found in Appendix A. From this, and the analysis of the overall distributions, we can conclude that the eccentricity is a more robust and reliable measurement of anisotropy for this type of data.

Influence of Age and Gender on Skin Properties
As discussed in the section "Bayesian Data Analysis", using a Bayesian approach we can examine the significance of the model covariates by examining the posterior distribution for each parameter. A model of the form Eq. 5 was built and the posterior distributions can be seen in Fig. 6.
Let us first consider the influence of age. We can see in Fig. 6a that the 0.95 HPDI (grey area) for the age coefficient does not contain 0. Thus, given the data, there is a 95% probability that age has a significant effect on the degree of anisotropy as measured by the eccentricity. The age coefficient is positive and centred around a posterior mean value of 0.007 (see Table 3). Therefore, there is a positive effect of the age on the log of the eccentricity, i.e. the eccentricity increases logarithmically with age. Furthermore, due to the shape of the log function we can conclude that the eccentricity (degree of anisotropy) increases as age increases with a steep increase from childhood into adulthood. This can also be seen by assuming no correlation between the outcome variables and plotting the age of the subject against the eccentricity (see Fig. 7).
In Fig. 6b we can see that the 0.95 HPDI (grey area) barely includes the value zero. Hence, according to the data, age may not have a significant effect on the average stiffness as measured by the ellipse area. However, it is worth noting that a considerable part of the posterior distribution is above the zero threshold, which is in line with the many reports in the literature indicating that skin stiffness increases with age. We must keep in mind that the area here refers to the "average stiffness" across all directions.
In Fig. 6c, we see the posterior distribution of the effect of age on stiffness in the direction of the Langer line (max stiffness) measured by the length of the semi-minor axis. We can see the age of the subject is significant according to the data. The age coefficient is negative and centered around a posterior mean value of − 0.391. Therefore, the data suggests there is a negative effect of age on the length of the semi-minor axis, i.e. as age increases, the length of the semiminor axis decreases (arrival time of the wave decreases in the direction of the Langer line indicating increased stiffness of the skin in that direction), as expected. Now, let us consider the influence of gender. We can see in Fig. 6a that the 0.95 HPDI (grey area) for the male coefficient contains 0. Thus, it is unlikely that the gender coefficient has a significant effect on the eccentricity (degree of anisotropy), i.e. according to the data, there is no significant Note that in order to capture the extreme outliers in the AR measure plot (right), the y axis was log transformed. It is clear that for low to high amounts of noise the eccentricity distributions are centered around (or are close to) the true value, whereas the AR distributions stray further and further away from the true value with increasing noise. Simulation Results for AR difference between the degree of anisotropy in males versus females. In Fig. 6b we can see that the data shows evidence of a difference between the average stiffness of male subjects versus female subjects. We can also see that, using the category female as the baseline, the coefficient for the shift in area due to the male category is negative and centred around a posterior mean value of − 16.725. Therefore, the area of the fit ellipse is on average smaller (average skin stiffness is higher) for males than females. Finally, in Fig. 6c we can see that the data provides evidence for a significant difference between the length of the semi-minor axis (arrival time in the direction of the Langer line) for male and female subjects. The male coefficient is negative and centered around a posterior mean value of − 15.174. Therefore, the skin stiffness in the direction of the Langer line is on average higher for males than females.  6 Posterior Distributions for the age, gender and configuration regression coefficients from Eq. 5. Note that the shaded region between the two vertical black lines represents the 0.95 HPDI and the vertical dashed red line denotes the location of 0. If 0 is within the 0.95 HPDI, there is not enough evidence to say the covariate has a significant effect on the outcome variable. The posterior mean and HPDI intervals for all distributions can be seen in Table 3.

Influence of Skin Tension
To explore the effect of skin tension on the speed of travelling surface waves and determine the potential efficacy of the Reviscometer in evaluating skin tension, an additional stretch (in the direction of the fastest wave) was manually applied to the skin and the same measurement procedure was carried out. Following this procedure, in theory, we would expect that: 1. The tilt angle should be conserved, provided the stretch is applied in the direction of the fastest traveling wave (along the Langer line).
2. The length of the semi-minor axis should decrease due to the increase in wave speed along the direction of the applied stretch. 3. The eccentricity should increase, as we would expect that the stretched data would be more anisotropic due to the additional applied stretch and decrease in semiminor axis length.
For an example of this behaviour see Fig. 8. Following the same analysis of age and gender in the section "Influence of Age and Gender on Skin Properties" we can also examine the effect that the additional applied stretch had on the length of the semi-major axis and the eccentricity by examining the posterior distribution of each parameter.
Note that as our covariate is a categorical variable with two outcomes (either natural or stretched), the "stretched coefficient" is the shift in outcome variable from the baseline natural configuration to the stretched configuration. From  Fig. 6, we can see that the data provides evidence of a significant difference between the natural and stretched configurations in the log eccentricity and the length of the semi-minor axis. The 0.95 HPDI (grey area) does not contain 0 and the posterior distributions are concentrated far from this value. Furthermore, we can see that the shift for the stretched configuration is positive for the log of the eccentricity outcome (posterior mean of 0.137) and negative for the length of the semi-minor axis (posterior mean of − 26.808). This demonstrates that the eccentricity of the fit ellipse increases and the length of the semi-minor axis decreases on average from the natural to the stretched configuration, as expected.
Finally, to examine the effect that the additional applied stretch had on the angle of the fit ellipse, which indicates the direction of the Langer line, a model was built of the form: where is the intercept, is the coefficient for the configuration and is the error. Note that as our covariate is a categorical variable with binary outcomes (natural or stretched), the intercept represents the baseline category (natural) and is the shift in angle for the alternative category (stretched).
Note that here we have a circular response variable and cannot build a regular linear model. Following the discussion in Ref. [41], we fit a circular regression model using the projected normal distribution within a Bayesian framework, using the MCMC methodology implemented in the function "bpnr" from the package "bpnreg" [42]. The posterior distributions for the natural and stretched configurations can be seen in Fig. 9. We see that the 0.95 HPDI overlap for the natural and stretched configurations. Therefore, in the data, there is no indication of a significant difference in the posterior circular means of the two configurations. This indicates that the configuration does not affect the angle of the fit ellipse, hence the direction of the Langer line, as expected. However, there is a degree of uncertainty as the posteriors do not overlap completely, so for some subjects there may be variation in angle due to the additional stretch, but the difference in configuration does not appear to be strong in the data.

Discussion
As discussed in the section "Anisotropy Measurement" we have demonstrated that the eccentricity of a fit ellipse is a better measure of the anisotropy, as it is more robust to outliers and noise, than the commonly used anisotropic ratio [3,21,24]. Robustness to noise and outliers is an important attribute for biological measurements and particularly in vivo skin measurements as we expect experimental error, patient variability and subject movement to impact data collection. Using our new measure of anisotropy, the relationship between the skin anisotropy and age was explored, see the section "Influence of Age and Gender on Skin Properties". We found that as age increases, the degree of anisotropy also increases, with a steep increase occurring from childhood to adulthood. Many early authors previously reported only a weak dependence of elastic wave speed with age [23][24][25], but as noted by Ruvolo et al., these studies employed poor resolution angular data (measurements taken every 45° only). More recent studies have noted an increase in anisotropy with age [43,44], but none of these studies included infants and therefore the steep increase in anisotropy which occurs from childhood into adulthood would not have been captured. Ruvolo et al. did include infants within their study and report an exponential increase in the anisotropic ratio with age, while in the current study, we have found a logarithmic increase in anisotropy with age. The reason for the different observations may be, in part, due to the use of the anisotropic ratio in Ruvolo et al. and also due to the fact that participants were divided into 5 age ranges for the purpose of the statistical analysis. In the current study, age was considered a continuous variable and hence provides a richer insight into the true relationship between age and anisotropy. This detailed information can provide evidence as to the validity of applying universal cosmetic surgery practises to very young or elderly patients, where we expect the level of anisotropy to vary significantly.
It has variously been reported in the literature that skin stiffness increases [45,46], decreases [25,3,21,43], or is not affected at all [24] as the age of the subject increases. In the section "Influence of Age and Gender on Skin Properties", we concluded that age does not appear to have a significant effect on the area of the fit ellipse (our measure of overall stiffness, independent of direction). However, by virtue of considering the area of the ellipse independent of its shape we are "averaging" over the known fact that the skin of older subjects is more anisotropic. And while the area of the ellipses does not seem to be affected by age, the length of the semi-minor axis (i.e. max stiffness) does appear to be affected (see Fig. 6c). This indicates that as the age of the subject increases, the stiffness of the skin increases in the direction of the Langer line only. In this context, our results agree with [44,46] who found that there is increased stiffness along the direction of Langer Lines with age. The current results are, however, in direct disagreement with Ruvolo et al. and Hermanns-Lê et al. who found that age influenced the maximum RRT (equivalent to the length of the semi major axis in the current study) but not the minimum RRT (equivalent to the length of the semi minor axis in the current study). It is possible that these differences may have arisen due to the method of data analysis e.g. in contrast to the length of the semi-major and semi-minor ellipses fit to the circular data, the use of maximum and minimum RRT values does not account for outliers and does not offer a robust measurement. Contrasting results with other studies may be as a result of reporting "average" stiffness which does not consider the significant anisotropy of skin (in particular elderly skin). We should acknowledge also a degree of uncertainty with the conclusion that the "average stiffness" (area of ellipse) is not affected by age, as 0 is just inside the 0.95 HPDI and a large portion of the posterior distribution is concentrated above this value (see Fig. 6b). Therefore, perhaps it is unsurprising that there is disagreement in the literature on this point.
Similarly, we acknowledge a degree of uncertainty with the conclusion that males have stiffer skin than females on average since for both the Area (average stiffness) and the length of the semi-minor axis (arrival time in the direction of the Langer line), 0 is just outside the 0.95 HPDI (See Fig. 6). This result is consistent with [46] who state that males have stiffer skin (albeit only within certain age ranges, 21-40 years) but contrasts with [47] who state that women have stiffer skin. These contrasting results found in the literature may be as a result of the complex interactions between the influence of age, gender, and body location on the skin stiffness.
On average, the angle of the ellipse is conserved from the natural to the stretched configuration; however, there is quite a wide range of individual behaviours (see Figs. 9 and 11). While some extreme behaviours like subjects going from close to 0° to close to 180° or vice versa can be explained due to the circular nature of the variable, other variations are likely due to experimental error. As discussed in the section "Data Collection", a fixture was used to take measurements every 10°. Therefore, the identified direction of fastest wave speed (i.e. Langer line) was to the nearest 10°, i.e. an accuracy of ±5 • . Once the direction of the Langer line was identified, tape was applied manually to stretch the skin in that direction, which also introduces a source of error. This in turn could affect the fundamental mechanics of the probe, introducing further uncertainty, and explaining some of the uncertainty in this result.
During the calibration of the Reviscometer ® the material being tested was assumed to be linearly elastic, isotropic and unstressed. While these assumptions are true for the reference materials used (see Table 1), the application of Eq. 1 to in vivo skin is not strictly valid, because of the pretension and the oriented collagen fibers. Future mathematical relations relating elastic wave speed to stiffness should seek to take these effects into account. Such a relationship may provide a means to explicitly determine both in vivo skin tension and skin stiffness using elastic waves.
A core assumption of this study was that the ellipse is a good fit for our data. We believe this is a reasonable assumption due to the flexibility of an ellipse fit and the data collection procedure which should result in axisymmetric measurements (due to the repeated measurements from 180° to 360°) about and perpendicular to the Langer lines. It should be noted, however, that the raw measurements from some subjects did not exhibit this symmetry and thus the ellipse fit was poor in those cases.
In this paper, an in vivo elastic wave technique was employed to investigate the role of age, gender and skin tension on both skin anisotropy and skin stiffness measurements on a sizeable population. By fitting an ellipse to angular data and reporting its eccentricity, we have proposed a more reliable, robust and informative metric of anisotropy than the classic "anisotropic ratio" favoured so far in the literature. Using a Bayesian approach, we have concluded that skin anisotropy increases logarithmically with age, with a steep increase occurring from childhood into adulthood. Furthermore, the maximum stiffness of skin increases linearly with age, but this increase is only seen along the direction of Langer Lines. We have also concluded that gender does not significantly influence the degree of skin anisotropy, but that both the average skin stiffness, and skin stiffness along the direction of Langer Lines is higher for males than females. Finally, we have also concluded that both the skin anisotropy and skin stiffness measurements are significantly influenced by the level of skin tension present. This suggests that in vivo elastic wave measurements may be a suitable method for inferring in vivo skin tension.
To the best of our knowledge, this is the first study which uses a sizeable sample of in vivo subjects and modern Bayesian statistical analysis to evaluate the effect of age and gender on in vivo skin anisotropy. This dataset will provide a useful reference to those wishing to evaluate the effect of subject specific parameters such as age and gender on the anisotropic response of skin.

Appendix A: Simulation Study-Examination of Poor Performance
In addition to examining Fig. 5, we can go one step further by analysing examples of poor performance for each measure. For example, lets choose the medium ellipse with medium noise scenario ( e = 0.7 and = 20 ) and examine the performance of both the eccentricity and the AR for two different simulations where: 1. Eccentricity is furthest from the true value (i.e. the simulation that produced the largest outlier in the eccentricity distribution). 2. AR is furthest from the true value (i.e. the simulation that produced the largest outlier in the AR distribution).
This allows us to explore the conditions in which each measure performs poorly and how the other measure performs in those conditions, see Fig. 10. As we have chosen the middle ellipse, the true eccentricity and AR values are 0.7 and 1.4003 respectively. In Simulation 668 (Fig. 10b), the eccentricity of the fit ellipse was 0.7196 and the AR was 3.6236. For this dataset the ellipse performs very well and the poor performance in the AR is driven solely by the two extreme points which are used to calculate the AR but are not representative of the data as a whole.
In Simulation 583 (Fig. 10a), the eccentricity of the fit ellipse was 0.4672 and the AR was 2.0644. For this dataset neither the AR or the eccentricity are close to their true values, the AR significantly overestimates the degree of anisotropy using two extreme values that are only 30 degrees apart (rather than 90° which would make physically intuitive sense). The poor performance in the eccentricity is because, by chance, the random noise pulled the points along the semi-major axis closer to the origin and pushed the points along the semi-minor axis away from the origin. Thus, even though we know the "true" eccentricity before noise to be 0.7, it could be argued that the calculated eccentricity of 0.4672 is a better measure of the anisotropy as it is representative of the data. This, along with the analysis in the section "Simulation Study Results" allows us to conclude that eccentricity is a more robust measurement of anisotropy for this type of data.

Appendix B: Individual Subject Behaviours
As well as looking at the average behaviours in Fig. 6 and Fig. 9 we can examine the specific behaviour of individual subjects using spaghetti plots (see Fig. 11). Most subjects exhibit the expected behaviour for the length of the semiminor axis and the eccentricity with only a small minority of subjects that have alternative behaviours (see Figs. 11b and 11c). For the angle, we can see quite a range of different behaviours (see Fig. 11a) but on average, most subject angles remain more or less the same. Note that some subjects go from close to 0 • to close to 180 • or vice versa, this behaviour is equivalent to a small increase/decrease in angle due to the circular nature of the measurement.  (c) Fig. 11 Spaghetti plots where the red circles denote the mean and the light red region is the 95% confidence interval of the mean. a There is a wide variety of behaviours for the angle but on average the angle remains constant. Note that some subjects go from close to 0 to close to 180 or vice versa, this behaviour is equivalent to a small increase/decrease in angle. b There is some variance in behaviour but the semi-minor axis decreases for the majority of subjects. c There is some variance in behaviour but the eccentricity increases for the majority of subjects from the natural to the stretched configuration.