Three-Dimensional Quantification of Collagen Microstructure During Tensile Mechanical Loading of Skin

Skin is a heterogeneous tissue that can undergo substantial structural and functional changes with age, disease, or following injury. Understanding how these changes impact the mechanical properties of skin requires three-dimensional (3D) quantification of the tissue microstructure and its kinematics. The goal of this study was to quantify these structure-function relationships via second harmonic generation (SHG) microscopy of mouse skin under tensile mechanical loading. Tissue deformation at the macro- and micro-scale was quantified, and a substantial decrease in tissue volume and a large Poisson’s ratio was detected with stretch, indicating the skin differs substantially from the hyperelastic material models historically used to explain its behavior. Additionally, the relative amount of measured strain did not significantly change between length scales, suggesting that the collagen fiber network is uniformly distributing applied strains. Analysis of undeformed collagen fiber organization and volume fraction revealed a length scale dependency for both metrics. 3D analysis of SHG volumes also showed that collagen fiber alignment increased in the direction of stretch, but fiber volume fraction did not change. Interestingly, 3D fiber kinematics was found to have a non-affine relationship with tissue deformation, and an affine transformation of the micro-scale fiber network overestimates the amount of fiber realignment. This result, along with the other outcomes, highlights the importance of accurate, scale-matched 3D experimental measurements when developing multi-scale models of skin mechanical function.


