Compressive Viscoelasticity of Freshly Excised Mouse Skin Is Dependent on Specimen Thickness, Strain Level and Rate

Although the skin’s mechanical properties are well characterized in tension, little work has been done in compression. Here, the viscoelastic properties of a population of mouse skin specimens (139 samples from 36 mice, aged 5 to 34 weeks) were characterized upon varying specimen thickness, as well as strain level and rate. Over the population, we observed the skin’s viscoelasticity to be quite variable, yet found systematic correlation of residual stress ratio with skin thickness and strain, and of relaxation time constants with strain rates. In particular, as specimen thickness ranged from 211 to 671 μm, we observed significant variation in both quasi-linear viscoelasticity (QLV) parameters, the relaxation time constant (τ1 = 0.19 ± 0.10 s) and steady-state residual stress ratio (G∞ = 0.28 ± 0.13). Moreover, when τ1 was decoupled and fixed, we observed that G∞ positively correlated with skin thickness. Second, as steady-state stretch was increased (λ∞ from 0.22 to 0.81), we observed significant variation in both QLV parameters (τ1 = 0.26 ± 0.14 s, G∞ = 0.47 ± 0.17), and when τ1 was fixed, G∞ positively correlated with stretch level. Third, as strain rate was increased from 0.06 to 22.88 s−1, the median time constant τ1 varied from 1.90 to 0.31 s, and thereby negatively correlated with strain rate. These findings indicate that the natural range of specimen thickness, as well as experimental controls of compression level and rate, significantly influence measurements of skin viscoelasticity.


Introduction
The skin plays a critical role in protecting the musculoskeletal system and internal organs and serves to detect external stimuli. The skin's mechanical properties greatly impact how these

Materials and Methods Overall
Uniaxial compression experiments on flat, cylindrically cut skin samples utilized controlled displacement ramped into the skin surface to collect time-force displacement data. For the purposes of analyzing the data, we generated material parameters of the quasi-linear viscoelasticity (QLV) model [24]. To decouple viscoelasticity from other factors such as material hyperelasticity and stimuli conditions, we obtain QLV parameters from a hyper-viscoelastic constitutive model and only examined its viscoelastic parameters, rather than comparing the force traces alone. The measured specimens were from 139 skin samples freshly excised from 36 mice, ranging 5.7-34.3 weeks in age, and from three skin sites on the hindlimb: distal (Distal), proximal on nerve trunk (NT), and off nerve trunk (OffNT). The three skin sites were selected due to differences in their thickness and underlying fascial structures [12].
Three independent variables were examined for their correlation with skin viscoelasticity: skin thickness (range from 211 to 671 μm, natural variation due to hair cycles over selected age), strain level (steady-state stretch λ 1 from about 0.2 to 0.8), and strain rate (median values of 0.06, 3.54 and 22.88 s −1 ). Strain level is defined as = |ln (λ)| = − ln (λ) in uniaxial compression, where denotes strain, λ denotes the stretch of material calculated from deformed thickness l divided by original thickness l 0 , l ¼ l l 0 . The strain level dependency was analyzed using stretch, which aligns with finite deformation theory [25] and negatively correlates with strain level in the case of compression. The rate of strain was defined as _ " ¼ d" dt . Finally, to validate that the viscoelastic parameters obtained in the skin compression experiments could be used to predict the behavior of the skin in a different context, we performed a secondary experiment with fresh mouse skin where we changed the stimulus, specimen size and different boundary condition. Finite element analysis was used to predict the results of this experiment.

Ethics Statement
All animal use was conducted according to the National Institutes of Health Guide for the Care and Use of Laboratory Animals and was approved by the Institutional Animal Care and Use Committee of Columbia University (protocol AC-AAAC1561).

Apparatus
Compression tests were conducted on a custom-built test machine (Fig. 1A), described in depth elsewhere [12]. Briefly, the test machine's components include a platen of aluminum (3 mm thick and 2.54 cm dia.) attached to a vertical load sled, which was driven by a motion controller (motion controller: Newport, Model ESP300, Mountain View, CA; linear stage: Newport, Model ILS100. Reaction force at the platen was measured by a loadcell (Honeywell, Miniature Model 31, Columbus, OH) with full capacity of 2.45 N mounted between the platen and vertical load sled, and its position was tracked by a laser displacement sensor (optoNCDT Model ILD 1402, Micro-Epsilon, Raleigh, NC) with a resolution of 1 μm. Both force and displacement were sampled at 1 kHz. The platen compressed the skin specimens against a rigid platform parallel to its surface, with sides of the cylindrical skin unconfined. The apparatus was equipped with a closed-loop control system integrated to maintain temperature of 32 degrees Celsius, consistent with prior works [26], using a BASIC Stamp microcontroller (Parallax Inc., Rocklin, CA) and associated electronic transistors and heating elements.

