Thermally damaged porcine skin is not a surrogate mechanical model of human skin

Porcine skin is considered a de facto surrogate for human skin. However, this study shows that the mechanical characteristics of full thickness burned human skin are different from those of porcine skin. The study relies on five mechanical properties obtained from uniaxial tensile tests at loading rates relevant to surgery: two parameters of the Veronda-Westmann hyperelastic material model, ultimate tensile stress, ultimate tensile strain, and toughness of the skin samples. Univariate statistical analyses show that human and porcine skin properties are dissimilar (p < 0.01) for each loading rate. Multivariate classification involving the five mechanical properties using logistic regression can successfully separate the two skin types with a classification accuracy exceeding 95% for each loading rate individually as well as combined. The findings of this study are expected to guide the development of effective training protocols and high-fidelity simulators to train burn care providers.


Scientific Reports
| (2022) 12:4565 | https://doi.org/10.1038/s41598-022-08551-z www.nature.com/scientificreports/ skin tissues are ~ 111 MPa after freezing to − 80 °C and at 100% relative humidity. The change in temperature due to freezing is found to affect the Young's modulus of stratum corneum and full thickness of human and porcine skin tissues in the exact opposite manner 15 . However, these comparisons are limited to unburned skin tissues. Further, there is no study that performed a statistical comparison of the mechanical properties of burned human and porcine skin tissues. The hypertrophic scar tissue from a burn injury is stiffer than normal human skin tissue 7,16 . Similarly, the mechanical properties of burned and unburned porcine skin tissues are found to be significantly different 17 . Porcine skin is known to soften with an increase in temperature along with a change in the transition region of the typical J-curve response of skin tissues. This softening behavior is more prominent at higher strain values 4 . Hardening of the porcine skin is also observed with an increase in strain rate 4 . Hence, it is important to quantify whether full thickness burned porcine skin can be considered as a surrogate for full thickness burned human skin at various strain rates.
To the best of our knowledge, this is the first paper to study statistical similarity between the mechanical properties of full thickness burned ex vivo porcine and human skin tissues. Five mechanical properties, namely, ultimate tensile stress, ultimate tensile strain, toughness, and two parameters of the Veronda-Westmann hyperelastic material model 18 , are calculated from stress-strain data collected from uniaxial tensile tests on full thickness burned human and porcine skin tissues. The uniaxial tests are performed at three loading rates: 0.3 mm/s, 2.0 mm/s, and 8.0 mm/s. The loading rate of 0.3 mm/s corresponds to a quasi-static loading condition, whereas the other two loading rates are consistent with the cutting speeds observed during skin surgery 19 . Univariate hypothesis tests are performed on each of the five mechanical properties to study the statistical similarity between the two types of tissues. In addition, multivariate statistical analysis is carried out to differentiate between two types of skin tissues based on the five mechanical properties.
The paper is organized as follows. The experimental setup for the uniaxial tensile tests and the statistical analysis on the material properties are described in "Materials and methods" section. The results of the univariate and multivariate statistical analysis are presented in "Results" section along with a discussion of the results in "Discussion" section. The summary of the study is presented in "Conclusion" section, along with future directions. The description of the human skin samples and the comparison of various material models used to fit the experimental stress-strain data are given in Appendix A in supplementary information. The details of the univariate and multivariate statistical analysis are provided in Appendix B in supplementary information.

