A model of human skin under large amplitude oscillatory shear

Skin mechanics is of importance in various fields of research when accurate predictions of the mechanical response of skin is essential. This study aims to develop a new constitutive model for human skin that is capable of describing the heterogeneous, nonlinear viscoelastic mechanical response of human skin under shear deformation. This complex mechanical response was determined by performing large amplitude oscillatory shear (LAOS) experiments on ex vivo human skin samples. It was combined with digital image correlation (DIC) on the cross-sectional area to assess heterogeneity. The skin is modeled as a one-dimensional layered structure, with every sublayer behaving as a nonlinear viscoelastic material. Heterogeneity is implemented by varying the stiffness with skin depth. Using an iterative parameter estimation method all model parameters were optimized simultaneously. The model accurately captures strain stiffening, shear thinning, softening effect and nonlinear viscous dissipation, as experimentally observed in the mechanical response to LAOS. The heterogeneous properties described by the model were in good agreement with the experimental DIC results. The presented mathematical description forms the basis for a future constitutive model definition that, by implementation in a finite element method, has the capability of describing the full 3D mechanical behavior of human skin.


Introduction
Skin is the largest organ of the human body and functions as a barrier against the external environment. It protects against excessive water loss from the aqueous interior, the ingress of pathogens and ultraviolet light. Human skin is composed of several layers, each with a unique composition and function. The outermost layer, the epidermis, consists of the stratum corneum and the viable epidermis. The stratum corneum is composed of terminally differentiated non-viable keratinocytes, called corneocytes. It has a thickness varying from 10 to 20 μm depending on body site and is the most important protection barrier. The viable epidermis consists of differentiating keratinocytes. Originating from the basal layer, keratinocytes differentiate while migrating upwards to the stratum corneum. By doing so, over time they lose their nuclei, become flattened, discharge lipids and are cornified. Eventually they form the stratum corneum. Via the dermo-epidermal junction, a basal membrane, the epidermis is connected to the underlying dermis. This connective tissue is mainly composed of collagen and elastin fibers embedded in an extrafibrillar matrix. The dermis can be divided in two layers. The papillary dermis is the uppermost layer, composed of loose areolar connective tissue forming fingerlike projections called papillae, which extend towards the epidermis. It contains terminal networks of blood capillaries and nerve endings. The underlying reticular dermis is a thicker layer, composed of densely packed collagen and elastin fibers. The presence of these fibers are important factors in the mechanical properties of human skin, providing both strength and elasticity. The fiber alignment is responsible for the presence of anisotropy, well known as the Langer's lines (Langer, 1978). Other components of the dermis are hair follicles, sebaceous glands and sweat glands. Underneath the dermis lies the subcutaneous tissue layer, mainly composed of adipose cells used for fat storage, and is generally not regarded as part of the skin.
The mechanical integrity of the skin can be threatened by diseases, trauma and medical or cosmetic treatments. Therefore, skin mechanics is important for various fields of research. This includes research into the development of pressure ulcers and the interaction between skin and devices or materials such as shaving appliances, prosthetic liners, bed linens and medical devices. A constitutive model capable of describing the mechanical behavior of human skin with meticulous experimental validation is thus needed. Such a model can be used to optimize the design of materials and devices involved in skin interaction, but also to improve the understanding of the mechanical behavior of human skin. The intricate structure of human skin, as described before, results in complex mechanical behavior. The absence of satisfactory constitutive skin models is partly due to difficulties in measuring this behavior. The mechanical response of human skin is highly dependent on the loading modality and can vary many orders of magnitude between shear, tensile and indentation behavior (Leyva-Mendivil et al., 2015). Since it is impossible to capture the full mechanical response of human skin in one single experiment, different types of experiments have been performed to partially capture this response. Most studies involved experiments on in vivo human skin, comprising uniaxial tensile (Mahmud et al., 2012;Delalleau et al., 2008b;Gambarotta et al., 2005), multiaxial tensile (Flynn et al., 2011(Flynn et al., , 2011(Flynn et al., , 2013Khatyr et al., 2004;Kvistedal and Nielsen, 2009), suction (Khatyr et al., 2006;Delalleau et al., 2008aDelalleau et al., , 2011Delalleau et al., , 2012Hendriks et al., 2003) and indentation (Groves et al., 2012;Delalleau et al., 2006;Kumar et al., 2015) experiments. Although it is preferable to measure skin in its natural environment, in vivo measurements have limitations for several reasons. Firstly, it is nearly impossible to separate the response of the skin from that of the underlying subcutaneous tissue. Secondly, because of the nature of the tests (the area of interest is not isolated from the environment) it is almost always necessary to perform a quite complicated mechanical analysis to account for the environment and the often complex boundary conditions in such a test. This means that the mechanical analysis has to be done by numerical methods and an inverse approach is needed to determine material properties. The success of such an approach strongly depends on the quality of the constitutive model, but it is very difficult to formulate and improve constitutive models based on these rather complicated tests.
Tests on ex vivo skin can be done in a much more controlled way with very clear boundary conditions, requiring less complicated analysis. However, it is difficult to obtain, prepare and preserve the tissue samples as the properties may change after the tissue is collected. The ex vivo tests can very well be used to develop the mathematical formulation of constitutive models and to determine initial estimates for material parameters. Then, using these constitutive models, in vivo tests could be used to fine-tune and individualize the material parameters.
Measurements on ex vivo animal skin were performed abundantly, mostly in uniaxial (Groves et al., 2013;Li and Xiaoyu, 2016;Lokshin and Lanir, 2009a;Bischoff, 2006;Karimi et al., 2015Karimi et al., , 2016Limbert, 2011) tests and a single multiaxial (Jor et al., 2011) tensile test. However, because of the well-known anatomical, biological and mechanical differences, it is difficult to translate these results to human skin. A justifiable compromise is thus to perform experiments on ex vivo human skin, since this allows for measuring the tissue of interest under highly controllable conditions. Various experiments on ex vivo human skin were used for evaluation and validation of existing constitutive skin models, including uniaxial tensile (Groves et al., 2013;Li and Xiaoyu, 2016;Lapeer et al., 2010), indentation (Geerligs et al., 2011a) and shear (Geerligs et al., 2011b) tests.
In the current paper we focus on the skin behavior under shear loading. Holt et al. (2008) modeled the shear response of ex vivo human skin to physiologically relevant frequencies using a 1D viscoelastic modified Kelvin-Voigt model. However, due to experimental limitations, they only were able to describe the response for very small strains up to 0.005. Geerligs et al. (2011b) modeled the shear response of the stratum corneum and viable epidermis of ex vivo human skin samples up to strains of 0.01 measured with a controlled rheological setup. Using a linear mixed model definition, the dynamic shear moduli of these layers were described taking into account the significant influence of both temperature and relative humidity. However, studies presenting a constitutive description of the shear response of human skin are all limited to the linear regime.
Experimentally, Gerhardt et al. (2012) and later (Lamers et al., 2013) studied the response of ex vivo human skin using large amplitude oscillatory shear (LAOS) showing highly nonlinear viscoelastic behavior for shear strains up to 0.1. Additionally, complex mechanical phenomena such as intra-cycle strain stiffening and inter-cycle shear thinning were observed. Using digital image correlation (DIC), heterogeneity was observed as gradually decreasing shear moduli for increasing skin depth. However, to the authors knowledge, no constitutive model presented in literature is capable of describing this complex mechanical response of full-thickness human skin to shear deformation.
With this study we present an improved experimental method, based on the previous work of Gerhardt et al. (2012) and Lamers et al. (2013), capable of measuring the nonlinear viscoelastic mechanical response of full-thickness human skin to LAOS for strains up to 0.2. This was combined with DIC to assess tissue heterogeneity. More importantly, the aim is to capture this mechanical response to LAOS with a proposed 1D mathematical description, including all observed complex phenomena, i.e. intra-cycle strain stiffening, inter-cycle shear thinning, strain softening and nonlinear viscous dissipation. This mathematical description forms the basis for a future constitutive model definition that, by implementation in a finite element method, has the capability of describing the full 3D mechanical behavior of human skin.

Experimental setup
Human skin was obtained from abdominoplastic surgery at the Catharina Hospital, Eindhoven, The Netherlands, according to Dutch guidelines of secondary used materials. Patients (aged 18-65 years) gave informed consent for the use of their skin for research purposes. The skin was processed immediately after surgery. First it was stretched on a stainless steel plate using forceps (Geerligs et al., 2011b). Fullthickness skin slices, thickness 1.5-2 mm, were obtained using an electric dermatome (D42, Humeca, The Netherlands) and stored at − 30°C until further use. After thawing at room temperature, 8 × 8 mm 2 skin samples were punched out and stored in phosphate buffered saline (PBS). Large amplitude oscillatory shear measurements on a single skin sample were performed as described by Lamers et al. (2013). In short, a skin sample was adhered eccentrically between the bottom plate and a 50 × 8 × 5 mm 3 custom made steel bar, radius (R) 25 mm, on a strain controlled rheometer (ARES LS-LC, Rheometric Scientific, USA) using glue (Bison Tipper, Bison International B.V., Goes, The Netherlands). Using a Thermo Neslab RTE-10 Digital Plus refrigerated bath (Artisan Technology Group, IL, USA) and internal Peltier element the temperature was maintained at 34°C. The skin sample was placed in the mechanical setup by first gluing the epidermal side onto the lateral end of the top bar. By placing the sample eccentrically the sensitivity of the measurement was increased. After attaching the bar to the rheometer, it was lowered onto the glue containing bottom plate with a normal force between 10 and 30 mN. LAOS measurements were performed by applying oscillatory shear deformation with strain amplitudes (γ 0 ) of 0.01, 0.05, 0.1, 0.15 and 0.2 at a frequency (f) of 1 Hz. Each strain amplitude was measured for 40 s, of which only the last two cycles were used for the analysis. After initial stress relaxation the response reached steady state at that time (see Fig. 3). Raw strain and torque signals were obtained by connecting the analog outputs of the rheometer to a stand-alone pc via an analog-to-digital converter. A sinusoidal strain γ t ( ) was applied to the bottom plate, described as: with γ 0 the strain amplitude and ω the angular frequency (2 πf ). The actual shear stress τ was calculated as a function of the measured torque M via: with M the measured torque, l the width of the sample (8 mm) and R the radius of the top bar (25 mm). This nonlinear shear response upon sinusoidal deformation was then filtered by describing it as a Fourier series (Ewoldt et al., 2008;Hyun et al., 2011;Wilhelm, 2002): where ′ G n and ″ G n represent the Fourier coefficients that are related to the viscoelastic moduli. In the most simple case of a linear viscoelastic regime, only ′ G 1 and ″ G 1 remain which are more commonly known as the storage and loss modulus, respectively.
Elastic and viscous Lissajous curves, where the stress is plotted versus the imposed shear strain or strain rate, respectively, were used to describe the nonlinear viscoelastic response of human skin to LAOS. Linear viscoelastic behavior appeared as an elliptic Lissajous curve, whereas deviation from this elliptic shape indicated nonlinear behavior. The focus of this study was to find a mathematical description for the observed complex mechanical response of human skin in the LAOS experiment that can be used as a basis for a constitutive model describing the full 3D mechanical behavior. Taking into account the rather similar response of different skin samples, the results of one representative skin sample were used in this study Fig. 1.

Digital image correlation
One cross-sectional area of the skin sample in the rheometer was visualized using a monochromatic CCD-camera (Imaging Source) mounted on a stereo microscope (SZ11, Olympus, Germany). Heterogeneity as previously observed by Lamers et al. (2013) was assessed on this cross-sectional area. To enhance contrast a speckle pattern was manually applied to this area using graphite spray (Graphit 33, CRC Industries Kontakt Chemie, Iffezheim, Germany) and a cotton bud. At a frame rate of 60 fps, 8 bit gray-value (640 × 480 pixels) movies were acquired. With a resolution of 3.5-4.3 μm/pixel a field of view (FOV) between 2.3 and 2.8 mm in width and 1.7-2.1 mm in height was produced. Gerhardt et al. (2012) proved that the dimensions of the skin sample captured with the current FOV are sufficient to represent the deformation of the total skin sample. These movies were converted into an image sequence using FIJI software (ImageJ, La Jolla, USA). Determination of local displacement was performed via DIC in Matlab 2014b (Mathworks, Natick, MA, USA) using a strain tracking code kindly provided by Victor Barocas' group from the University of Minnesota (Raghupathy and Barocas, 2010). A region of interest was defined by specifying a rectangular grid of 21 × 11 tracking points in xand y-direction, respectively (see Fig. 2). This results in 10 layers of 20 elements each. Local strains were calculated by tracking the displacement of each point over time. The shear strain in each of the 10 layers was calculated by averaging over all 20 elements. Heterogeneity was defined as the shear strain distribution over all 10 layers in y-direction at maximum displacement amplitude. The skin behavior was considered to be isotropic during LAOS and the resulting oscillatory shear strain in each layer was described with a sine function (see Fig. 3). In order to quantitatively compare experiment and model, the 10 layers were subdivided in epidermis (1-2), papillary region (3-4) and reticular dermis (5-10). Low resolution in the graphite speckle pattern in the top layer (see Fig. 2) led to instability in the determination of local displacements with DIC. Therefore only values of the second layer were used for the epidermis.

Constitutive model build-up
The skin is modeled as a one-dimensional layered structure, with = … i N 1, 2, , layers layers (see Fig. 4). Every sublayer behaves as a nonlinear Standard viscoelastic model (Fig. 5). The stiffness of the parallel spring is defined with k 2 , the spring in the viscoelastic element with k 1 and the viscosity of the dashpot with .
The shear strain γ i and shear strain rate γ˙i in each sublayer were calculated as follows: with u i and h i the displacement and thickness of each sublayer, respectively. Nonlinearity in the LAOS response was observed in the experimental results (see Fig. 6) by Lamers et al. (2013) as intra-cycle strain stiffening, indicated as a transition from concave to convex form in the elastic Lissajous curve. Following van Kempen et al. (2015) this is included in the mathematical model by assigning nonlinear parameters to both elastic components as follows:    with heterogeneity implemented by assigning an initial stiffness distribution to k i 02, over all sublayers i, as experimentally determined with DIC. The parameter x i describes a strain history dependent change in material properties and is explained in detail later. The observed nonlinear strain dependent viscous dissipation was described by assigning similar nonlinearity to the viscous element: The model response per sublayer can be described with the following equation: with τ the shear stress (Pa). By implementing the following boundary conditions these equations were solved using a recursive scheme: As shown by Lamers et al. (2013), the stiffness of the skin decreases during multiple deformation cycles of the LAOS experiment. The increasing amplitude of the cyclic deformation during LAOS leads to a lower stiffness. Subsequent decrease in strain amplitude results in an increase in the stiffness, making the effect fully reversible. This softening effect is best illustrated by the minimum-strain modulus ′ G M (Ewoldt et al., 2008), defined as the local derivative of the shear stress with respect to strain at γ = 0: It decreases with increasing γ, as shown in Fig. 6, indicating the softening effect. Analogous to Münster et al. (2013) and van Kempen et al. (2015), this effect was quantified with the introduction of a network state parameter (NSP), x, which describes the evolution of the mechanical properties based on the deformation history: i i were experimentally determined for the total skin response at each γ 0 . This definition was generalized with respect to γ i by describing the NSP as an exponential function through x τ γ ( , ) i i : with fitting parameters a and b. The softening effect is a material property which is dependent on the strain history γ max in each individual layer i, in this case of LAOS defined as: with t ω is currently set as the period of a single oscillation (1 s). Recovery of the NSP x is included with the dependency on γ max i , . The constitutive skin model definition was implemented in Matlab 2014b.

Iterative parameter estimation
An initial estimation of the six model parameters, k 01 , k 02 , η 0 , m 1 , m 2 and m 3 was performed by manually fitting the simulated and experimental data. Additionally, an adaptive sparse generalized polynomial chaos expansion sensitivity analysis (Quicken et al., 2016) was performed to investigate the influence of parameter value variation on the model output. Subsequently, the error ̰ ξ between the experimental and simulated data was minimized using an iterative parameter estimation where ̰ m is a vector containing experimental data. The vector ̰ h contains the two sets of simulation output corresponding to boundary conditions ̰ u and parameter values ̰ θ. First, the total shear strain in each increment, resulting from the sinusoidal displacement as described in Eq. (8)  Both were determined at strain amplitudes γ 0 = 0.01-0.2 in all increments (N inc ).
To minimize the influence of startup effects in the simulation, five cycles were simulated and only the last was used for the iterative parameter estimation. The error was minimized by means of the Jacobian J, which was defined as follows: where the weighting matrix V assigns a weight to the errors, and was defined as: where δ d and δ F are the uncertainty in displacement and force, respectively, caused by the limited accuracy of the measurement systems. They were set to = ⋅ − δ 3.5 10 mm where ̰ e j is a P-dimensional column of which the j th entry is one and all other entries are zero. A choice for Δθ j , a small variation of parameter P j , was made as: Δ rel is set to − 10 3 and Δ abs to a value of − θ 10 j 3 (0) , with θ j i ( ) the estimate of the parameter for each iteration i. The step size per iteration ̰ δ θ i ( ) was calculated with a so-called Gauss-Newton linearization equation: which has quadratic convergence in the neighbourhood of the optimal solution. This resulted in an iterative scheme to update the parameters ̰ θ i ( ) : This iterative procedure was continued until the change in the parameter updates was smaller than a critical value δ θ , set at − 10 4 : The normalized step size ̰ δ θ i ( ) was calculated as follows: Here the factors δ θ and δ j prevented numerical problems if θ j i ( ) approached zero. During the entire iteration process, no model parameters were allowed to become negative.

LAOS: overall results
The overall skin response altered from linear at = γ 0.01 0 to increasingly nonlinear stress response for strain amplitudes up to = γ 0.2 0 , as shown in Fig. 6. In accordance with Lamers et al. (2013), two phenomena were observed for increasing strain amplitudes: intracycle strain stiffening, shown as a concave to convex deviation for high strains, and inter-cycle shear thinning, shown as a deviation from elliptical to non-elliptical in the viscous Lissajous curve. In accordance with van Kempen et al. (2015) two more phenomena were observed: a Fig. 7. Shear strain profile through full-thickness skin averaged over 5 oscillations at γ = 0.1, with the average (solid) and standard deviations (dashed). The red dash-dotted vertical line denotes the computed average strain across the entire skin thickness. softening effect, indicated by a decreasing minimum-strain modulus ′ G M (see Eq. (9)) with increasing γ 0 , and (nonlinear) strain dependent viscous dissipation, shown by the intra-cycle broadening of the elastic Lissajous curves with increasing strain amplitude.

Digital image correlation
DIC was performed on images with strain amplitudes ranging from = − γ 0.01 0.2 0 at a frequency of 1 Hz. The DIC software was able to capture the local displacements of the skin at the defined grid points, see Fig. 2. Heterogeneity over the depth of the skin was assessed and it was shown that the skin stiffness gradually decreases with depth (see Fig. 7), i.e. the dermis is softer than the epidermis. These results are similar to the findings of Gerhardt et al. (2012). The resulting, manually fine-tuned, values for the heterogeneity of k i 2,0 as mentioned in Eq. (5) are: with subdivision in epidermis, papillary region and dermis indicated by the dots. As larger strains were present in the dermis it is dominating the nonlinear shear response of the skin. Indicated by the relative small area enclosed by the elastic Lissajous curves, the elastic component of skin dominates the viscous component.

Constitutive model
The iterative parameter estimation method had a high rate of convergence and resulted in a unique solution for a wide range of initial parameter values. Since all material parameters were estimated simultaneously it is most likely that the found solution is a global minimum. The resulting set of parameter values are shown in Table 1. The numerical results for the LAOS sequence are shown in Fig. 8 and compared with both the elastic and viscous Lissajous curves from the experiments. All complex mechanical phenomena observed in the experimental results, i.e. intra-cycle strain stiffening, inter-cycle shear thinning, the softening effect and nonlinear viscous dissipation, are captured with the constitutive model. This includes the maximum stress values for all γ 0 .
The NSP x describes the change in material properties as a function of the strain history γ max in each layer i, as defined in Eq. (11). The evolution of γ max i , and x i during the LAOS sequence, with applied shear strain γ 0 = 0.01-0.2, is shown in Fig. 9. Due to the strain dependent viscoelasticity in the constitutive model, the distribution of γ i evolves upon each step in γ 0 . This results in an initial underestimation of γ max i , in the epidermis and papillary region and a corresponding overestimation in the dermis, gradually reaching equilibrium during following oscillations. The values of the NSP x i , as defined in Eq. (11), show corresponding but inversely proportional behavior. The resulting softening effect was captured accurately, indicated by the decreasing ′ G M with increasing γ 0 in a similar manner to the experimental results. However, there is a noted deviation in slope at ≥ γ 0.15 0 . Finally, the nonlinear viscous dissipation is in good agreement with the experimental results which is also indicated by the excellent fit of the viscous Lissajous curve which includes the observed shear thinning effect.
With the decomposition of the skin in epidermis, papillary region and dermis, the results of the constitutive model and experiment were quantitatively compared (see Fig. 10). Overall, the model results for all three sublayers are in good agreement with the experiment. For all γ 0 the stiffness decreased with increasing depth and lowest strains were thus observed in the epidermis, with mainly linear elastic behavior for all γ 0 . At γ 0 = 0.01 all three sublayers are in the linear regime, for higher γ 0 the behavior becomes increasingly nonlinear viscoelastic for the deeper layers.

Discussion
This study aimed to develop a new constitutive material model for human skin that is capable of describing the nonlinear viscoelastic mechanical response of ex vivo human skin under LAOS deformation. To the best of our knowledge this is the first time this type of behavior has been captured in a constitutive model of human skin. Holt et al. (2008) were able to model the mechanical response of ex vivo human skin to low magnitudes of shear using a Kelvin-Voigt model. With this approach they were able to describe the shear creep response, the modelling of oscillatory data had more limited success. Especially the determination of ″ G , representing the viscous portion of the response to oscillatory shear, was inaccurate. Geerligs et al. (2011b) described the linear shear response of the stratum corneum and viable epidermis with a mixed model definition. With our model both the elastic and viscous Table 1 Parameter values as determined by the iterative parameter estimation method. behavior can be described as shown with the Lissajous curves in Fig. 8. In addition, we were able to capture the nonlinear response resulting from large shear deformations up to γ 0 = 0.2. Compared to Lamers et al. (2013) improvements in image acquisition and especially the determination of local skin displacements have been made by increasing the frame rate of the camera and the use of novel DIC software (Raghupathy and Barocas, 2010). This enabled stable local displacement analysis at large strain amplitudes and associated strain rates. Although strains above γ 0 = 0.2 were not used in this study due to potential skin damaging effects, the experimental method proved to be stable up to γ 0 = 0.6 at a frequency of 1 Hz. The use of ex vivo human skin under highly controllable laboratory conditions has several benefits. First, the skin was isolated and no unknown influences of surrounding tissue were present, which is always the case with in vivo experiments. This allowed for focussing exclusively on the mechanical response of skin and made the experiments highly reproducible. Secondly, the boundary conditions could be defined accurately. Obviously, the downside of this method is that the Fig. 9. Applied γ 0 during the LAOS sequence, with the corresponding strain history γ max i , in the epidermis (solid), papillary region (dashed) and dermal layer (dotted) and therefrom calculated evolution of the NSP x in every layer. results are not entirely representative for healthy in vivo skin. Although the preparation, preservation and handling of skin samples has been performed under optimized conditions, analogous to Geerligs et al. (2011b), the mechanical response may be affected. Possibly even more important is the origin of the skin samples, namely from patients that underwent abdominoplastic surgery. Although readily available, the strain history of the expanded skin may include changes in its structural properties when compared to healthy skin. Moreover, it is not representative for all skin as the structural and mechanical properties vary with anatomic location.
Final assessment of heterogeneity was done by manual fitting through the calculated average shear strain of all 10 layers of the DIC grid at maximum displacement amplitude. By definition this leads to a stepwise stiffness decrease over the skin layers. However, in line with the findings of Lamers et al. (2013) this can also be approximated with a gradually decreasing function. Despite the fact these layers are biologically distinct, this gradual decrease can also be rationalized from a biological perspective since the composition of the epidermis and dermis, both cellular and fibrous, changes with depth (Venus et al., 2011). First, the epidermis consists of several distinct strata which eventually merge in the dermal-epidermal junction. Moreover, the connectivity of the dermis decreases with increasing depth, lowering the stiffness. It is therefore justifiable to describe the heterogeneity with a gradually decreasing exponential function, in line with previous findings from Lamers et al. (2013).
Since the fibers are assumed to have an in-plane orientation, their (anisotropic) contribution in response to pure shear loading was neglected. This is in agreement with Ní Annaidh et al. (2012), who state that the preferred orientation of the three-dimensional fiber network lies parallel to the surface. The three-dimensional nature of these fiber networks, however, implies the presence of out-of-plane fiber orientations which do contribute to the mechanical response to shear loading. Their separate contribution is neglected, but is part of the described total shear response of the skin.
The model is able to represent the viscous hysteresis under cyclic loading of the LAOS experiment, acting on a short characteristic time scale. However, viscosity on a larger characteristic time scale is also observed in the mechanical response of human skin in the form of stress relaxation. Karimi et al. (2016) and Lokshin and Lanir (2009b) used one-and two-term Prony series, respectively, to describe this viscous behavior under uniaxial deformation. In the current model it can be implemented kinetically in the NSP x, in analogy with van Kempen et al. (2015). Subsequent validation could be performed with relaxation experiments at constant step shear strain.
The LAOS experiments have been performed at a constant frequency of 1 Hz, implying the strain rate to increase with the applied strain. As mentioned by Lokshin and Lanir (2009b), the skin response is strain rate dependent and thus changes over the range of γ 0 = 0.01-0.2. The current constitutive model proved to be capable of describing the mechanical response at this specific frequency. However, it is unknown if the model is capable of describing the skin response at other frequencies.
In the current study we developed a constitutive model of human skin that is capable of describing the heterogeneous, nonlinear viscoelastic mechanical response to LAOS. However, this does not mean the current model is capable of describing the mechanical response under different types of loading. It for instance lacks fiber constituents providing anisotropic properties. The future goal is therefore to develop a constitutive model that is capable of describing the full 3D mechanical behavior of human skin. This can be achieved by the addition of an anisotropic fiber contribution to the current model, preferably implemented in a finite element method. Unambiguously this requires multiaxial experimental validation, i.e. biaxial tensile testing, on ex vivo human skin. For the field of skin research, having a constitutive model that is able to describe the full 3D mechanical behavior of human skin would be invaluable. The currently developed constitutive model forms a solid basis for such a 3D model.

Conclusions
The new constitutive model presented in this study proved to be capable of accurately describing the complex mechanical response of human skin to LAOS. With the experimental setup combining LAOS and DIC on ex vivo human skin the heterogeneous, nonlinear viscoelastic response was determined. Using an iterative parameter estimation method the model parameters were optimized. The model captured the experimentally observed strain stiffening, softening and increasing viscous dissipation, which were all present in the mechanical response to LAOS. The heterogeneous properties described by the model were in very good agreement with the experimental DIC results.