Animals and dissection
The animal preparation and dissection protocol has been described previously [12]. Skin samples were obtained using a 6-mm diameter punch (Acuderm Inc., Ft. Lauderdale, FL) after skin specimens were dissected from the mouse hindlimb. Three sampling sites (Fig. 1B) were chosen at distal end of the hindlimb (Distal), and the proximal end of the hindlimb both off (OffNT) and on (NT) the saphenous nerve trunk because these sites appear to be categorically differentiable in terms of thickness. Specimens were harvested from a total of 36 mice, at ages ranging from 5.7 weeks to 34.3 weeks and body weights ranging from 15.9 grams to 61.4 grams. A total of 139 skin samples including 46 from Distal, 46 from OffNT and 47 from NT areas were harvested.

Skin test procedure
For all specimens, we first set the starting position of the platen to ensure it was positioned above the flat skin surface by observing the reaction force. Next, displacement-controlled compression was applied with a ramp-up phase at a velocity to achieve the desired strain rate, a hold phase at the maximum load position for 6 seconds (note that only the first 5 seconds were used in analysis to avoid analyzing the ramp-off response), and an unloading phase of the same rate as the ramp-up. Multiple repetitions of same loading protocol were applied to the specimen, where the 6th run for each skin specimen was analyzed and the first 5 repetitions were used as pre-conditioning to minimize the variance due to stress history [27]. Fig. 2 demonstrates a typical experimental procedure, where strain rates are varied. Synthetic interstitial fluid (SIF) [28] was added via eye dropper to prevent drying of the skin. The reaction force at the platen was measured by a loadcell and platen position measured by a laser displacement sensor. The recorded force trace was then used to determine the point of contact (Fig. 2C). A "light-contact point" of the platen to the specimen surface was determined at the moment when reaction force on the platen exceeded 0.1 N. After that, the "contact point" was defined as distance from the platen to the rigid table at the "light-contact point", timed by a correction coefficient of 1.3. Then, the specimen thickness (l 0 ) and deformed thickness (l) was each defined as the distance from the platen to the rigid table at and after the "contact point".
Three different experimental paradigms were used to elucidate viscoelastic properties of multiple skin samples under compression. First, we measured skin under the same level of strain from 44 skin samples of naturally varying thickness. Maximum indentation depths were determined by manually searching for an instantaneous reaction force around 2 N, which is  Example run of the compressive test procedure for one skin specimen when varying strain rates. A: Position of the compression platen over time, as measured by its distance from the fixed platform. B: Reaction force at the compression platen. C: Magnified view of reaction force and platen position for Compression 6, demonstrating that "skin thickness" was defined by "contact point" as determined from the force trace. The platen was moved into the skin with an acceleration of 0.06 s −1 for each of the first 10 repetitions. Then, 10 additional compressions were performed at 22.88 s −1 . The 6th compression was analyzed in each sequence of 10 compressions. the approximate magnitude to generate a level of stretch of 0.6, similar to indentation in neurophysiological studies [26]. The velocity of the compression platen was 1 mm/s to achieve a strain rate of about 3.54 s −1 . Second, we applied similar procedures to 41 skin samples, but linearly varied steady-state stretch for each specimen. The minimum stretch level was set when the responsive force recorded at the loadcell was above zero, determined by its magnitude being one standard deviation above background noise, and the maximum stretch level was set when the maximum responsive force reached 2.45 N. Any force above this level was avoided to prevent damage to the skin or instrumentation. Fig. 3 demonstrates the difference in stretch level between the first and second experimental paradigms. In the third experiment, strain rate was varied. Two rates (medians are 0.06 and 22.88 s −1 ) were applied to 54 specimens, and the data were analyzed together with the first experiment (3.54 s −1 ) to constitute three strain rates at different orders of magnitude. The low strain rate ( _ " ¼ 0:06 s À1 ) is comparable to activities with long relaxation times, such as standing and lying in bed, where one can still perceive the mattress even after several minutes. The medium strain rate ( _ " ¼ 3:54 s À1 ) corresponds to typical light-touch activities, such as typing on a keyboard. The high strain rate ( _ " ¼ 22:88 s À1 ) corresponds to impact loading, which one perceives to avoid imminent danger. The three strain rates correspond to 0.01, 1 mm/s and the fastest moving velocity of our test machine.

Constitutive laws
The QLV model [24] was used here to fit the data to be analyzed. Given that our test was uniaxial, we only considered the one-dimensional situation. The QLV model utilizes the stretch-time curve to calculate the response stress, shown as ! @s e ðlÞ @l @lðt 0 Þ @t 0 dt 0 : ð1Þ For more details of the constitutive laws please refer to Appendix I.

Fitting experimental data to constitutive model
To attain the parameters of the constitutive model, we fit the model to the stress-stretch-time curves calculated from the experimental data. The stretch value was calculated by dividing deformed thickness l over original thickness l 0 , i.e. l ¼ l l 0 (l, l 0 defined in Section 2.4). Recorded experimental force data were converted to stress values by dividing force over area, calculated from Equation (7), with the assumption that the volume of the specimen remained constant because it is nearly incompressible [24]: where A 0 is the area of the punch, a 6-mm diameter circle. For the detailed numerical algorithm used for fitting, please refer to the Appendix I. QLV model parameters were then adjusted to fit to the stress-time and stretch-time measurements (Fig. 4) taking the number of terms n = 1 and n = 2.
The stress and stretch history over the whole time window, including phases of both the ramp (from contact to peak stress) and hold (from peak stress to 5 seconds after), were used in the fitting to account for relaxation during loading, similar to Laksari et al. [29]. Because the number of data points in the hold phase were much greater than that in ramp phase, we used an R 2 value as the equally-weighted sum between the R 2 from fit of the ramp phase and hold phase. For each specimen, the weighted R 2 was maximized through a constrained nonlinear optimization (fmincon, in the MATLAB Optimization Toolbox) using the SQP algorithm. The reduced chi-square value (w 2 red ) and the residual standard deviations (σ res ) were also checked to assure the quality of fit.
Then, fitting of both one-and two-term models was performed in two steps: 1. We fit the experimental data with all free viscoelastic parameters (2 and 4 for one-and twoterm models, respectively) and all free hyperelastic parameters (2 for each model), with manually chosen initial values for the optimization algorithm.

We fixed all time constant parameters (τ) and initial shear modulus (μ) to median values found in
Step 1 between specimens tested under same strain rate. For the initial values for free parameters, corresponding median values found in Step 1 were used.

By fixing certain parameters in
Step (2), the total number of free viscoelastic parameters were reduced to 1 for one-term (G 1 ) and 2 for two-term models (G 1 , G 1 ) respectively, and also included only 1 hyperelastic parameter (α) in each. For analysis of the distribution of time constants (τ), fitting results from Step 1 were used, and for that of stress ratios (G) results from Step 2 were used.
The fitting of QLV model to experiment data for each skin specimen was performed and results were listed in Table 1, which shows high R 2 values, w 2 res close to 1 and low σ res (< 1 kPa, compared to peak stress of about 50 kPa in the 2 nd experiment), indicating a good fit. Data in Table 1 and Fig. 4 reveal the trade-off for increasing the number of free parameters from 1 to 2 was attaining only a small improvement in fit. Thus, we decided to use the one-term model so that comparisons between specimens were easier with only a single free parameter. More importantly, by strictly controlling the number of free parameters, we minimized the non-unicity of the fitting.    . Black line shows the modeled prediction, and gray data points show the experimental data. The average weighted R 2 value for the one-term case for the three strain-rates is 0.86 while the R 2 value using the two-term case is 0.93. Therefore, the tradeoff is that the number of free parameters increased from 1 to 2, versus attaining a slight improvement in the fitting, and for this reason we chose the one-term case. doi:10.1371/journal.pone.0120897.g004

Results
The parameters returned by fitting the one-term QLV model revealed that skin viscoelasticity is highly variable between specimens, yet correlates with the three independent variables. Specifically, the residual stress ratio G 1 positively correlates with skin thickness and stretch level, and the time constant τ 1 negatively correlates with strain rate.

Positive correlation between thickness and residual stress ratio
In the first experiment where the skin thickness naturally varied, residual stress ratio G 1 positively correlated with skin thickness, with Pearson correlation coefficient of 0.883 (Fig. 6A). Linear regression with residual stress ratio G 1 as a dependent variable was performed, which returned p < 0.001 for independent variable thickness l 0 , and G 1 = 9.997 × 10 −4 μm −1 Á l 0 + 0.077.

Strain-level dependency
In the second experiment where the change in stretch delivered accompanied skin thickness variation, the residual stress ratio G 1 was found to positively correlate with skin thickness (Fig. 6B), and moreover, also positively correlated with stretch level (Fig. 6C). Multilinear regression with residual stress ratio G 1 was also performed, which returns p < 0.001 for independent variable stretch λ 1 , p < 0.001 for independent variable thickness l 0 , and G 1 = 0.810 Á λ 1 + 4.25 × 10 −4 μm −1 Á l 0 − 0.074, indicating that both skin thickness and stretch level are positively correlated with residual stress ratio and thus contribute to variability. Correlations between skin thickness/stretch level and residual stress ratio (G 1 ). A: In the first experiment where only thickness varied, the steady-state residual stress ratio (G 1 ) correlates with increasing skin thickness, n = 44. Linear regression (solid line) with residual stress ratio G 1 as the dependent variable was performed, which returns p < 0.001 for independent variable thickness l 0 , and G 1 = 9.997 × 10 −4 μm −1 Á l 0 + 0.077. In the second experiment where both thickness and strain level varied, the residual stress ratio (G 1 ) correlates with both B: stretch and C: skin thickness. Note that the two correlations are independent from each other because there is no correlation between stretch and skin thickness. Multilinear regression with residual stress ratio G 1 was also performed, which returns p < 0.001 for independent variable stretch λ 1 , p < 0.001 for independent variable thickness l 0 , and G 1 = 0.810 Á λ 1 + 4.25 × 10 −4 μm −1 Á l 0 -0.074. Note that in B and C, solid lines are single-linear regressions for residual stress ratio with respect to stretch and thickness respectively.

Strain-rate dependency
In the third experiment, the strain rate largely varied (median _ " ¼ 0:06, 22.88) and combined with data from first experiment (median _ " ¼ 3:54), we found that the same magnitude of relaxation takes place at significantly shorter time constants at higher strain rates (Fig. 7A). Linear regression was performed with time constant τ 1 as the dependent variable and strain rate _ " as the independent variable, which yielded a significantly negative correlation (p < 0.001) between the strain rate and time constant. Using the same regression but replacing the dependent variable from time constant τ 1 with residual stress ratio G 1 , we found that the strain rate did not significantly affect residual stress ratio (p = 0.988 > 0.05). With a closer examination of the distribution of time constants and residual stress ratios (Fig. 7B-D), we noticed that the distributions of time constants notably skewed to the left as strain rate increased, whereas the distribution of residual stress ratios did not show systematic changes.

Discussion
This work shows, for the first time with mouse skin under compression, that skin's viscoelasticity is highly variable (relaxation time constant τ 1 = 0.19 ± 0.10 s and steady-state residual stress ratio G 1 = 0.28 ± 0.13) among a population of skin specimens (n = 139). However, we found systematic correlation in three cases: 1) the residual stress ratio G 1 positively correlates with skin thickness (p < 0.001), 2) residual stress ratio positively correlates with stretch level (p < 0.001), in other words, negatively correlates with strain level and 3) the time constant τ 1 negatively correlates strain rate (p < 0.001). Overall, these findings shed light on the natural range of between-specimen variance under compression, and reveal how experimental controls of strain level and rate can influence measurement of the same specimen.
A small, secondary experiment with fresh mouse skin was performed to validate that the viscoelastic parameters obtained in the skin compression experiments could be used to predict the behavior of the skin in a different context. In particular, using the viscoelastic parameters obtained with the flat plate, we sought to predict the force relaxation of a 1.5 mm probe indented into a skin specimen of different cut-out size (8 mm as opposed to 6 mm), for two indentation depths. This required a compression experiment with mouse skin, as well the use of a finite element model. As denoted in Appendix IV, the force relaxation predicted by the FE model well agrees with experimental data, with an average R 2 = 0.932.
We found that as thickness decreases, residual stress ratio decreases, which means the skin relaxes to a greater extent. This finding agrees with a study by Escoffier et. al [16], who reported that relaxation time decreases as people age, and we know that skin thickness decreases with aging [15]. Also, we identified that the residual stress ratio decreases with lower levels of stretch, i.e., higher strain levels, which echoes Funk et. al [21] who reported the same effect in ankle ligaments. The work herein is the first to report a decrease in time constant under a faster strain rate from biological measurements.
Although the dependency of the skin's mechanical properties on strain and strain-rate is constitutively defined as material non-linearity, the dependency on skin thickness indicates that skin specimens of varying thickness are essentially different materials. Additional analyses indicate that the dependency of skin viscoelasticity on thickness and strain level are neither from frictional edge effects (computational finite element analysis, Appendix II) nor from different dermis/epidermis thickness ratios (statistical regression, Appendix III).
Our results from mouse hindlimb skin are comparable to prior tests of compression with pig dorsal skin [11], exhibiting similar time constants within a 5-second time-scale (median τ 1 = 0.18 s from our one-term model fit compared to τ 1 = 0.57 s on pig skin) and residual stress ratios (median G 1 = 0.284 herein, compared to G 1 = 0.234 on pig skin). However, if we compare the reduced relaxation functions of skin under compression with those of rat skin under tension [4], the compression curves are clearly distinguishable by their significantly smaller residual stress ratio G 1 (Fig. 8). Another key difference compared with that prior work is our use of skin from the hindimb, instead of dorsal skin, which is more commonly measured. The measurement of hindlimb skin is vital for studies of the sense of touch [26], known to be dependent on skin mechanical properties [18]. In particular, slowly adaptive type I (SAI) mechanosensitive afferents, essential for our ability to discriminate edges and curvature [30], display firing rate decay under constant displacement stimuli. This phenomena is known as adaptation and is dependent, in part, on the skin's viscoelastic relaxation [31]. We chose a hold phase at the maximum load position of 5 seconds to align with such adaptation and the typical length of neurophysiological recordings from SAI afferents [26]. Therefore, one would need to be careful in extrapolating the conclusions of this work outside of the chosen time window.
Furthermore, the results of this work give important insights into issues currently being examined in the field of tactile mechanotransduction. SAI adaptation may carry information about a mechanical stimulus, for example, an object's compliance. Since thinner skin relaxes more than thicker skin, these data predict that the neural response from a population of SAI afferents in thin skin might adapt their firing rates to a greater extent than a similar population in thicker skin. This could negatively affect the ability of those with thin skin (e.g., the elderly population) to accurate assess tactile stimuli. In concordance with this, it is known that tactile acuity decreases with age [32]. Studies investigating changes in tactile sensation with aging or after injury usually focus on neuronal causes, but our results suggest skin mechanics might also contribute to changes in tactile sensation. Our understanding of such mechanical propertiesat the level of macro-scale compression-is important to develop realistic models of touch stimuli for haptic technology [18]. The solid line shows median data from work presented here on mouse hindlimb skin fitted to the one-term model, in the first experiment with median strain rate " _ ¼ 3:54 s À1 . The dotted-dash line gives a measurement from pig dorsal skin [11]. The dashed and dotted lines are both from rat skin, but the dashed line function is attributable only to collagen elements in the skin while the dotted function is only elastin elements [4]. Note that the skin in compression relaxes more than skin in tension. The results presented herein are based the assumption of a spatially homogeneous constitutive model; however, the skin is a heterogeneous and anisotropic material, and it is yet unclear what microscopic mechanisms underlie the nonlinear viscoelasticity we observe at the bulk level. Sub-micron studies have begun to suggest that individual skin layers indeed exhibit different degrees of viscoelasticity [9]. This may indicate that viscoelastic nonlinearity at the bulk level are dominated by one or more specific layers, such as the dermis, or a specific constituent, such as the interstitial fluid.
It is worth noting that there are some anatomical differences among various types of skin. The structure of skin differs between mouse hairy skin, our testing site, and glabrous skin. Hairy skin is composed of a thin epidermis that involutes deep into the dermis to form hair follicles. By contrast, glabrous skin, which lacks hair follicles, has a thick epidermis with undulating ridges at the dermal-epidermal junction. Human skin comprises the same fundamental layers as mouse skin with different thickness for each layer, with the exception that the muscular layer of panniculus carnosus in mouse skin does not exist in most areas of human skin [33]. In both species and both types of skin, the density and structure of the layers changes over the course of an animal's life, as the dermal papillary ridges flatten with age [34] and hair follicles undergo growth cycles [35]. While we have standard testing data for murine skin [4,12], the existing literature on human specimens covers only in-vivo viscoelastic measurement with complex stress fields, for example, Krueger et al. [13] investigated how viscoelasticity changes with aging using a Cutometer. Future work on human skin specimens are needed to provide hyperviscoelastic constitutive parameters, and the contribution of each layer to the skin's viscoelastic nonlinearities and the changes in these properties with age is yet to discover for both species, in order to be used for numerical simulations to better aid clinical practice.
Our work suggests that normal features of the neuronal response could be mediated by skin mechanics. In particular, we hypothesize that SAI afferents may adapt their firing rates more quickly to strong stimuli than to weak stimuli, since the skin relaxes more under high-strain conditions. Such changes in neuronal firing could be one mechanism by which the nervous system gains information about stimulus properties. Furthermore, SAI afferents may adapt their firing rates more quickly to faster stimuli than to slower stimuli, since the skin relaxes more quickly under higher strain rates [31,36]. This said, one must also note that intrinsic neuronal properties play a role in the overall adaptation of the mechanosensitive response independent of the skin's response. These results suggest a need to carefully control stimulus magnitude and velocity in performing electrophysiology experiments with tactile stimuli [18].

Appendix I: Details of constitutive model selection and numerical implementation
Hyper-and visco-elastic models have been adopted to fit the behavior of the experimental data, as previous efforts [11] have shown that skin under compression is hyper-viscoelastic. On a macro-scale, most biological tissues are viscoelastic [24] and have well-developed material models depending on the deformation level. Under small deformation, various spring-dashpot models have been used, including the most commonly used Kelvin-Voigt model, a standard linear solid model and generalized Maxwell model (i.e., Maxwell-Wiechert model) [37]. As biological tissue often undergoes finite deformation, these linear models must be modified to incorporate hyperelastic components. Two of the most popular models are the quasi-linear viscoelastic (QLV) model [24] and parallel-network viscoelastic (PNV) model [38]. Although the PNV model yields accurate and stable strain-energy outputs, the QLV model is more popular because the parameters are typically easier to interpret and it has a longer history (Fig. 9). For the QLV model, a convolution integral was used to calculate stress from strain data, Gðt À t 0 Þ @s e ðlÞ @l @lðt 0 Þ @t 0 dt 0 ; ð3Þ where t and λ denote time and stretch. The instantaneous elastic function of material is σ e (λ), where herein we utilized a 1 st -order Ogden form of the hyperelastic strain energy function [39], with μ and α being the material constants, μ also known as instantaneous elastic modulus. G(t) is defined as the reduced relaxation function, and is in the form of a Prony series, where τ i were the time constants associated with stress relaxation ratio G i , and G 1 was the residual stress ratio at the steady state. At time t = 0 the value of G(t) was defined as unity, The QLV models presented here include an Ogden elastic representation and a Prony series relaxation function utilizing one and two terms, (where n = 1 and n = 2 specifically), that were referred to as "one-term QLV" and "two-term QLV", respectively. There are 2 independent viscoelastic parameters for the one-term model (τ 1 , G 1 ; note that G 1 = 1 − G 1 and therefore is not independent) and 4 for the two-term model (τ 1 , τ 2 , G 1 , G 1 ) as shown in Equation (7) and Here, in addition to one solid chain, models including one and two chains with viscous components are evaluated and denoted as one-term and two-term models. B: Illustration of how stimuli with a low strain-rate may lose information from low time-constant QLV chains. The solid line is the response of a typical two-term viscoelastic solid with time constants τ 1 and τ 2 under a step load. Two dashed lines represent the response of the material under slower strain-rates. For the slowest strain-rate, " _ 2 , stress relaxation properties of the faster chain (τ 1 ) may not show up because its relaxation for the faster chain takes place within its ramp phase. Therefore, this will not be captured by curve fitting, which simply characterizes the material as a one-term QLV solid and only calculates the slower chain (with time constant τ 2 ). Thus, we can eliminate the two extra parameters (G 1 , τ 1 ) if we only care about low strain-rate situations. In other words, for low strain-rate cases, the single term model is sufficient and therefore more appropriate than the two-term models because of the reduced number of free parameters. (8), not counting the two hyperelastic parameters (μ, α). Since the one-term model had fewer free variables (at the cost of inability to predict response for high strain-rates, Fig. 9B) it was therefore preferred given similar goodness of fit to the two-term model.
In contrast to the previous recursive method [40], the implemented numerical algorithm is designed to be in vector form so for-loops can be avoided, thus it is much easier to implement in numerical packages like MATLAB and NumPy. We first start from Equation (3), which can be split by applying Equation (5): Where By assuming no stress history before t = 0 and identifying the constants in Equation (10), we have And this is now ready for numerical implementation, as where the superscript k means the value at k th point in time, i.e., σ k means stress at time t k . The summation starts from l = 2 because we assert the stress change is zero at time zero. Also, from Equation (4) we calculate instantaneous stress as Similarly, we can obtain ðs l e À s lÀ1 e Þ: ð15Þ And the final stress evolution can be computed from Thus Equation (13)(14)(15)(16) completes the numerical implementation of QLV model.

Appendix II: Finite element analysis
To exclude the case where the frictional boundary conditions might confound the trend between thickness, stretch levels and viscoelasticity, numerical experiments of finite element (FE) analysis were performed using the commercial FE software package ABAQUS Standard, version 6.12 (Dassault Systèmes, Vélizy-Villacoublay, France). Numerical experiments of skin specimens were performed using axisymmetric models (Fig. 10), in which the geometry was the same as the biological specimen (6 mm dia. cylindrical skin piece). The material parameters assigned were from Table 1, using the two-term QLV model for higher accuracy. Three frictional coefficients (μ f ) between the skin and compression platen/table were tested, namely 0 (frictionless), 0.3 (between human finger and metal tip [41]), and 1 (rough). The rough friction coefficient also accounts for the cohesive force between skin and the metal, given that our boundary condition enforces no separation after contact [42]. Contact behaviors were defined in both a) tangential behavior, where isotropic friction was specified with penalty friction formulation, and b) normal behavior, where "hard" contact pressure-overclosure was used. A high Poisson's ratio ν = 0.475 for the skin were used. First, in line with the skin compression experiments, models of skin thickness from 200 to 800 μm with 100 μm increments was constructed. CAX4RH elements 50 μm in edge length were used. The stretch level and ramp time were both derived from the median values in the compression experiment where skin thickness varied, namely λ 1 = 0.63, t ramp = 0.129 s. The aforementioned model with skin thickness of 400 μm was then modified for a second experiment on variability in stretch level, where steady-state stretches were varied from 0.2 to 0.8 with an increment of 0.1. After the analyses were completed, the reaction force and displacement at the compression platen were extracted and processed in the same manner as the data from the skin compression experiments (described in Section 2.6), and the viscoelastic parameters were then compared to those obtained from the compression experiments.
In Fig. 11 we showed FE simulations on same skin thickness (400 μm) but under extreme frictional conditions (frictionless and rough), for different stimulus magnitudes (λ 1 from 0.5 to 0.7). Reduced relaxation functions obtained from fitting the hyper-viscoelastic constitutive model was also plotted in Fig. 11C. This shows that although changes in frictional conditions result in different force responses, the viscoelastic reduced relaxation functions are not impacted. The final outcome of the FE analysis, testing whether the effect of thickness on viscoelasticity is caused by frictional edge effects or the innate property of the skin, showed that the edge effect negligibly influences the outcome (Fig. 12A), independent of three levels of friction coefficients. Similarly, the frictional edge effects negligibly influence the outcome caused by strain level as well (Fig. 12B).
Appendix III: Role of dermis/epidermis thickness ratio In addition to the absolute thickness of a skin specimen, another independent variable that might contribute to viscoelastic variability is the relative thickness ratios between skin layers, as this changes between skin sites of Distal, OffNT and NT. The thickness ratios between dermis and epidermis, as previously obtained for different skin sites [12], were used as the independent variable here. Value of this dermis/epidermis thickness ratio, denoted as r, is listed in Table 2.
During the analysis of the first experiment when skin thickness varied, a multi-linear regression was performed with residual stress ratio G 1 as the dependent variable, and returned p < 0.001 for the independent variable thickness l 0 but p = 0.63 > 0.05 for the independent   Table 2. Value of dermis/epidermis thickness ratio r previously measured [12].