Materials and methods
The uniaxial tensile tests are conducted on the full thickness burned ex vivo porcine skin and debrided human skin tissues at three different loading rates, that is, 0.3 mm/s, 2 mm/s, and 8 mm/s, until the rupture of the specimen. Force and displacement data collected from these experiments are used to calculate ultimate tensile stress, ultimate tensile strain, toughness, and the parameters of Veronda-Westmann hyperelastic material model 18 . These five material properties are used to compare the mechanical characteristics of porcine and human skin tissues.
The experimental protocol for the uniaxial tensile testing of full thickness burned human and porcine skin samples are presented in "Human sample testing protocol" and "Porcine sample testing protocol" sections, respectively. The calculation of the material properties from the data collected from the uniaxial tensile tests is given in "Characterization of material properties" section. The methodologies for the univariate and multivariate statistical analysis are given in "Statistical analysis" section.
Human sample testing protocol. The human skin samples were tested at the University of Washington in Seattle, WA, and de-identified measurement data were provided for this study. The human samples were collected as debrided/discarded patient tissues that had undergone full thickness or deep partial burns, at the Harborview Medical Center Burn Unit in Seattle, Washington, in accordance with the protocol approved by the University of Washington Institutional Review Board (IRB). The informed consent requirement is waived by the IRB committee for the access and use of medical records and debrided/discarded skin tissue specimens. Deep partial burns can be classified as full thickness burns for this study as the deep partial burns can develop eschar, which requires escharotomies just as full thickness burns do. The anatomical location of the burned skin and the number of post-burn hours after the tissue had been excised were recorded, along with the demographics of each patient such as age and gender. After the burned skin tissue was excised, it was submerged in 1X phosphate-buffered saline (PBS) (pH 7.4) to prevent loss of hydration and then stored in refrigeration. The tissue sample was transferred to the University of Washington. The subcutaneous tissue was removed, and the tissue was processed into standard dog-bone shape specimens using an ASTM D638 Type V die, and the gauge length was recorded using tissue dye and a stencil. The uniaxial tensile tests were performed on a TA Instruments ElectroForce Test-Bench (TA Instruments, MN) instrument in a 37 °C PBS bath at loading rates of 0.3 mm/s, 2 mm/s, and 8 mm/s until failure. The loading rate of 0.3 mm/s corresponds to a quasi-static loading condition, whereas 2 mm/s and 8 mm/s are consistent with the cutting speeds observed during skin surgery 19 . The samples were preloaded to 0.10 N and the cross-sectional area was measured using laser micrometers IG-028 (Keyence, IL). Force and displacement data are collected from the testing machine. The material characterization study described above is performed in accordance with the guidelines and regulations of the IRB protocol.
For the loading rates of 0.3 mm/s, 2 mm/s, and 8 mm/s a total of 107, 105, and 108 burned skin samples excised from 14, 14, and 15 human subjects, respectively, are used in this study. The samples were obtained from 13 male and 2 female subjects in the age range 38 to 78 years. The time between the burn injury and the incision of the skin sample varies between 23 and 68 hours. In four extreme cases, skin tissues were excised after 11, 120, 139, and 360 hours. The mean thickness of the human skin samples is 2.1 ± 1.0 mm. The anatomical location of the burn injury for the 15 subjects is given in supplementary information (Appendix A). in this study. The tissue was stored in refrigeration before the skin was separated from the subcutaneous and underlying muscles using a scalpel. A commercial griller (Cuisinart ® GR-300WS Griddler Elite Grill, Conair Corporation, NJ) was set to 450 °F, which was confirmed using an infrared camera (InfraCam, FLIR, OR), and the tissue was burned for 30 s to induce full thickness burns. This burn condition is chosen as a surrogate of burn depth, which is known to result in full thickness burn 17,20 . The tissue was submerged in 1X PBS solution until it reached room temperature, and it was then processed into standard dog-bone shape specimens using an ASTM D568 Type V die. This was done to avoid shrinkage of skin tissue during the contact heating process. The uniaxial tensile tests were performed on an Instron ® MTS system (INSTRON, MA) at the loading rates of 0.3 mm/s, 2 mm/s, and 8 mm/s until failure. Force and displacement data are collected from the testing machine. A total of 38, 39, and 39 samples were used for testing at three loading rates, respectively. The mean thickness of the porcine skin samples is 1.5 ± 0.7 mm.
Characterization of material properties. The nominal stress σ N = F/A 0 and nominal strain ε N = �L/L 0 are computed using force F and displacement L measurements obtained from the uniaxial tensile tests of the dog bone samples. Here, A 0 is the initial cross-sectional area of the gauge section and L 0 is the initial length of the sample, as shown in Fig. 1. The ultimate tensile stress, ultimate tensile strain, and toughness were calculated from the nominal stress-strain curve, as shown in Fig. 1. The ultimate tensile stress describes the maximum nominal stress that the sample can withstand before it ruptures under uniaxial tensile loading. The ultimate tensile strain corresponds to the nominal strain at the ultimate tensile stress. Toughness is calculated as the area under the nominal stress-strain curve, i.e., the amount of energy per unit volume that a sample can absorb before rupturing. The stress-strain constitutive response of the full thickness burned human and porcine skin tissues is described by the Veronda-Westmann hyperelastic material law, with the strain energy density function given as where µ and γ are two material coefficients, I 1 and I 2 are the first and second deviatoric strain invariants, respectively. For the uniaxial tensile test, the nominal stress for the Veronda-Westmann material is then given as where U = 1 + ε N is the principle stretch along the uniaxial loading direction and ε N is the nominal strain under uniaxial loading conditions. The parameters, μ and γ, are calculated by nonlinear least squares fit of the  Statistical analysis. Univariate and multivariate statistical analysis is carried out on the five material properties of the full thickness burned human and porcine skin tissues, viz., ultimate tensile stress, ultimate tensile strain, toughness, and the parameters of the Veronda-Westmann hyperelastic material model, μ and γ. The methodology for the univariate and multivariate statistical analyses is presented in "Univariate statistical analysis" and "Multivariate statistical analysis", respectively.
Univariate statistical analysis. The distributions of the five material properties of the full thickness burned human and porcine skin tissue samples are compared between the two types of tissues at each of the three loading rates, described in the following steps.
1. First, the Shapiro-Wilk normality test 21 is performed on the sample data of each material property for human and porcine skin tissues and each loading rate. 2. If we can assume that a pair of sample data (human and porcine for a particular parameter and loading rate) was drawn from normal distributions: a. The two-sample F-test for equal variances 22 is carried out. b. Depending on whether the variances of the two samples are equal or unequal, a two-sample t-test for equal or unequal variances 22 is conducted, respectively.
3. If we cannot assume that a sample pair was drawn from two normal distributions, a. the two-sample Kolmogorov-Smirnov test 23 is performed. b. If the two samples were drawn from continuous distributions that have the same shape, the Wilcoxon rank-sum test 22 is performed. c. If the two samples were not drawn from continuous distributions that have the same shape, a two-sample t-test for unequal variances is performed, as per the recommendation in Ruxton 24 .
The significance of all the hypothesis tests is chosen to be 0.01. A flowchart of the univariate hypothesis testing and the details of the null hypothesis for each test used in the study are defined in supplementary information (Appendix B).
For a post-hoc analysis to verify that the sample sizes exceeded the required minimum, the effect size for each sample pair is estimated by computing Cohen's d statistic. This statistic is a numerical representation of the strength of the relationship between two variables in a population. An effect size of d = 0.2, is considered small, d = 0.5 is considered medium, and d = 0.8 is considered large 25 . The sample size estimations are calculated using G*Power 3.1.9.7 26 , with the significance level and power chosen to be 0.01 and 0.99, respectively.
Multivariate statistical analysis. Although univariate hypothesis testing can highlight differences in the mean or median between two classes, it cannot qualitatively describe the degree of such differences. However, this is essential since simply rejecting the null hypothesis does not offer an insight into the significance of the differences between the mean or median of the classes. Hence, a multivariate classification analysis 27 is carried out by combining the five variables: ultimate tensile stress, ultimate tensile strain, toughness, and the parameters of the Veronda-Westmann hyperelastic material model, μ and γ.
Binary classification is carried out for each loading rate of 0.3 mm/s, 2 mm/s, and 8 mm/s based on pairwise samples (one for burned human and one for porcine specimen). The classification is performed using logistic regression 28 . A leave-one-out cross-validation is applied in order to guarantee an independent assessment of the classifiers. This is a holdout strategy, where one data point is removed from training the classifier model and later used to assess the classification accuracy of the trained model. The procedure is concluded when each data point is removed once, starting with the first and finishing with the last data point. The smaller the misclassification error, i.e., the percentage of incorrectly classified observations, the more significant is the difference between the two samples.
In addition to a binary classification based on the three different loading rates, the three datasets for human, and the three datasets for porcine are also combined. This yields two data sets that allows conducting a direct, rate independent, binary classification to assess the difference in the material properties between full thickness burned human and porcine skin tissues.