INTRODUCTION
Skin is a complex organ that provides a protective external barrier for the underlying musculoskeletal tissues of the human body during daily living. Understanding the biomechanics of skin is important for preventing injury, planning clinical procedures, engineering wearable technologies, improving skin care products, and understanding aging (Aubert et al., 1985;Sen et al., 2009;Hussain et al., 2013;Alberto et al., 2018;Blair et al., 2020). The mechanical response of skin has been extensively studied in the past. Through studies involving mechanical testing and histology, collagen fibers within the dermis have been identified as the primary load bearing structure in skin (Fackler et al., 1981;Blair et al., 2020). Additionally, the relative organization of collagen fibers plays a key role its anisotropic mechanical behavior (Billiar and Sacks, 1997;Lake et al., 2012). To better understand the relationship between collagen fiber organization and the mechanical function of soft tissues, substantial advancements involving computational modeling that incorporate the collagen microstructure have been made over the last 20 years (Fung, 1967;Sander et al., 2009;Fang and Lake, 2016;Limbert, 2017). The development of these computational models has outpaced experimental measurements of the collagen microstructure under mechanical loading, and many assumptions regarding the fiber and tissue kinematics at the microscale remain untested, which limits our understanding of tissue structure-function relationships (Fang and Lake, 2016).
Historically, histological staining of tissue has been the gold standard for evaluating skin microstructure and its relationship to biomechanics (Ridge and Wright, 1966;Lanir, 1979;Jor et al., 2013;Pensalfini et al., 2018). However, this assessment of microstructure is inherently destructive due to the sample preparation process and cannot be used to track tissue or fiber kinematics during mechanical loading. Digital image correlation (DIC) has previously been used in various studies of skin to observe interesting macro-scale tissue kinematics, such as a large Poisson's ratio and a decrease in volume during uniaxial mechanical stretching, but these techniques cannot be used to quantify fiber kinematics (Hepworth et al., 2001;Pissarenko et al., 2019;Wahlsten et al., 2019;Pissarenko and Meyers, 2020). To understand the relationship between collagen fiber kinematics and the tissue's mechanical response, imaging methods such as quantitative polarized light imaging (QPLI) and small angle light scattering (SALS) have been integrated with mechanical testing Tower et al., 2002;Quinn and Winkelstein, 2008;Sander et al., 2009;Woessner et al., 2019). Instead of visualizing individual fibers, however, these methods infer the average orientation and degree of alignment of thousands to millions of fibrils with macro-scale resolutions on the order of millimeters. These approaches provide a dynamic two-dimensional (2D) summary of fiber kinematics and organization during mechanical loading. However, collagen fibers are organized in a threedimensional (3D) matrix that can vary as a function of depth, and skin and other soft tissues have differences in organization among all three-dimensions (Roeder et al., 2002;Limbert, 2017).
To better understand the reorganization of collagen fibers during mechanical loading, 2D distributions of the micro-scale collagen fiber reorganization in human tendon (Lake et al., 2012), bovine pericardia (Sacks, 2003), and porcine mitral valve leaflets (Lee et al., 2015) have been found to share an affine relationship with tissue deformation. When similar experiments were performed on porcine (Hepworth et al., 2001) and mouse (Bancelin et al., 2015) skin, the 2D micro-scale distribution of collagen fiber reorientation did not match the distributions one would predict if the fiber kinematics were affine with the macro-scale tissue deformation. Furthermore, previous 3D measurements of multi-scale strain in collagen gels revealed that the magnitude of micro-scale strain in the direction of stretch was lower than the amount of applied strain at the macro-scale, indicating a need to quantify both 3D fiber kinematics and organization at the micro-scale in skin (Roeder et al., 2009). In summary, a more detailed 3D description of the multi-scale mechanical behavior of skin, particularly with respect to its microstructure, is still needed to inform, validate, and refine appropriate computational and constitutive models.
Depth-resolved optical microscopy approaches offer both 3D imaging capabilities and sub-micrometer resolution capable of non-destructive visualization of individual fibers within a tissue. One such imaging approach, reflectance confocal microscopy, has been used to characterize multiscale deformation and fiber kinematics in 3D collagen gels during mechanical stretch, showing that the relative magnitude of micro-scale strain can be quite different from the amount of strain experienced at the bulk tissue scale (Roeder et al., 2002(Roeder et al., , 2004(Roeder et al., , 2009Voytik-Harbin et al., 2003). However, soft biological tissues are more complex than collagen gels, with more heterogeneity in matrix structure and composition and a greater density of extracellular matrix (ECM) proteins. Due to the high density of matrix components in native tissues, reflectance confocal microscopy has limited depth penetration and lacks the ability to discriminate compositional elements in skin. Multiphoton microscopy (MPM) is a non-linear optical microscopy technique that has greater depth penetration than confocal microscopy due to the use of near-infrared light (Yasui et al., 2013). Through MPM, collagen fibers can be visualized at the micro-scale with high specificity based on a strong second harmonic generation (SHG) signal, and two-photon excited fluorescence (TPEF) from endogenous cellular fluorophores and other ECM proteins can provide additional context during MPM imaging (Quinn et al., 2016). Recently, SHG has been utilized to measure 2D micro-scale fiber kinematics and organization of mouse skin during mechanical testing, but these experiments were performed only on the dermal layer of the skin, and have yet to characterize fiber kinematics and tissue deformation in 3D (Bancelin et al., 2015;Lynch et al., 2017a,b).
The goal of this study was to quantify 3D fiber organization and kinematics at the micro-scale, as well as characterize tissue kinematics across multiple length-scales during tensile loading of skin. To this end, a uniaxial mechanical testing device was integrated with a multiphoton microscope to quantitatively assess macro-and micro-scale kinematics simultaneously during mechanical stretch. 3D analysis of micro-scale fiber kinematics and organization was performed using established analysis techniques. Macro-scale kinematics were calculated through analysis of fiduciary markers and tissue kinematics were compared across length-scales. This multi-scale analysis of skin provides new insights into how the 3D microstructure responds to tensile loading, which is critical for developing and refining computational models that relate tissue structure and mechanical function.

Sample Preparation
Full thickness ventral skin from male C57BL/6J mice (4 months old, n = 7) was collected, snap-frozen, and stored at −80 • C until mechanical testing. Prior to testing, tissue samples were thawed and prepared by cleaning the skin with an alcohol wipe and carefully resecting the hypodermis. Traditional dog-bone shapes were cut from the skin in the medial-lateral direction using a precision-machined steel stamp (26 mm × 8 mm total shape size, 6 mm × 2 mm gauge region) (Sander et al., 2014), and four fiduciary markers were added to the tissue gauge region using a fine-tipped ink pen to enable macro-scale analysis of tissue kinematics.

Data Acquisition
Tissue samples loaded in the uniaxial mechanical testing device were submerged in a bath filled with phosphate buffered saline (PBS) to ensure hydration during the duration of the test. Prior to loading, the tissue geometry was optically measured using the macro-scale imaging system and pre-stressed to 10 kPa, which was considered the undeformed configuration for all samples (Quinn and Winkelstein, 2008). Tissue samples were preconditioned to ∼5% strain (20 cycles, 4 mm/min) to ensure a repeatable mechanical response (Giles et al., 2007;Quinn and Winkelstein, 2011), and the allowed to rest for 5 min to allow for hydrostatic recovery. Representative SHG and TPEF image volumes at three discrete locations spanning approximately the length of the gauge region were then collected to assess the undeformed microstructure. During mechanical testing, tissues samples were stretched in 1 mm increments at a quasi-static rate of 0.2 mm/min (Figure 2A). The initially undeformed ROI taken within the center of the gauge region was visually tracked during mechanical stretching via SHG imaging (∼1.2 Hz), and at every 1 mm increment, the mechanical stretching was stopped to allow for collection of a high-resolution SHG image z-stack of the ROI (Bancelin et al., 2015).

Analysis of Deformation
At the macro-scale, fiduciary markers were tracked within the camera images using an automated normalized cross-correlation algorithm, and the deformation gradient (F) was computed at every 1 mm increment of displacement using isoparametric mapping (Quinn and Winkelstein, 2009). The deformation  To enable the collection of high-resolution micro-scale SHG image volumes, the skin was stretched in 1 mm increments (red circles) until mechanical failure. (B) Tissue kinematics were restricted to 2D at the macro-scale but were quantified in 3D at the micro-scale. (C) Macro-scale tissue deformation measured from fiduciary markers showed a substantial compressive principal strain orthogonal to the direction of stretch (scale bar represents 2.5 mm). (D) 3D micro-scale tissue kinematics measured from SHG features indicated compressive principal strains were occurring in both directions orthogonal to the loading direction (scale bar represents 100 µm). For panels (C,D), white polygons indicate tracked features, and arrows represent the direction and magnitude of principal strains. gradient was then used to calculate the 2D Green strain tensor (E) (Figure 2B), where I is the identity tensor. Principal strain magnitudes and directions were computed as the eigenvalues and eigenvectors of the 2D strain tensor, respectively ( Figure 2C; Quinn and Winkelstein, 2008;Oñate, 2009;Nikishkov, 2010). Using the 2D strain tensor measurements, the incremental Poisson's ratio (ν), was calculated for each tissue sample during loading. Due to the Poisson's ratio being a material property, all incremental Poisson's ratio measurements for a particular sample were averaged together.
To calculate micro-scale tissue kinematics, trackable collagen fiber patterns were first identified using a 3D extension of the Harris corner detection algorithm (Harris and Stephens, 1988). Pixels surrounding each fiber location selected by the corner detection algorithm were stored (21 × 21 × 11 pixels), and at every 1 mm increment of displacement, locations were tracked in 3D. Accurate and precise tracking of 3D features consisted of first computing the position of current features (x t , y t , z t ) in the next sequential image volume (x t+1 , y t+1 , z t+1 ) using a 3D normalized cross-correlation algorithm. Pixels surrounding the (x t+1 , y t+1 , z t+1 ) position were then collected and used to recalculate the feature position in the current image volume (x t ', y t ', z t '). Feature displacement was then calculated as the average displacement between the (x t+1 , y t+1 , z t +1 ) and (x t ', y t ', z t ') vectors, and the amount of tracking error was calculated as the Euclidean distance between the (x t , y t , z t ) and (x t ', y t ', z t ') vectors. To filter out features that did not correctly track among all image volumes, only features with a maximum tracking error of <5 pixels for the entire mechanical test were considered for analysis (Quinn and Winkelstein, 2010). Since at least four features are required to calculate micro-scale deformation in 3D, only image volumes that contained at least four correctly tracked features were used for analysis. To get a 3D measurement of micro-scale deformation, 3D tetrahedrons were generated from all possible combinations of correctly tracked features, and their volumes were calculated in the undeformed state. The four feature locations that made up the largest possible tetrahedron based on volume were used to calculate micro-scale tissue kinematics ( Figure 2D). Similar to the macro-scale analysis, the deformation gradient at each 1 mm increment was calculated using 3D isoparametric mapping (Nikishkov, 2010), and a 3D Green strain tensor was calculated at each 1 mm increment ( Figure 2B). Principal strain magnitudes and directions were computed as the eigenvalues and eigenvectors of the Green strain tensor, respectively ( Figure 2D). To quantify how the micro-scale volume changed as a function of stretch, a volume ratio was calculated from the deformation tensor at each 1 mm of displacement,

Analysis of Fiber Organization and Volume Fraction
Three-dimensional collagen fiber orientation for SHG image volumes was measured using a weighted vector summation algorithm as previously described (Quinn and Georgakoudi, 2013;Liu et al., 2017). Calculation of pixel-wise 3D fiber orientation was based on a surrounding 7 × 7 × 7 pixel window and resulted in a pixel-wise map of spherical coordinates, θ (X-Y plane orientation) and ϕ (Z axis inclination) (Figures 3A,B). Fibers within the image voxel were then isolated by first computing an average SHG intensity for the entire image volume, and all pixels with intensity greater than or equal to the average voxel intensity were used to create a collagen-positive volumetric mask. To assess fiber alignment during stretch, directional variance based on the 3D orientation of all fiber-containing pixels within an image volume was calculated. Directional variance is bounded by values of zero and one, which correspond to an aligned and random fiber network organization, respectively ( Figure 3C; Liu et al., 2017). Additionally, a fiber volume fraction was calculated as the ratio between the number of collagen-positive pixels and the total number of pixels within an image volume.
To assess whether fiber alignment during stretch was affine with the local tissue deformation, comparisons were made between the measured fiber directional variance at peak strain and the theoretical affine fiber kinematics at peak strain. The measured 3D deformation gradient tensor was applied to the undeformed image volume using builtin MATLAB functions (affine3D and imwarp) to produce theoretical image volumes at peak strain under an affine assumption (Fan and Sacks, 2014;Jayyosi et al., 2017). As was the case with the actual fiber kinematics, maps of 3D fiber orientation were computed from the theoretical affine image volumes. To compare fiber organization for both measured and theoretical deformed image volumes, 3D directional variance was calculated for all fiber-containing pixels within the 3D tetrahedron formed from the tracked features. Additionally, 2D directional variance was calculated for both fiber orientation (θ) and inclination (ϕ) (Figures 3A,B) to determine whether the fiber kinematic response was affine with the local deformation in all planes.
Collagen fiber directional variance has previously been found to vary as a function of the length scale over FIGURE 3 | SHG image volumes were used to compute a pixel-wise map of 3D micro-scale orientation (A) and inclination (B). The two 3D orientation and inclination maps and collagen-positive mask were then convolved with a range of cubic averaging kernels varying in size to determine how regional measurements of directional variance (C) and fiber volume fraction (D) change as a function of the length scale over which they are measured (convolution with 40 µm 3 kernel shown here, scale bar represents 100 µm).
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org which it is measured (Quinn et al., 2015). For this study, fiber directional variance and volume fraction at the microscale was measured at multiple length scales within the image volume by convolving the pixel-wise fiber orientation and inclination maps, as well as the fiber-containing mask with cubic averaging kernels of varying sizes (4 3 , 10 3 , 20 3 , 40 3 , 70 3 , and 100 3 µm 3 ). This convolution process results in a regional summary of directional variance and fiber volume fraction, which was then averaged and compared among length-scales (Figures 3C,D

Statistical Analysis
At both the macro-and micro-scale, the tissue kinematics were compared to the expected outcomes of an incompressible material (ν = 0.5 for Poisson's ratio, V/V 0 = 1 for volume ratio) through a one sample t-test. Comparisons of strain across the two length scales, as well as between undeformed and deformed (at peak strain) micro-scale organization and fiber volume fraction were performed using paired student's t-tests. Additionally, comparisons of directional variance between theoretically and mechanically deformed image volumes were performed using paired student's t-tests. To assess the length scale dependencies at the micro-scale, comparisons were performed with one-way ANOVAs and post hoc Tukey's HSD tests to determine significant differences. Significant differences were defined as a p-value < 0.05, and all statistical analysis was done using JMP Pro 13 (SAS Institute; Cary, NC, United States).

Tissue Deformation
At the macro-scale, skin under uniaxial tension (in the x-direction) exhibited a significant inward orthogonal contraction (in the y-direction) (Figure 2C), producing an average Poisson's ratio (1.57 ± 0.65) that was significantly greater than that of an incompressible material (p = 0.004) (Figure 4A). When compared across length-scales, micro-and macro-scale strain measurements did not significantly differ in either the E xx or E yy directions (p = 0.26 and p = 0.44, respectively) ( Figure 4B). Interestingly, the 3D deformation gradient measured via MPM image feature tracking revealed a substantial decrease in tissue volume ratio (V/V 0 ) with increasing strain that was significantly different from the expected response of an incompressible material (p = 0.0003) or other rubber-like material models historically used to model soft biological tissue (Figure 4C; Fung, 1967;Crichton et al., 2013).

Micro-Scale Fiber Organization and Volume Fraction
The volume fraction and alignment of collagen fibers can provide some context for the large changes in tissue volume observed during loading. At the micro-scale, initial collagen fiber directional variance for the entire image volume (0.654 ± 0.107) was significantly higher than at peak strain (0.524 ± 0.100; p = 0.002), indicating increased fiber alignment with stretch ( Figure 5A). Despite significant reductions in the macro-scale tissue volume, no significant difference between the undeformed  Frontiers in Bioengineering and Biotechnology | www.frontiersin.org fiber volume fraction (28.6 ± 4.95%) and that at peak strain (30.1 ± 7.35%; p = 0.524) (Figure 5B) was found.
In order to determine whether fiber reorganization was affine with stretch or not, the experimentally measured fiber kinematics in the mechanically deformed images were compared to the kinematics resulting from an affine transformation of the same images (Figures 6A-C). The measured 3D directional variance within the tetrahedron in mechanically deformed image volumes (0.457 ± 0.101) was significantly higher (p = 0.012) than the variance under a theoretical affine transformation (0.372 ± 0.100) (Figure 6D). These differences were conserved along multiple anatomical planes. A higher average directional variance was measured for orientations in the x-y plane (θ) during deformation (0.380 ± 0.091; p = 0.015) compared to the theoretical affine transformation (0.314 ± 0.099) (Figure 6E). A higher average directional variance was also measured (0.467 ± 0.141) for fiber inclination out of the x-y plane (ϕ) compared to the theoretical affine transformation (0.322 ± 0.086), but these changes did not meet the criteria for statistical significance (p = 0.056) ( Figure 6F). When these results are viewed as a change in directional variance relative to the undeformed configuration (Figures 6G-I), it becomes apparent that the actual fiber realignment in the direction of loading is less than what was calculated for an affine transformation.
Dermal collagen is often described as having a basket-weave like organization, which has previously been found to vary in histological sections as a function of the length scale over which it is measured (Fisher et al., 2008). To quantify this variability in 3D SHG image volumes, the average and standard deviations of regional fiber volume fraction and directional variance were computed over different scales within each undeformed SHG image volume. When the average regional directional variance was compared across length scales, it was found to significantly increase with increasing length scale (p ≤ 0.01) until it plateaued around a scale of 40 µm 3 (Figures 3C, 7A). Interestingly, the standard deviation of regional directional variance was not significantly different between volumes of 4 3 and 10 3 µm 3 (p = 0.99), but it did decrease significantly as the length scale increased from this range (p ≤ 0.001) (Figure 7C). The average regional fiber volume fraction remained constant from 4 3 to 40 3 µm 3 (p ≥ 0.09), followed by a slight increase at larger length scales (p < 0.001) (Figures 3D, 7B). As length scale increased, the standard deviation of regional fiber volume fraction decreased (p ≤ 0.01) ( Figure 7D). FIGURE 6 | Collagen fiber kinematics are not affine with the local skin deformation. Micro-scale image volumes were compared among the undeformed configuration (A), actual tissue response at peak strain (B) and the theoretical affine transformation (C) based on the measured deformation gradient. The white tetrahedron in panels (A-C) indicates the position of 3D image features that were tracked during deformation (scale bar represents 100 µm). The computed 3D directional variance (D) as well as the measured directional variance for individual 3D orientation and inclination components [θ and ϕ, panels (E,F), respectively] were found to be higher than the theoretical affine transformation. Similarly, the change in overall directional variance (G) and the 3D orientation and inclination components (H,I) indicate the measured fiber re-alignment in the direction of loading is less than the theoretical affine response. Points are colored to indicate paired comparisons among groups. *represents p < 0.05. FIGURE 7 | In the undeformed configuration, regional averages of directional variance (A) and fiber volume fraction (B) increased as a function of the length scale over which they were measured. Alternatively, the average regional standard deviations of directional variance (C) and fiber volume fraction (D) decreased at larger length-scales, indicating the importance of length-scale when implementing measurements into computational models. Points are colored to indicate paired comparisons among groups. *represents p < 0.05.

DISCUSSION
In this study, 2D macro-scale and 3D micro-scale tissue kinematics were quantified simultaneously during mechanical stretch. Historically in biomechanics, it has often been assumed that soft tissues are nearly incompressible (Lees et al., 1991). In this study, however, deformation at the macro-scale revealed a very large Poisson's ratio for mouse skin. Recently others have reported large Poisson's ratios both computationally (Reese et al., 2010;Ban et al., 2019) and experimentally in tendon (Lynch et al., 2003;Lake et al., 2010;Reese and Weiss, 2013), cartilage (Elliott et al., 2002), and skin (Pissarenko et al., 2019;Wahlsten et al., 2019). These observations may be related to a contractile force orthogonal to the direction of stretch caused by collagen fiber realignment (Figure 5A; Ban et al., 2019). The large Poisson's ratio phenomenon has been identified to be mechanically anisotropic (Elliott et al., 2002) and location dependent (Pissarenko et al., 2019). Although Poisson's ratio can be estimated during uniaxial stretching, the material is unconstrained in the axis orthogonal to the direction of stretch, and biaxial mechanical stretching of skin is required to fully characterize the force-coupling that exists between multiple axes during stretch. Interestingly, an examination of the tissue's micro-scale kinematics revealed a substantial decrease in volume ( Figure 4C) and significant compression along both directions orthogonal to the direction of stretch. Although previous studies have described a decrease in the micro-scale volume of collagen gels during uniaxial stretch (Roeder et al., 2004), and macro-scale volume decreases in murine and human skin (Wahlsten et al., 2019), a decrease in the micro-scale volume of skin has not been explicitly measured before. Volume loss during deformation is thought to be due to interstitial fluid flowing out of the dermal ECM (Roth and Mow, 1980;Wahlsten et al., 2019), but more research is needed to verify this hypothesis and determine if other mechanisms are involved.
A significant decrease in micro-scale volume was measured during uniaxial stretch in this study, but the fraction of voxels that contain fibers did not change significantly during stretch ( Figure 5B). This lack of a change in volume fraction indicates that the collagenous regions of the dermis are decreasing in volume in a similar amount as the non-collagenous regions of the skin during stretch. Given that the mechanical properties of various skin components are different, this result is surprising. However, it is important to consider the resolution of these volume fraction measurements. MPM imaging of collagen fibers via SHG is a diffraction-limited microscopy technique that is also limited by the digital resolution of individual voxels during image acquisition (Campagnola and Dong, 2011). It remains quite possible that the actual density of collagen fibers within the SHG-positive voxels may be changing even if the overall volume fraction of SHG-positive voxels does not change. Previous experiments have determined SHG intensity is sensitive to changes in sub-voxel fiber density due to the relationship between measured intensity and fiber diameter (Bancelin et al., 2014). However, we cannot definitively assess fiber density within the SHG-positive voxels based on their SHG intensity in intact skin and other thick samples for a number of reasons. Photon scattering complicates the relationship between SHG intensity and fiber diameter/density (Campagnola and Dong, 2011). Other experiments have also determined that the organization of fibrils within the fiber directly affects the measured SHG intensity due to various phase-matching conditions (Campagnola et al., 2002). Therefore, additional future work is needed to provide a complete assessment of how the density of collagen fibers may change at various scales in intact skin during mechanical stretching.
Previous work that focused on micro-scale deformation of mouse skin under uniaxial stretch was based on the relative position of hair follicles in 2D within SHG image volumes. Additionally, the skin used in those studies (Bancelin et al., 2015;Lynch et al., 2017a,b) was excised from the back of the mouse, stripped of both the hair and epidermis, and was mechanically stretched as a strip of tissue rather than as the traditional dogbone shape used in our current study. Although micro-scale volume changes were not computed in previous studies by Lynch et al., they measured a much smaller contraction in the transverse direction compared to our current study (Figure 2D), which could be due to differences in excisional location, tissue orientation relative to loading direction, and tissue geometry. Additionally, by utilizing features within the collagen network rather than using hair follicles as reference points for computing the micro-scale deformation, our strain measurements are more reflective of the kinematics and deformation of the collagen microstructure. Multi-scale 3D strain measurements during uniaxial tensile loading of soft tissue are relatively rare, but previous work involving multi-scale strain measurements of collagen gels showed that the amount of micro-scale strain about the axis of stretch was less than at the macro-scale (Roeder et al., 2004). In contrast, no significant differences were found between macro-and micro-scale strain measurements in this study ( Figure 4B). However, skin is more complex and denser than typical collagen gels, and multi-scale transmission of strain may be aided by the more complex and dense collagen organization in the dermis or interactions with other layers of the skin (Blair et al., 2020).
To provide context for the changes in tissue strain at the micro-scale, the 3D fiber alignment was quantified during mechanical stretching. During stretch, 3D directional variance measurements in this study showed a significant increase in fiber alignment (Figure 5A), similar to previous findings from 2D fiber organization measurements within skin (Lynch et al., 2017a;Ducourthial et al., 2019). Most computational models of soft tissue kinematics assume that the micro-scale fiber kinematics are primarily driven by tissue-fiber interactions, resulting in an affine relationship between tissue deformation and fiber kinematics (Billiar and Sacks, 1997;Sacks, 2003;Lake et al., 2012;Fan and Sacks, 2014;Lee et al., 2015;Quinn et al., 2016). In skin however, previous findings based on 2D fiber orientation measurements indicated that this relationship is nonaffine, and that the affine assumption results in an overestimate of fiber re-organization (Hepworth et al., 2001;Jayyosi et al., 2017). Results from this study ( Figure 6D) agree with those findings, showing that the 3D fiber kinematics shares a nonaffine relationship with tissue deformation and that the affine assumption overestimates the fiber alignment in the direction of stretch. Our results uniquely show this non-affine fiber behavior in 3D and demonstrates it occurs both in and out of the x-y plane of the skin (Figures 6E,F). The non-affine relationship between tissue deformation and fiber kinematics may be due to fiberfiber interactions or interactions between the multiple layers of skin or its appendages. To elucidate the key factors driving this non-affine fiber response and understand other complex tissue behaviors, computational modeling will be critical in the future. Modeling efforts have greatly increased in complexity from neo-Hookean strain energy functions proposed decades ago (Fung, 1967) to more recent Fung-type hyperelastic models (Meador et al., 2020) and modern multi-scale simulations involving fiber networks (Sander et al., 2009;Fang and Lake, 2016;Ban et al., 2019). Future work will focus on integrating detailed 3D experimental results with multi-scale image-based computational models of skin mechanical behavior.
In summary, our 3D micro-scale measurements of tissue deformation under uniaxial tension revealed a significant decrease in tissue volume that is contradictory to traditional neo-Hookean or Mooney-Rivlin models previously used in skin applications (Fung, 1967;Weiss et al., 1996). Our 3D measurements of the orientation of individual fibers in this study suggest that fiber kinematics are not affine to the local deformation ( Figure 6D). This indicates that models allowing for non-affine kinematics, such as those involving fiber network simulations may be necessary to accurately model the mechanical behavior of skin. We found that summary statistics of collagen fiber alignment and volume fraction are also very dependent on the length scale over which the data are summarized (Figure 7). Due to this length scale dependence, experimental results used to build and test image-based fiber network models need to have appropriate matching length scales (Sander et al., 2009). When interpreting these experimental findings and integrating them with computational modeling, several limitations should be considered. This study focused on the kinematics of collagen fibers, but characterization of other ECM constituents may be important in understanding mechanical phenomena such as the change in tissue volume. To facilitate 3D imaging in this study, an incremental, quasi-static loading scheme used, which does not allow for exploration of the strain-rate dependence of the kinematic response. In addition, skin is multi-axially loaded in vivo and attached to the underlying fascia, conditions that undoubtedly put additional constraints on fiber kinematics not present in the uniaxial in vitro testing scheme. Finally, it is also important to note that these findings apply to mouse skin, which differs in many notable ways from human skin (Wei et al., 2017). Nonetheless, this study provides unique quantitative insights into skin kinematics at different scales that can not only aid in establishing more biofidelic computational models, but also provide a strong foundation to assess how changes in tissue structure can affect the mechanical properties of the tissue and clinical outcomes.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.