Results
The univariate and multivariate statistical analysis of the five material properties, i.e., the ultimate tensile stress, ultimate tensile strain, toughness, and the parameters of the Veronda-Westmann material model, of full thickness burned human and porcine skin is carried out to compare the behavior of these skin types at different loading

Univariate statistical analysis. The samples of the full thickness burned human and porcine skin tissues
with at least one of the five properties greater than three times the interquartile range of the distribution of the property were removed from the analysis. The number of samples used in the statistical analysis after removing the outliers for human and porcine skin tissues for the 0.3 mm/s, 2.0 mm/s, and 8.0 mm/s loading rates is given in Table 1. The Cohen's d effect size for the human and porcine skin tissue samples based on the five material properties is found to be 1.4. An allocation ratio of 2.5 is chosen, which is close to that used in this study. The sample size estimation for human and porcine tissues based on the Cohen's d value and the allocation ratio is 46 and 18, respectively. This indicates that the sample size used in this study is appropriate for the univariate hypothesis tests. The box plots of the five material properties of the human and porcine skin tissues at the three loading rates are shown in Fig. 2.
The results of the univariate hypothesis tests are detailed in supplementary information (Appendix B). The hypothesis tests show that the means or medians of ultimate tensile stress, toughness, and the two parameters, μ and γ, of the Veronda-Westmann model are not equal for all loading rates between the human and porcine skin tissue samples. The difference in the median and variance of the properties between the human and porcine skin tissues at all loading rates is also observed in Fig. 2. Based on the Kolmogorov-Smirnov test, we reject the claim that the shape of the two distributions is identical for the two tissues for the ultimate tensile stress at loading rate of 8 mm/s and for the two parameters, μ and γ of the Veronda-Westmann model at all three loading rates. Although, the median of the ultimate tensile strain of human and porcine skin tissues is equal for the loading rate of 0.3 mm/s based on the Wilcoxon rank-sum test, it is not equal for the loading rates of 2 mm/s and 8 mm/s.  www.nature.com/scientificreports/ Multivariate statistical analysis. The multivariate analysis using logistic regression statistical model is carried out to differentiate between full thickness burned human and porcine skin tissues based on the five material properties, namely, the ultimate tensile stress, the ultimate tensile strain, toughness, and the parameters μ and γ of the Veronda-Westmann material law. Binary classification is performed to separate human and porcine skin tissues for each of the loading rates, i.e., 0.3 mm/s, 2 mm/s, and 8 mm/s. In addition, human and porcine skin samples are classified by combining samples from all three loading rates into respective classes. Table 2 summarizes the performance metrics, such as traditional accuracy, sensitivity, specificity, F1-score, area under the receiver operating characteristic curve (ROC-AUC), as well as Matthew's correlation coefficient (MCC) 29 , Fowlkes-Mallow's index (FMI) 30 , and Adjusted Rand index (ARI) 31 , computed from the confusion matrix obtained from the leave-one-out cross-validation. The confusion matrices for the binary classifications are provided in supplementary information (Appendix B). The human and porcine skin types can be separated with a high degree of accuracy (> 0.95), sensitivity (> 0.96), and specificity (> 0.92) for all loading rates (Table 2). This can also be confirmed by the high value of the F1-score, as well as the MCC, FMI, and ARI metrics further indicates that the two tissue types can be accurately classified into their respective classes based on five mechanical properties, irrespective of loading conditions. A factor analysis of the classifier is finally performed to quantify the contribution of each mechanical property in differentiating between the two types of skin tissues. This is to examine how influential each of the five material properties is in differentiating between the two skin types. Figure 3 shows that the parameter μ of the material model and the ultimate tensile strain, respectively, are the most and the least discriminative factors among the five properties for all loading conditions. The material model parameter γ only contributes to differentiating the two skin types for loading rate 2 mm/s.

Discussion
The univariate and multivariate statistical analyses using five mechanical properties, i.e., ultimate tensile stress, ultimate tensile strain, toughness, and two parameters of the Veronda-Westmann material model, reveal that the mechanical characteristic of full thickness burned human skin is significantly different from that of porcine skin tissues, irrespective of the three loading conditions tested in this study. These findings are consistent with previous studies that indicate that the tangent modulus of human skin is different from the porcine skin tissues, even at similar anatomical locations 13 . However, it contradicts the assumption that porcine skin is a de facto surrogate for human skin to evaluate the mechanical characteristics 3,4,9 . The difference in the mechanical behavior of the full thickness burned human and porcine skin tissues can be attributed to the structural, compositional, and functional differences between these tissue types. The elastic fiber content in dermis is different between the two types of skin tissue 8 . Phenomenological models further show that the elastic modulus of the skin of animals and humans is a function of the thickness of each skin layer, which are not similar between humans and porcine 32 . In addition, the volume fraction of moisture content in healthy porcine skin tissue is more than that of human skin. The corneocytes of the stratum corneum retain moisture when skin tissue is not hydrated, while the corneocytes and lipids in the stratum corneum swell with hydration and thereby significantly alter the tangent modulus of the tissue. These layers of corneocyte and lipid content are different between human and porcine tissues 15 . The hydrogen bonding states of water and the natural moisturizing factor concentration is also lower in porcine skin compared to human skin which affects the permeability of the two skin types 33 . The natural moisturizing factor is composed of water-soluble compounds such as free amino acids and its derivatives, inorganic salts, sugars, lactic acid, urea etc. 34 . Further, the water loss in a healthy pig is mainly through rapid evaporation and it is found to be similar in in vivo and ex vivo porcine skin tissues. Whereas the water loss in human skin at higher temperatures is significantly more than that in porcine skin 35 . Hence, it is expected that the water content in a full thickness burned porcine skin is significantly different from that of burned human skin. The viscoelastic behavior of skin tissue depends on the interaction between collagen, proteoglycans, and water molecules 4 . A difference in the water content between the two types of burned skin severely affects this interaction among the constituents of skin tissue which results in the dissimilar mechanical properties.
Further, the lipid content and keratin type and its structure in the stratum corneum of the unburned skin tissue is different between human and porcine 33 . It is observed that the disruption of the lipid and keratin structures in the stratum corneum causes significant changes to the permeability and mechanical response of burned skin tissue 36 . Hence, it is expected that the difference in the lipids and keratin content and structure between burned human and porcine skin tissue may lead to a significant difference in mechanical properties, as observed in this study.
The effect of the anatomical location of the full thickness burned skin tissue on the mechanical properties could not be studied due to the insufficient human skin tissue samples available for different anatomical locations. It is unlikely that the difference in mechanical characteristics of human and porcine skin stems from the difference in the experimental conditions and anatomical location of the tissue samples. The mechanical properties of skin tissues are known to change with temperature due to collagen denaturation 4 . Collagen and elastin content contribute to the elasticity of the skin tissue 1 . Collagen fiber orientation is responsible for the variation in elasticity and strength of the skin tissue along different orientations at various anatomical locations of the skin 13 . However, when the skin tissue undergoes full thickness burn, the collagen fibers undergo denaturation 37 . The fibers transform from an organized crystalline structure to a random gel-like state 38 . Hence, the difference in burn temperature and collagen fiber orientation in samples from anatomically different locations is unlikely to contribute to the statistically significant difference in the mechanical properties of the two skin tissue types.
In the study, in vivo for human skin samples are used, whereas it is ex vivo for the porcine skin samples. However, both human and porcine skin tissue samples have undergone full thickness burn. For the human skin tissue, the burn condition is confirmed by experienced burn care personnel at the Harborview Medical Center Burn Unit, Seattle, WA. For the porcine skin tissue, the contact burn at 450 °F and for 30 s is known to cause full thickness burn 17,20 . Hence, we expect similar burn-related thermal damage in both tissue types. Similar comparisons of in vivo human skin and ex vivo porcine skin have been studied in the literature 33 . It should also be noted that mechanical properties of ex vivo porcine skin at high temperatures are studied in literature to understand the mechanical behavior of human skin at high temperatures 4 .
The conclusions of the present study are limited to the mechanical characteristics of full thickness burned human and porcine skin. Further research needs to be conducted to explore the differences between human and porcine tissues for other properties under thermal damage.

Conclusion
Uniaxial tensile tests are conducted on full thickness burned ex vivo human and porcine skin tissues at three loading rates of 0.3 mm/s relating to quasi-static loading, 2 mm/s, and 8 mm/s relating to the rates in skin surgery. The ultimate tensile stress, ultimate tensile strain, toughness, and the two parameters of the Veronda-Westmann material model are calculated for the human and porcine samples at the three loading rates using the force-displacement data from the tensile tests. Univariate hypothesis tests indicate a significant difference in all five material properties between the two types of skin tissues for all loading rates. The classification accuracy from the multivariate analysis between the human and porcine skin tissues is more than 95% when considering samples from each loading rate separately and together using logistic regression statistical model. Hence, the full thickness burned porcine skin tissue is statistically not similar to that of the human skin tissue. This observation is significant for the development of training protocols and high-fidelity simulators to train burn care providers. www.nature.com/scientificreports/ This study could be further extended to investigate the effect of loading rates on the rate-dependent viscoelastic response of the severely burned skin tissues.

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.