Home | About | Journals | Submit | Contact Us | Français |

**|**Proc Math Phys Eng Sci**|**PMC5549575

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Basic structural anatomy and biophysical characteristics of the skin
- 3. Modelling of the skin
- 4. Constitutive models of the skin
- 5. Application: skin wrinkles
- 6. Perspective and conclusion
- Supplementary Material
- References

Authors

Related links

Proc Math Phys Eng Sci. 2017 July; 473(2203): 20170257.

Published online 2017 July 26. doi: 10.1098/rspa.2017.0257

PMCID: PMC5549575

email: ku.ca.notos@trebmil.g

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3823390.

Received 2017 April 9; Accepted 2017 June 21.

Copyright © 2017 The Author(s)

Published by the Royal Society. All rights reserved.

The objective of this paper is to provide a review on some aspects of the mathematical and computational modelling of skin biophysics, with special focus on constitutive theories based on nonlinear continuum mechanics from elasticity, through anelasticity, including growth, to thermoelasticity. Microstructural and phenomenological approaches combining imaging techniques are also discussed. Finally, recent research applications on skin wrinkles will be presented to highlight the potential of physics-based modelling of skin in tackling global challenges such as ageing of the population and the associated skin degradation, diseases and traumas.

Besides the brain, no other organ of the human body plays such a central role in our everyday biological *and* social life than the skin. The skin is the first line of defence of our body against the external environment, and therefore, acts as a primary and essential physical interface. This interface controls many types of exchanges between our inner and outside worlds which take the form of mechanical, thermal, biological, chemical and electromagnetic processes. These processes typically do not operate in isolation but are rather parts of a very dynamic system featuring complex nonlinear feedback mechanisms [1,2]. At a secondary, but nonetheless very important exchange level, the skin is not only a sort of social sign post that tells a story about us, but also a fundamental active social instrument that is pivotal in how we socially interact with our fellow humans, through conscious and unconscious communication cues [3]. The skin is also a biochemical plant that synthesizes vital compounds like vitamin D, hosts vital immunologic biochemical and cellular processes, and contains a rich sensory biophysical network that informs us in real-time of any haptic cues or potentially threatening physical insults and noxious agents [1,3].

Considering the place of the skin in our life and its multiple physiological functions, from genesis to death, understanding its complex physiology and biophysics in health, disease and trauma has become, particularly in the last two decades, a broad research arena. Interdisciplinary scientists, engineers and clinicians join forces to improve health, quality of life, design products that are easier and more comfortable to use or help us to slow down the signs of ageing by developing innovative surgical cosmetic procedures and skin care pharmaceutical compounds.

As in many branches of life sciences, be it in academia or industry, physics-based modelling of the skin has become an integral part of research and development. Nowadays, advanced numerical simulations, typically relying on finite element [4] and/or meta-modelling techniques [5], are used in the rational design of products intended to interact with the skin, from razors and incontinence nappies through respiratory masks and medical surfaces to running shoes, vehicle interiors and tactile electronic surfaces. At a more fundamental level, and as hypothesis-driven research tools, mathematical and computational models of the integument are developed to shed light on the biophysical complexity of skin physiology and to unravel particular mechanobiological aspects associated with diseases and the ageing process. This focus is particularly relevant for medical, pharmaceutical and cosmetic sciences.

Jor *et al.* [6] provided an excellent review discussing the current (as of 2013) and future challenges of the computational modelling and experimental characterization of skin mechanics. These authors highlighted the need for tight integration of modelling, instrumentation and imaging. They also recommended that focus should be directed towards the development of *in vivo* quantitative characterization techniques so that research could be more seamlessly translated into the clinical setting. Recently, Li [7] conducted a short review of modelling approaches for the description of *in vitro* biomechanical properties of the skin. He concluded that major research efforts should be devoted to the development of constitutive models of skin capable of accounting for viscoelastic effects, damage and fracture. Although Jor *et al.* [6] and Li [7] reviewed some popular constitutive models of the skin, no attempt was made to provide an extensive, detailed and *unified* review of formulations capable of representing the finite strain and anisotropic behaviour of skin, from elasticity through anelasticity to damage and growth.

Here, it is proposed to address this current shortcoming. The objective of this paper is to provide *a* review on some aspects of the mathematical and computational modelling of skin biophysics with special focus on theories based on the powerful constitutive framework offered by nonlinear continuum mechanics [8–10]. As in any review paper, there is a natural bias towards the topics covered which stems from the author's personal research but, here, it is hoped that the treatment of the subject is sufficiently general and high level to appeal to a broad range of scientists and engineers, particularly those involved in interdisciplinary skin research. The ambition is to present *a* selection of state-of-the-art skin models that are presently used in academic and industrial research, inform about the attractive prospects offered by biophysical modelling, highlight less well-known areas where the modelling of skin can address important scientific questions and practical engineering challenges. Ultimately, the aim is to stimulate cross-fertilizing of ideas and techniques, to encourage researchers from diverse scientific background to engage dialogue within their communities and outside, particularly with clinicians and biologists. More generally, the review should also be of interest to researchers involved in the modelling of biological soft tissues.

In mammalian species, the skin is part of the integumentary system that includes the skin itself as the largest organ of the body in terms of surface, and various appendages such as hair, nails and hooves, nerve receptors and glands. In humans, the skin—which can be considered as a membrane—accounts for up to 16% of an adult's total body weight while covering an average surface area of about 1.6m^{2} [11]. As alluded to in Introduction section, the skin is a very complex biological system featuring a multitude of coupled physical processes acting in concert, or sequentially, as for instance, in the case of wound healing where a cascade of biochemical and mechanobiological events is triggered by an injury [12]. From the mechanical and material science point of view, the skin is primarily a multi-phasic and multi-scale structure which, as a result, encompasses a rich set of mechanical properties and constitutive behaviours [13,14]. The biological nature of this structure renders these properties very dynamic, particularly over a lifetime, and like most biological tissues, there is a strong variability according to body site, individuals, age, sex, ethnicity and exposure to specific environmental conditions [15,16]. At the meso-macroscopic level, skin is generally considered as a multi-layer assembly made up of three main distinct structures: the epidermis, dermis and hypodermis (figure 1) [11,14,17,18].

Histological section of haematoxylin and eosin stained back human skin sample obtained from a 30-year-old healthy White Caucasian female volunteer following biopsy (10 × magnification, image resolution: 1600×1200 pixels, **...**

The epidermis—which is avascular—is a terminally differentiated stratified squamous epithelium about 200µm thick. 95% of the cells contained in the epidermis are keratinocytes which undergo mitosis at the 0.5–1µm thick epidermal basement membrane (also known as *basal lamina*) [19]. From this basement membrane which separates the epidermis from the dermis, the cells subsequently migrate from towards the skin surface, forming four main well delineated layers during their transit, namely the *stratum basale* (also called the *stratum germinativum*), the *stratum spinosum*, the *stratum granulosum* and *the stratum corneum*. This latter layer consists of a one to three cell-thick layer of dead keratinocytes representing a thickness of 15–30µm. The epidermis includes other cells such as melanocytes, Langerhans's cells and Merkel cells [11]. The complex formed by the living epidermis at the exclusion of the *stratum corneum* is called the living or viable epidermis.

The *stratum corneum* is the prime line of defence against external threats be they mechanical, thermal, chemical, electromagnetic or biological. In particular, the mechanical properties of the *stratum corneum* are key in conditioning transmission of loads and subsequent deformations of the other underlying skin layers across several spatial scales [20,21]. These purely mechanical aspects are essential for certain biophysical processes such as the stimulation of mechano-receptors that transduce mechanical energy into neural signalling (e.g. tactile perception [22] and pain [23]) or mechanobiological transduction involved in metabolic processes (e.g. homoeostatic regulation of the skin barrier function [24]).

Any variation in the mechanical properties of the *stratum corneum*—like those occurring daily as a results of fluctuations in internal and external environmental conditions (e.g. relative humidity, temperature) [25]—is likely to affect the material mechanical response and also the subsequent altered external surface topography. This has obvious consequences for the tribological response of the skin [21].

The living epidermis is connected to the underlying dermis through a three-dimensional (3D) interlocking wavy interface, called the dermal–epidermal junction (DEJ) which is the *basal lamina* evoked earlier. Papillae are the protrusions of the papillary dermis into the epidermis (figure 1). These finger-like structures increase the contact surface area between the reticular dermis and the living epidermis, and are thus believed to favour biochemical mass exchanges between these layers e.g. transport of oxygen or nutrients [19,26,27]. The DEJ controls the transit of biomolecules between the dermis and epidermis according to their dimension and charge. It allows the passage of migrating and invading cells under normal (i.e. melanocytes and Langerhans cells) or pathological (i.e. lymphocytes and tumour cells) conditions [27]. The behaviour of keratinocytes is influenced by the DEJ via modulation of cell polarity, proliferation, migration and differentiation.

As principally constituted of a dense array of crimped stiff collagen fibres, the dermis is the main load-bearing structural component of the skin when subjected to tension-inducing loads (e.g. in-plane tension, out-of-plane indentation and suction). It is 15–40 times as thick as the epidermis [11] and much thinner on the eyelids than on the back. The dermis can be decomposed into three main layers: the papillary layer juxtaposed to the epidermis, the sub-papillary layer underneath and the reticular layer which is connected to the underlying subcutaneous tissue. The papillary layer is defined by the rete ridges protruding into the epidermis and contains thin collagen fibres, a rich network of blood capillaries, sensory nerve endings and cytoplasms. The sub-papillary layer which is the zone below the epidermis and papillary layer features similar structural and biological components to the papillary layer. Besides a dominant content of types I and III collagen (respectively, 80 and 15% of total collagen content), the reticular layer is innervated and vascularized, contains elastic fibres (e.g. elastin) and the dermal matrix made of cells in the interstitial space. Cells present in the reticular dermis include fibroblasts, plasma cells, macrophages and mast cells.

Collagen fibres account for approximately 70% of the weight of dry dermis. Collagen fibres in the papillary and sub-papillary dermis are thin (because of their low aggregate content of fibrils) and sparsely distributed while reticular fibres are thick, organized in bundles and densely distributed. Fibrils are typically very long, 100–500nm in diameter featuring a cross striation pattern with a 60–70nm spatial periodicity. The diameter of thick collagen bundles can span 2–15µm. Birefringence techniques are promising tools to characterize the orientation and supramolecular organization of collagen bundles in skin [28].

The dermal matrix is composed of an extracellular matrix, ground substance which is a gel-like amorphous phase mainly constituted of proteoglycans and glycoproteins (e.g. fibronectin) as well as blood and lymph-derived fluids which are involved in the transport of substances crucial to cellular and metabolic activities. Proteoglycans are composed of multiple glycosaminoglycans (i.e. mucopolysaccharides) interlaced with back bone proteins. Dermal fibroblasts produced glycosamine which is rich in hyaluronic acid and therefore plays an essential role in moisture retention.

Unlike collagen fibres, elastic fibres are extremely elastic and can fully recover from strains in excess of 100% [29]. Their diameter ranges between 1 and 3µm. Their mechanical entanglement with the collagen network of the dermis is what gives the skin its resilience and recoil ability. This is evidenced by the correlation between degradation of elastin/abnormal collagen synthesis associated with ageing and the apparent stiffening of the dermis [30]. The diameter of elastic fibres in the dermis is inversely proportional to their proximity to the papillary layer where they tend to align perpendicular to the DEJ surface.

The subcutaneous tissue is the layer between the dermis and the fascia which is a band of connective tissue, primarily collagen, that attaches, stabilizes, encloses, and separates muscles and other internal organs. The thickness of subcutaneous tissue is highly variable intra- and inter-individually. This layer is mainly composed of adipocytes. Its role is to provide mechanical cushioning, heat generation and insulation as well as a reserve of nutrients.

A very important aspect of skin mechanics is that, *in vivo*, skin is in a state of complex in-plane heterogeneous tension patterns which depend on individuals, their age, body location and position. This was first evidenced by French anatomist and military surgeon Baron Guillaume Dupuytren (1777–1835) who, incidentally, also gave his name to a proliferative connective tissue disorder, the Dupuytren's contracture. In 1834, Dupuytren examined the chest of a man who had stabbed himself with a round-tipped awl in a failed suicide attempt, and noted that the inflicted wounds were elliptical and not circular [31–33]. Later, in 1838, French surgeon and medical historian Joseph-François Malgaigne (1806–1865) reported that this type of elliptical wound patterns were oriented differently according to body location. Although these two surgeons had discovered the existence of tension lines *in vivo*, it was the Austrian anatomist Karl Langer (1819–1887) [31–33] who provided a sound explanation and comprehensive anatomical basis by conducting a large number of punctures on cadavers, using a round awl, and meticulously mapping the direction and dimensions of the created elliptical wounds on the body. Because these experiments were conducted on cadavers with extremities in extension, Langer observed that tension lines varied with the position of the cadaver. Real *in vivo* conditions are also likely to differ from cadaveric conditions because of *rigor mortis* which stiffens muscles and joints. Moreover, tension lines are also dependent upon skin dysfunctions and diseases such as the Elhers–Danlos syndrome.

In the light of the skin anatomy described above, and, as will be discussed in the next section in more details, depending on the spatial scale considered, the skin can be viewed as a structure or material (i.e. a homogeneous structure). Here, in this brief review of mechanical properties, we will mainly focus on the latter. Owing to its characteristic dimensions and interfacial nature, the skin is often described as a membrane in the general sense. From the mechanical point of view, it would be more accurate to describe it as a multi-layer shell structure because it possesses a non-negligible bending stiffness despite mainly sustaining membrane strains.

Owing to its significant thickness compared with that of the *stratum corneum* and viable epidermis, combined with its high content in collagen fibres, the dermis is the main contributor to the tensile mechanical properties of the skin. The intricate fibrous collagen architecture of the dermis and its close mechanical connectivity to elastin fibres lead to anisotropic and nonlinear macroscopic mechanical properties which is consistent with the strain stiffening effect typically observed in biological soft tissues under tension [8].

Here, it is important to point out and recognize that the *measured* mechanical anisotropy of the skin is not only due to the structural characteristics of the skin layers and their microstructural constituents but is also the result of its mechanical interplay with the Langer lines. In-plane anisotropy of the skin is correlated with the distribution and orientation of Langer lines [34] while out-of-plane (or across-the-thickness) anisotropy is due to the distinct mechanical properties and complex 3D structure of the skin layers. This presents a number of challenges for the experimental characterization of the mechanical properties of skin [6,35] as well as their mathematical and computational realization. Intuitively, one would expect that the inclusion or not of pre-stress/pre-strain in a mathematical or computational model of the skin would significantly influence its mechanical response. Flynn *et al*. [36,37] determined a relation between the *in vivo* relaxed skin tension lines on a human face and the directional dependency of the skin stiffness using a combination of contact measurement techniques and inverse finite-element methods. These authors demonstrated the need to account for these tension lines in the characterization of the anisotropic properties of the human skin. Recently, Deroy *et al.* [38] developed a non-destructive and non-invasive experimental protocol based on elastic surface wave propagation to determine the orientation of skin tension lines. The method was validated on canine cadaveric specimens. The Cutiscan CS 100® (Courage-Khazaka Electronic GmbH, Köln, Germany) is a device offering the possibility of qualitatively measuring mechanical anisotropy of the skin *in vivo* [39].

A variety of methods have been used to measure the mechanical properties of skin: e.g. uniaxial and biaxial tensile tests [34,40–42], multi-axial tests [13,43], application of torsion loads [44], indentation [45], suction [39,46–49] and bulge testing [50]. More recent techniques have focused on the experimental characterization of the mechanical properties of the epidermis [51–53] which are particularly relevant for cosmetic and pharmaceutical applications. Beside intra- and inter-individual biological variability as well as sensitivity to environmental conditions, the nature of these experimental techniques combined with their operating spatial scale (e.g. macroscopic or cellular levels) are the main reasons for such a wide variability in mechanical properties reported in the literature. Differences for the mechanical properties of the skin or those of its individual layers can span several orders of magnitude ( Table 1 in [20]).

Given its complex hierarchical structure, and like most soft tissues, the skin exhibits a wide range of viscoelastic phenomena including creep, relaxation, hysteresis [16,40,54,55] and strain-rate dependency [42].

The review of the structural and macroscopic mechanical properties of the skin in the previous section has highlighted three main characteristics, namely, nonlinear behaviour, structural non-homogeneity and mechanical anisotropy. The first task in any constitutive modelling endeavour is to assess the ‘why’ of the constitutive model. Why do we need a model? What are the answers we are looking for? What are the hypotheses we want to assess? Do we want to be *predictive* or simply *descriptive*? What are the operating conditions we need to account for? Can we practically design a physical experiment or a series of experiments that will inform the mathematical formulation of the constitutive equations and allow us to extract constitutive parameters so we can exploit the constitutive model?

The formulation of any constitutive model must be designed in the light of the intended application, which means that the model will have to rely on a number of (simplifying) assumptions that will discard certain aspects of the physical system (e.g. heterogeneous microstructure, time-dependent behaviour, large strains) while accounting for others. For example, for studying the transport of drug through the skin it would be highly relevant to account for the geometry and porous nature of the microstructure of the relevant skin layers. If biochemical aspects such as molecular water binding to collagen were of interest, then models and techniques capable of representing these atomic/molecular scale phenomena would be appropriate. In these conditions, using a finite strain model would be ‘overkill’ as only very small strains would be expected. If one were to investigate the failure of skin in traumatic situations (e.g. vehicle accident) it would be essential to avail oneself with a model capable of describing the large strain and failure behaviour of the tissue.

Mathematical and computational models of the skin can be classified into three main categories: phenomenological, structural and structurally based phenomenological models.

Phenomenological models are based on the assumption that the skin is a homogeneous material where the microstructure and multiple phases as well their associated mechanical properties are ignored These phenomenological models aim to capture the overall—generally macroscopic—behaviour of the tissue without accounting for the individual behaviour of its elemental constituents and their mutual interactions [56]. Typically, if one considers mechanical behaviour only, a phenomenological model is a set of mathematical relations that describe the evolution of stress as a function of strain [8]. Provided the formulation is appropriate, it is generally always possible to fit such a constitutive law to a set of experimental data. However, the main drawback of this approach is that the resulting constitutive parameters often do not have a direct physical interpretation, and the model effectively acts as a ‘black box’ without the flexibility of integrating and studying mechanistic structural effects.

Structural models of the skin represent the tissue as a composite material made of an *explicitly* defined assembly of key microstructural elements (e.g. collagen fibres arranged in bundles with a certain degree of crimp and dispersion, within a matrix mainly composed of proteoglycans). The way these structural elements interact can also be specified by developing appropriate equations (e.g. mutual shear interactions of collagen fibrils and fibres or small-range electromagnetic interactions between proteoglycans and collagen fibres [57]).

In this approach, not only the mechanical properties of the individual basic structures need to be determined or known but also the way they are geometrically arranged to form the macroscopic tissue, and how they interact mechanically, thermally or through any other type of physics. The overall mechanical properties of the tissue are the result of this—generally nonlinear—coupling between geometry and mechanics. Structural models can be viewed as geometrical assemblies of phenomenological models. For this reason, one could argue that, strictly speaking, structural models are not fully structural, and that classification as a structural or phenomenological model is a matter of spatial scale. Only if models are built *ab initio* from the first principles of quantum chemistry could one talk about structural models. In that respect, the constitutive parameters of phenomenological models describing the behaviour of elemental microstructural components may also not have a direct physical interpretation but, typically, they do. For example, if collagen fibres are part of the formulation, their elastic modulus, waviness and degree of dispersion around a main orientation are indeed constitutive parameters with have a direct physical interpretation. One of the main drawbacks of structural models lays in the necessity to have accurate information about the geometrical and material characteristics of each elemental building block as well as their mutual spatial and interfacial arrangement. This presents obvious experimental characterization challenges which, however, are progressively overcome as technologies in this field improve, and new techniques emerge. Moreover, in a computational finite-element environment, structural models spanning several orders of magnitude in terms of length scale require significantly high mesh density to explicitly capture the geometry of the structural constituents. Inevitably, despite tremendous advances in the computing power of modern processors, this can lead to prohibitively computationally expensive analyses and lengthy run times.

A judicious compromise between strictly phenomenological and structural models is a third-class of models, formulated by combining certain characteristics of phenomenological models to those of structural ones. These models could be denoted under the general appellation of *structurally based phenomenological models* or *structurally based continuum models*. The continuum composite approach of Spencer is a good example of that idea [58]. In this type of constitutive formulation for fibre-reinforced composite materials, fibres are not explicitly modelled at the geometric level but their mechanical contribution to the overall (i.e. continuum) behaviour is implicitly accounted for by strain energy density terms directly related to microstructural metrics (e.g. stretch along the fibre direction). Nowadays, models adopting this type of constitutive hypotheses dominate the published literature as they have been very successful at representing the biomechanics of many biological tissues [9,10]. The mathematical formulation of these models can be viewed as multi-scale in the general sense as it accounts for structural geometrical/mechanical features/effects arising from a smaller spatial scale than that at which the gross physical response is represented in a continuum sense. Multi-scale formulations in the homogenization theory/computational mechanics sense are also possible but require more advanced finite-element formulations for practical calculations [59,60].

It appears that over the last 15 years there has been a significant thrust in the biomechanics community to delve into this under-researched topic that is skin biophysics. This possibly stems from the academic potential to deliver novel research with high impact which is intimately correlated to the drive to address growing research needs in many industrial sectors such as personal care and consumer products, cosmetics and biomedical devices. Like what happened in the communities of orthopaedic, dental and cardiovascular clinicians, the new generations of reconstructive surgeons and dermatologists are progressively acknowledging and embracing the potential and key role that computational modelling can play in improving and optimizing treatment procedures. An attractive feature of computational models is the ability to quantify physical parameters required to obtain a specific outcome. For example, Prof. Kuhl's group at Stanford University has pioneered the development of mechanobiological adaptation models for skin to optimize the outcome of reconstructive surgery procedures in children [61–67] (see electronic supplementary material, §S1.10).

Given the hierarchical nature of the structure of biological tissues, especially the skin, capturing certain key aspects of their microstructure in computational models opens up the possibility of conducting hypothesis-driven research about the links between their structure and function [8] in a very systematic and mechanistic way. One can isolate and study the influence of certain microstructural and/or material parameters on the global response of the system [20,21,68] to formulate hypotheses and theories that can further our fundamental understanding of skin physiology and also assist in the design of physical experiments.

According to the classification of models presented in the previous section, image-based models could be best described as computational (micro)structural models where the constitutive model used for representing the material behaviour of each substructure could be a phenomenological, structural or structurally based phenomenological model. Image-based modelling has become ubiquitous in many research domains where capturing complex multi-phasic structures at various length scales is essential for conducting physics-based numerical analyses [69]. In the context of biological tissues and structures, any type of imaging modalities (1D, 2D, 3D, 4D) can be used including digital optical photography [20], computed tomography [70], magnetic resonance imaging [71] and laser confocal microscopy [72].

Developing an image-based computational model is the process of acquiring the geometry of a physical system through a single or a combination of imaging modalities, bring it into the digital world under the form of a picture or a series of pictures, conduct image processing on these files to identify and isolate regions of interest (i.e. segmentation) (e.g. distinct skin layers) and mesh these digital regions into finite elements so that the partial differential equations governing the physics considered could be solved using the finite-element method (figure 2).

Image-based microstructural finite-element models of the skin can provide unique insights into the interplay between material and structural properties of the skin even when considering 2D geometries. Leyva-Mendivil *et al*. [20] unveiled and quantified structural folding mechanisms driven by the mechanics of the *stratum corneum* in combination with its intricate surface topography when the skin is subjected to plane-strain compression (figure 3). It was shown that *macroscopic* strains may significantly be modulated at the *microscopic* (i.e. local) level. The convoluted shape of the *stratum corneum* can effectively acts as a strain reducer. These authors subsequently used a similar approach to show that the *macroscopic* coefficient of friction between the skin and a rigid slider moving across its surface is noticeably higher that the *local* coefficient of friction applied as an input parameter to the finite-element analyses [21].

This demonstrated that the deformation-induced component of skin friction is significant unlike what has been widely assumed in the tribology community [73–77] where adhesion-induced friction is deemed to be the dominant contributor to macroscopic friction. Similar observations were made in a similar computational contact homogenization study by Stupkiewicz *et al.* [78]: geometrical effects alone can have a significant impact on the macroscopic frictional response of elastic contacts. To date, despite much experimental and modelling studies investigating shear stress at the surface of the skin in relation to skin injuries and pressure ulcers [51,53,79], very little efforts have been made to develop methodologies to gain a more accurate and mechanistic understanding of how shear stresses are induced at the level of skin microrelief asperities and how they propagate from the skin surface to the deeper layers where they are likely to mechanically stress living cells. Ultimately, excessive stress or strain can lead to cell damage and death, which, at a meso-macroscopic level translates into tissue damage and loss of structural integrity. An image-based finite-element study was recently conducted by Leyva-Mendivil *et al*. [68] to investigate these aspects and demonstrated again the importance of accounting for skin microstructure in this type of study.

Here we present what are considered state-of-the-art mechanical constitutive models of the skin. The focus is on nonlinear elasticity and viscoelasticity theories. In the electronic supplementary material section, plastic effects [80], softening and damage [81], recent coupled thermomechanical [82] as well as mechanobiological models for growth [64] will be presented. Other types of constitutive laws such as chemo-mechanobiological models of wound healing [12,83–85] are beyond the scope of the present review.

Nonlinear continuum mechanics provides a flexible mathematical framework to model a wide range of constitutive behaviours [86,87] from nonlinear elasticity and finite strain plasticity to thermoelasticity and biological growth.

Accounting for finite deformations is essential for many applications in skin biophysics, and it is therefore sensible to develop constitutive equations *a priori* valid for finite kinematics. To this end, in this section, essential definitions of kinematic and kinetic quantities are provided, and, without loss of generality, particularized for Cartesian coordinate systems. Because of the popularity of their use in the constitutive modelling of biological tissues [9], particularly those featuring anisotropic mechanical properties, tensor invariants and structural tensors are also introduced. Invariant formulations typically postulate the existence of a strain energy function depending on a set of tensorial invariants of a given deformation or strain measure [87–90]. Introducing a dependency of the strain energy to tensor agencies characterizing local microstructural features (e.g. local collagen fibre bundle orientation) is a very efficient way to incorporate microstructural information within a macroscopic continuum constitutive formulation. This class of material models corresponds to that of the microstructurally based models defined in the previous section. In order to ensure a direct physical interpretation of constitutive parameters it is judicious to select a set of tensor invariants that characterize particular deformation modes the tissue is known to be subjected to and that also can be physically measured [91,92]. For collagenous tissues, a standard assumption is to consider the tissue as a *continuum* composite material made of one or several families of (oriented) collagen fibres embedded in a highly compliant isotropic solid matrix composed mainly of proteoglycans. A preferred fibre alignment is defined by the introduction of a so-called *structural* or *structure tensor* which appears as an argument of the strain energy function [88,90]. These notions are detailed in the next sections.

A key basic kinematic entity in continuum mechanics is represented by the deformation gradient **F** defined as:

$$\mathbf{\text{F}}(\mathbf{\text{X}}):=\frac{\mathrm{\partial}\phi (\mathbf{\text{X}})}{\mathrm{\partial}\mathbf{\text{X}}}=\sum _{i,I=1}^{3}\frac{\mathrm{\partial}{\phi}_{i}}{\mathrm{\partial}{\mathrm{X}}_{I}}{\mathbf{\text{e}}}_{\mathbf{\text{i}}}\otimes {\mathbf{\text{E}}}_{\mathbf{\text{I}}},$$

4.1

where **X** is the position of a material point in the *Lagrangian*—or *reference*—configuration, while **x**=(**X**) is its material placement in the *Eulerian*—or *current*—configuration. ${\{{\mathbf{\text{E}}}_{I}\}}_{I=1,2,3}$ and ${\{{\mathbf{\text{e}}}_{i}\}}_{i=1,2,3}$ are fixed orthonormal bases in the *Lagrangian* and *Eulerian* configurations, respectively. ‘.’, ‘:’, $\otimes $ and ‘T’, respectively, denote the scalar product of second-order tensors, double-contraction tensor product, tensor outer product and transpose operators. **F** maps infinitesimal line vectors from the material to spatial configuration while cofactor(**F**)=*J*·**F**^{−T} and *J*=determinant(**F**), respectively, maps oriented infinitesimal surface and volume. The right (material) and left (spatial) Cauchy–Green deformation tensors are, respectively, defined as **C**=**F**^{T}·**F** and **b**=**F**·**F**^{T}. Because these two tensors only contain information about change of length, and therefore exclude local material rotations, they are appropriate to define valid constitutive equations that satisfy the Principle of Material Frame Indifference [87]. From the deformation tensor **C**, one can define the Green–Lagrange strain tensor **E**=(**C**−**I**)/2, where **I** is the second-order identity tensor. The classical three principal deformation invariants of **C** which define the isotropic (or non-directional) response of a given material are defined as follows [86]:

$${I}_{1}={\textstyle \text{trace}}(\mathbf{\text{C}})=\mathbf{\text{C}}:\mathbf{\text{I}};\phantom{\rule{1em}{0ex}}{I}_{2}=\frac{1}{2}[{I}_{1}-{\mathbf{\text{C}}}^{2}:\mathbf{\text{I}}]\phantom{\rule{1em}{0ex}}{\textstyle \text{and}}\phantom{\rule{1em}{0ex}}{I}_{3}={\textstyle \text{determinant}}(\mathbf{\text{C}}).$$

4.2

To characterize a general orthotropic symmetry, one can introduce three unit vectors ${\mathit{n}}_{0}^{i}$ associated, respectively, with the principal material directions *i* in the reference configuration. Because these material directions are *signed* directional properties one introduces the concept of structural tensors ${\mathfrak{B}}_{0}^{i}$ [87–90] which are even functions of these unit vectors:

$${\mathfrak{B}}_{0}^{i}={\mathit{n}}_{0}^{i}\otimes {\mathit{n}}_{0}^{i}\phantom{\rule{1em}{0ex}}\{i=1,2,3\}({\textstyle \text{no summation on}}i).$$

4.3

The spatial counterparts of these structural tensors are defined as:

$${\mathfrak{B}}_{t}^{i}=\frac{1}{{\lambda}_{\hspace{0.17em}i}^{2}}\mathbf{\text{F}}\cdot {\mathfrak{B}}_{\hspace{0.17em}0}^{i}\cdot {\mathbf{\text{F}}}^{\mathrm{T}}=\frac{1}{{\lambda}_{\hspace{0.17em}i}^{2}}\mathbf{\text{F}}\cdot ({\mathit{n}}_{0}^{i}\otimes {\mathit{n}}_{0}^{i})\cdot {\mathbf{\text{F}}}^{\mathrm{T}}={\mathit{n}}_{t}^{i}\otimes {\mathit{n}}_{t}^{i}\phantom{\rule{1em}{0ex}}\{i=1,2,3\}({\textstyle \text{no summation on}}i),$$

4.4

where ${\mathit{n}}_{t}^{i}\hspace{0.17em}\{i=1,2,3\}$ are the unit vectors in the spatial configuration and *λ _{i}* {

For material *isotropy* all material directions are equivalent so that:

$${\mathfrak{B}}_{0}^{1}={\mathfrak{B}}_{0}^{2}={\mathfrak{B}}_{0}^{3}={\mathit{n}}_{0}^{1}\otimes {\mathit{n}}_{0}^{1}={\mathit{n}}_{0}^{2}\otimes {\mathit{n}}_{0}^{2}={\mathit{n}}_{0}^{3}\otimes {\mathit{n}}_{0}^{3}=\mathbf{\text{1}}\otimes \mathbf{\text{1}}.$$

4.5

For *transverse isotropy*, if the preferred material direction is given by ${\mathit{n}}_{0}^{3}$:

$${\mathfrak{B}}_{0}^{3}={\mathit{n}}_{0}^{3}\otimes {\mathit{n}}_{0}^{3}\phantom{\rule{1em}{0ex}}{\textstyle \text{and}}\phantom{\rule{1em}{0ex}}{\mathfrak{B}}_{0}^{1}={\mathfrak{B}}_{0}^{2}={\textstyle \frac{1}{2}}(\mathbf{\text{I}}-{\mathfrak{B}}_{0}^{3}).$$

4.6

For *orthotropic symmetry*, the preferred material directions are orthogonal to each other.

In the context of soft tissue mechanics, depending on the spatial scale of investigation, a single unit vector ${\mathit{n}}_{0}^{i}$ can represent the local orientation of a single collagen microfibril or fibril, that of a single fibre bundle or that of a family of fibres. This gives rise to a *transverse isotropy* symmetry if one assumes that these fibres are embedded in an isotropic matrix. Two additional tensorial invariants characteristic of transverse isotropy symmetry can be defined [87,88,90]:

$${I}_{4}^{i}=\mathbf{\text{C}}:{\mathfrak{B}}_{0}^{i}={\lambda}_{\hspace{0.17em}i}^{2}\phantom{\rule{1em}{0ex}}{\textstyle \text{and}}\phantom{\rule{1em}{0ex}}{I}_{5}^{i}={\mathbf{\text{C}}}^{2}:{\mathfrak{B}}_{0}^{i},$$

4.7

where ${I}_{4}^{i}$ is a very convenient invariant with a direct physical interpretation as it is the square of the stretch along the corresponding fibre direction. It is worth pointing out that **C** and **b** have the same invariants.

As a general procedure to formulate constitutive equations for hyperelastic materials, one can postulate the existence of a strain energy density *ψ*, *isotropic* function of its deformation or strain invariant arguments. Stress and elasticity tensors are obtained by, respectively, first- and second-order differentiation [87]. Hyperelastic formulations can serve as a basis for inelastic formulations such as finite strain plasticity, finite strain viscoelasticity and growth, or any combination of these constitutive behaviours.

The Lagrangian second Piola–Kirchhoff stress tensor, **S**, is readily obtained by differentiation of the strain energy density function with respect to the right Cauchy–Green deformation tensor while applying the chain rule of differentiation for the *n* deformation invariants *I _{i}*:

$$\mathbf{\text{S}}=2\sum _{i=1}^{n\hspace{0.17em}\mathrm{invariants}}\left(\frac{\mathrm{\partial}\psi}{\mathrm{\partial}{I}_{i}}\frac{\mathrm{\partial}{I}_{i}}{\mathrm{\partial}\mathbf{\text{C}}}\right).$$

4.8

The Cauchy stress tensor, often referred as *true* stress tensor, is obtained by push-forward operation of the second Piola–Kirchhoff stress tensor **S** from the reference to the current configuration [87]:

$$\mathit{\sigma}=\frac{1}{J}({\phi}_{\ast}\mathbf{\text{S}})=\frac{1}{J}\mathbf{\text{F}}\cdot \mathbf{\text{S}}\cdot {\mathbf{\text{F}}}^{\mathrm{T}}=\frac{2}{J}\mathbf{\text{F}}\cdot \left(\frac{\mathrm{\partial}\psi}{\mathrm{\partial}\mathbf{\text{C}}}\right)\cdot {\mathbf{\text{F}}}^{\mathrm{T}}=\frac{2}{J}\sum _{i=1}^{n\hspace{0.17em}\mathrm{invariants}}\left(\frac{\mathrm{\partial}\psi}{\mathrm{\partial}{I}_{i}}\mathbf{\text{F}}\cdot \frac{\mathrm{\partial}{I}_{i}}{\mathrm{\partial}\mathbf{\text{C}}}\cdot {\mathbf{\text{F}}}^{\mathrm{T}}\right),$$

4.9

while the (volume ratio)-weighted Cauchy stress is the Kirchhoff stress tensor:

$$\mathit{\tau}=J\mathit{\sigma}=2\frac{\mathrm{\partial}\psi}{\mathrm{\partial}\mathbf{\text{b}}}\cdot \mathbf{\text{b}}=2\sum _{i=1}^{n\hspace{0.17em}\mathrm{invariants}}\left(\frac{\mathrm{\partial}\psi}{\mathrm{\partial}{I}_{i}}\frac{\mathrm{\partial}{I}_{i}}{\mathrm{\partial}\mathbf{\text{b}}}\right)\cdot \mathbf{\text{b}}.$$

4.10

Remarks about the associated material and spatial elasticity tensors in the context of finite-element procedures are provided in electronic supplementary material, §S1.1. The particular form of the tangent stiffness required for the implementation of constitutive models in the commercial software package Abaqus/Implicit (Simulia, Dassault Systèmes, Providence, RI, USA) is also derived.

In this section, a selection of some of the seminal and/or current state-of-the-art nonlinear elastic models of skin are presented. As mentioned earlier in this review, the anisotropic mechanical properties of the skin are critical for most applications, so the focus of this review is, therefore, on constitutive models capturing these characteristics.

A very significant microstructurally based continuum model of skin was first proposed by Lanir [93] who assumed that skin was a continuum composite made of an isotropic ground substance matrix in which oriented collagen and elastin fibres were embedded. These fibres obeyed an angular distribution ${\mathcal{R}}_{k}(\theta )$ (see electronic supplementary material, §S1.2). The total strain energy function of the material is defined as:

$$\psi =(1-{\mathcal{V}}_{\mathrm{c}}-{\mathcal{V}}_{\mathrm{e}})\frac{\mu}{2}({I}_{1}-3)+\sum _{k}^{\hspace{0.17em}}{\mathcal{V}}_{k}\sum _{\theta}^{\hspace{0.17em}}{\mathcal{R}}_{k}(\theta )\frac{{K}_{k}}{2}{(\lambda -1)}^{2},$$

4.11

where *μ* is a shear modulus defining a neo-Hookean hyperelastic potential for the matrix, ${\mathcal{V}}_{k}$, ${K}_{k}$ and ${\mathcal{R}}_{k}$ are, respectively, the volume fraction, fibre stiffness and fibre orientation probability density function of each fibrous phase *k* (collagen: *k*=c and elastin: *k*=e).

Collagen and elastin fibres are assumed to be linear elastic and unable to sustain compression along their axis:

$$\begin{array}{rl}{f}_{k}(\lambda )& \hspace{0.17em}={K}_{k}(\lambda -1),\phantom{\rule{1em}{0ex}}\lambda \ge 1\\ {\textstyle \text{and}}\phantom{\rule{2em}{0ex}}{f}_{k}(\lambda )& \hspace{0.17em}=0,\phantom{\rule{1em}{0ex}}{\textstyle \text{otherwise}}.\end{array}\}$$

4.12

In Lanir's model, collagen fibres are assumed to be undulated and can only resist loads when fully straightened. The degree of crimp of collagen fibres oriented in the direction *θ* is not uniform and follows a Gaussian distribution *D*(*x*), where *x* is the stretch required to straighten an uncrimped fibre. The Cauchy stress in skin is given as:

$$\mathit{\sigma}=\frac{1}{J}\sum _{k}{\mathcal{V}}_{k}{\int}_{0}^{\pi}\left[{\mathcal{R}}_{k}(\theta ){f}_{k}^{\ast}(\lambda )\mathbf{\text{F}}\frac{\mathrm{\partial}\lambda}{\mathrm{\partial}\mathbf{\text{E}}}{\mathbf{\text{F}}}^{\mathrm{T}}\right]\hspace{0.17em}{\textstyle \text{d}}\theta -p\mathbf{\text{1}},$$

4.13

where ${f}_{k}^{\ast}(\lambda )$ is the force per unit undeformed cross-sectional area of an individual fibre at stretch *λ* and *p* is a hydrostatic pressure that represents the mechanical contribution of the isotropic matrix. For collagen fibres ${f}_{k}^{\ast}(\lambda )$ is expressed as follows:

$${f}_{k=\mathrm{c}}^{\ast}(\lambda )={f}_{\mathrm{c}}^{\ast}(\lambda )={\int}_{1}^{\lambda}D(x){f}_{\mathrm{c}}\left(\frac{\lambda}{x}\right){\textstyle \text{d}}x.$$

4.14

As elastin fibres are assumed to be undulation-free, ${f}_{\mathrm{e}}^{\ast}(\lambda )={f}_{\mathrm{e}}(\lambda ).$

Although Lanir [93] did not identify its model with experimental data, Meijer *et al.* [94] exploited it by characterizing in-plane mechanical properties of human forearm skin using an hybrid numerical-experimental approach. Collagen fibre stiffness and mean undulation were estimated from the physical tests while the other constitutive parameters were extracted from the literature. Very low values of collagen fibre stiffness ranging from 51 to 86MPa were determined from the identification procedure. As an explanation, the partial mechanical recruitment of fibres for the specific limited strain range applied during the physical test was invoked. Lanir's model was used as a constitutive framework for an *in vitro* porcine skin model in a study by Jor *et al.* [95]. Constitutive parameters were determined through numerical optimization, and, as it is often the case in identification procedures, multiple comparable solutions were obtained. The mean orientation of collagen fibres was determined to lay within the 2–13 degrees range which is consistent with observations. Young's modulus ranges for the ground substance matrix and collagen fibres were, respectively, 5–12kPa and 48–366MPa. Fixing structural parameters that can be experimentally measured (e.g. fibre splay distribution) is a way to accelerate numerical identification and eliminate feasible but unrealistic solutions, thus reducing the issue of non-unicity in parameter sets.

Fibre dispersion can be accounted for by means of two main modelling approaches [96]. The first one, termed the ‘*angular integration approach*’ is due to Lanir [93] while the second approach, known as the ‘*generalized structure tensor*’ approach is due to Gasser *et al.* [97]. For a more detailed discussion on these topics and a new approximation of the π-periodic von Mises distribution, see electronic supplementary material, §S1.2.

This class of models is based on the notion of entropic elasticity and the micromechanics of macromolecule mechanical networks [98]. These concepts are briefly reviewed in electronic supplementary material, §S1.4. Bischoff *et al.* [99] extended the eight-chain model of Arruda & Boyce [98] developed for polymer elasticity to orthotropy by considering distinct values $a=\parallel \mathbf{\text{a}}\parallel $, $b=\parallel \mathbf{\text{b}}\parallel $and $c=\parallel \mathbf{\text{c}}\parallel $ for the three characteristic dimensions of the original cuboid microscopic unit cell defined by Arruda and Boyce (see electronic supplementary material, figure S8). Earlier, Bischoff *et al.* [100] had demonstrated that, despite being mechanically isotropic, the Arruda–Boyce model could reproduce load-induced anisotropy. However, the model could not capture the mechanical response of rabbit skin under biaxial extension [54]. Bischoff–Arruda–Grosh (BAG)'s strain energy function per unit volume is decomposed into chain and bulk energies:

$$\psi ={\psi}_{\mathrm{chain}}+{\psi}_{\mathrm{bulk}}$$

4.15

and

$${\psi}_{\mathrm{chain}}={\psi}_{0}+\frac{n\mathfrak{K}\theta}{4}\left\{N\sum _{i=1}^{4}\left[\frac{{\overline{r}}_{i}}{N}{\beta}_{i}+\mathrm{ln}\left(\frac{{\beta}_{i}}{\mathrm{sinh}{\beta}_{i}}\right)\right]-\frac{{\mathcal{L}}^{-1}(1/\sqrt{N})}{\sqrt{N}}\mathrm{ln}({\lambda}_{\mathbf{\text{a}}}^{{a}^{2}}{\lambda}_{\mathbf{\text{b}}}^{{b}^{2}}{\lambda}_{\mathbf{\text{c}}}^{{c}^{2}})\right\},$$

4.16

where *ψ*_{0} is a reference energy of the chain in the reference configuration, calculated so that the reference state is stress free, *n* is the volumetric chain density, ${\overline{r}}_{i}={r}_{i}/d$ is the end-to-end length of deformed chain normalized by the chain link length *d*, $\sqrt{N}$ is the length of each *undeformed* chain:

$$N={\textstyle \frac{1}{4}}({a}^{2}+{b}^{2}+{c}^{2}),$$

4.17

where *λ*** _{i}** are the stretches along the principal material axes (

$${\lambda}_{\mathbf{\text{i}}}=\sqrt{\mathbf{\text{C}}:(\mathbf{\text{i}}\otimes \mathbf{\text{i}})},\phantom{\rule{1em}{0ex}}(\mathbf{\text{i}}=\mathbf{\text{a}},\mathbf{\text{b}}{\textstyle \text{or}}\mathbf{\text{c}})$$

4.18

and

$${\beta}_{i}={\mathcal{L}}^{-1}\left(\frac{{\overline{r}}_{i}}{N}\right).$$

4.19

Finally, the bulk energy is defined as:

$${\psi}_{\mathrm{bulk}}=\frac{\kappa}{{\alpha}^{2}}\{\mathrm{cosh}[\alpha (J-1)]-1\},$$

4.20

where *κ* is the bulk modulus and *α*, a parameter that controls the shape of the pressure-*J* curve. Bischoff *et al.* [99] calibrated their constitutive model using data from biaxial testing on rabbit skin [54] and obtained an excellent fit with only seven parameters. The parameters were *n*=3.75×10^{22} (m^{−3}), *N*=1.25, *κ*=50 (kPa), *α*=1 and {*a*, *b*, *c*}={1.37, 1.015, 1.447}.

Flynn and McCormack developed a series of skin models based on BAG's formulation to investigate the wrinkling behaviour of skin [101], scar tissue contraction [102] and wrinkling/ageing of the skin [103]. Kuhl *et al.* [104] particularized BAG's model to the case of transverse isotropy by setting two of the unit cell dimensions equal and fitted their model to the ubiquitous rabbit skin data from Lanir & Fung [54] but the fit was not as good as that provided by the original BAG's model.

A polyconvex anisotropic strain energy function for soft tissues was developed by Limbert & Middleton [105] and independently formulated by Itskov *et al.* [106] shortly after. The constitutive framework was based on the generalized structural tensor invariant formulation proposed by Itskov & Aksel [107]. Rabbit skin biaxial tensile test data from Lanir & Fung [54] were used by Limbert and Middleton to identify sets of constitutive parameters. It was shown that a three-term series formulation was sufficient to obtain a very good fit between the experimental measurements and the predictions of the model. The starting point of the formulation which features novel invariants (compared to those described in §4.1.1) is the definition of a generalized structural tensor (indexed by *k*) as the weighted sum of three-mutually orthogonal structural tensors [106]:

$${\mathfrak{B}}_{0[k]}^{i}=\sum _{i=1}^{3}{w}_{k}^{i}{\mathfrak{B}}_{0}^{i}={w}_{k}^{1}{\mathfrak{B}}_{0}^{1}+{w}_{k}^{2}{\mathfrak{B}}_{0}^{2}+{w}_{k}^{3}{\mathfrak{B}}_{0}^{3},\phantom{\rule{1em}{0ex}}\{k=1,2,\dots ,n\},$$

4.21

where ${\omega}_{k}^{i}\hspace{0.17em}(i=1,2,3;\hspace{0.17em}\hspace{0.17em}k=1,2,\dots ,n)$ are non-negative scalars dependent on the principal directions. The generalized structural tensors must satisfy the normalization condition so that:

$$\text{trace}}({\mathfrak{B}}_{0[k]}^{i})=1\phantom{\rule{1em}{0ex}}{\textstyle \text{if}}\sum _{i=1}^{3}{\omega}_{k}^{i}=1.$$

4.22

Itskov & Aksel [107] proposed the following two sets of invariants to describe the generalized orthotropic behaviour of hyperelastic materials:

$${I}_{[k]}=\sum _{i=1}^{3}{\omega}_{k}^{i}({\mathfrak{B}}_{0[k]}^{i}:\mathbf{\text{C}})$$

4.23

and

$${J}_{[k]}={I}_{3}[{\omega}_{k}^{1}{\textstyle \text{trace}}({\mathbf{\text{C}}}^{-1}{\mathfrak{B}}_{0[k]}^{1})+{\omega}_{k}^{2}{\textstyle \text{trace}}({\mathbf{\text{C}}}^{-1}{\mathfrak{B}}_{0[k]}^{2})+{\omega}_{k}^{3}{\textstyle \text{trace}}({\mathbf{\text{C}}}^{-1}{\mathfrak{B}}_{0[k]}^{3})].$$

4.24

The original strain energy function proposed by Itskov & Aksel [107] was designed to model the mechanics of transversely isotropic calendered rubber sheets at high strains which it did very well:

$$\psi =\frac{1}{4}\sum _{k=1}^{n}{\mu}_{k}\left[\frac{1}{{\alpha}_{k}}({I}_{[k]}^{{\alpha}_{k}}-1)+\frac{1}{{\beta}_{k}}({J}_{[k]}^{{\beta}_{k}}-1)+\frac{1}{{\chi}_{k}}({{I}_{3}}^{-{\chi}_{k}}-1)\right].$$

4.25

This function was subsequently modified by Limbert & Middleton [105] to capture the typical exponential behaviour of the toe region of the stress–strain curve of biological soft tissues:

$$\psi =\frac{1}{4}\sum _{k=1}^{n}{\mu}_{k}\left[\frac{1}{{\alpha}_{k}}({{\textstyle \text{e}}}^{{({I}_{[k]}-1)}^{{\alpha}_{k}}}-1)+\frac{1}{{\beta}_{k}}({J}_{[k]}^{{\beta}_{k}}-1)+\frac{1}{{\chi}_{k}}({{I}_{3}}^{-{\chi}_{k}}-1)\right].$$

For both strain energy functions *ψ*, the polyconvexity condition is fulfilled if the material coefficients ${\mu}_{k},{\alpha}_{k}.{\beta}_{k}$ and ${\chi}_{k}$ satisfy the following inequalities:

$${\mu}_{k}\ge 0,\phantom{\rule{1em}{0ex}}{\alpha}_{k}\ge 1,\phantom{\rule{1em}{0ex}}{\beta}_{k}\ge 1,\phantom{\rule{1em}{0ex}}{\textstyle \text{and}}\phantom{\rule{1em}{0ex}}{\chi}_{k}\ge -\frac{1}{2}\phantom{\rule{1em}{0ex}}\{k=1,2,\dots ,n\}.$$

4.27

Although the function developed by Limbert and Middleton was able to fit very well the data from Lanir & Fung [54], it featured 12 parameters with no direct physical interpretation which, depending on the intended application, might be a limiting factor.

Gasser–Ogden–Holzapfel (GOH)'s constitutive formulation [97] was designed to capture the orthotropic hyperelastic behaviour of arterial tissues while accounting for fibre splay around two main directions corresponding to each collagen fibre families, characteristic of arterial microstructure. It is an extension of the original model developed by Gasser and co-workers [108] and is based on the introduction of a new structural tensor **H*** _{i}* accounting for fibre dispersion of fibre family

$${\mathbf{\text{H}}}_{i}={\kappa}_{i}\mathbf{\text{1}}+(1-3{\kappa}_{i}){\mathit{n}}_{0}^{i}\otimes {\mathit{n}}_{0}^{i}={\kappa}_{i}\mathbf{\text{1}}+(1-3{\kappa}_{i}){\mathfrak{B}}_{0}^{i},$$

4.28

where *κ _{i}* is a measure of fibre dispersion. The deceptively simple form for

$$\psi =\frac{\mu}{2}({I}_{1}-3)+\sum _{i=1}^{2}\frac{{k}_{i1}}{{k}_{ki2}}(\mathrm{exp}\{{k}_{ki2}{[{\mathbf{\text{H}}}_{i}\hspace{0.17em}:\hspace{0.17em}\mathbf{\text{C}}-1]}^{2}\}-1),$$

4.29

where *μ*, *k _{i}*

The GOH model was used in the context of skin mechanics by Ní Annaidh *et al.* [34] who conducted a series of physical uniaxial tensile tests to failure using digital image correlation (DIC) to measure stress–strain characteristics and mean collagen fibre distribution [41]. The physical characterization was conducted on cadaveric human skin specimens at various body locations and along several orientations (defined with respect to Langer lines).

GOH's model was fitted to the experimental tensile stress–strain curves and implemented into the finite-element environment of Abaqus (Simulia, Dassault Systèmes, Providence, RI, USA). The physical experiments were replicated by means of finite-element analysis and demonstrated very good performance of the numerical model.

Combining 3D DIC and bulge tests, Tonge *et al.* [50] determined in-plane anisotropic mechanical properties of post-mortem human skin using cyclic full-field measurements. Two main directions of anisotropy were considered and a series of full skin samples (located on the back torso) obtained from six male and female donors (43, 44, 59, 61, 62 and 83 years old) were used. The effect of preconditioning and humidity of the sample on the stress–strain response was investigated and was shown to be negligible. Age of the donors had a significant effect on the stiffness and directional properties of the skin. Specimens for older donors exhibited a stiffer and more isotropic response compared to those of younger donors. The authors also found that the bulge test method was limited by its inability to accurately determine stress and material parameters due to significant bending effects. In a companion paper, Tonge *et al.* [109] analysed the results of their bulge tests [50] using an analytical method based on thin shell theory which considered the effects of bending stiffness of the skin. The method accounts for through-the-thickness linear strain gradients. These authors fitted the shell version of GOH's model featuring a fully integrated fibre distribution to their experimental data. Two cases were considered for the GOH model—2D and 3D fibre distributions—while the fully integrated fibre model was restricted to 2D planar fibre orientation. It was found that both the 2D and 3D GOH model were unable to capture the anisotropic mechanics of skin from bulge tests unlike the 2D fully integrated fibre model which was shown to capture it very well. Tonge *et al.* [109] attributed the differences between their results and those of Ní Annaidh *et al.* [34,41] mainly to the younger age of their donors, lower strain level considered in their tests and their assumption about fibre orientation. Tonge *et al.* [109] considered only one fibre family aligned with the principal stretch direction while Ní Annaidh *et al.* [34,41] assumed that skin was made of two fibre families symmetric about the loading axis. In the context of surgical simulation procedures, the GOH model was used by Buganza-Tepole *et al.* [63] to model the mechanics of skin. In that study, computed stress profiles in skin flap were used as a surrogate measure of likelihood of tissue necrosis, following reconstructive surgery.

To address the issue of computationally expensive integration methods required for microstructural models of soft tissues featuring statistical distributions of material or structural characteristics (e.g. integration of fibre splay angular distribution over the unit sphere, distribution of engaged fibres) Flynn *et al.* [110] proposed a new model termed ‘discrete fibre icosahedral structural model’. The model relies on the use of a *discrete* rather than *continuous* fibre orientation distribution kernel which allows closed-form solutions for strain energy and stress. Six discrete directions are considered: they are oriented parallel to the lines connecting opposing vertices of a regular icosahedron. Six unit vectors ${\mathit{n}}_{0}^{i}\{i=1,\dots ,6\}$ corresponding to distinct fibre bundles are thus defined (electronic supplementary material, §S1.5). Flynn *et al.* [110] introduced a strain energy function per unit mass of the collection of six fibre bundles which features separate contributions for collagen and elastin fibres:

$$\psi \mathbf{\text{=}}\sum _{i=1}^{6}{w}_{i}[{\psi}_{\mathrm{c}}({\lambda}_{i})+{\psi}_{\mathrm{e}}({\lambda}_{i})]+U(J),$$

4.30

with

$${\psi}_{\mathrm{e}}({\lambda}_{i})=\frac{1}{{\rho}_{0}}\sum _{i=1}^{6}{w}_{i}\frac{{E}_{\mathrm{e}}}{4}\u27e8{\lambda}_{i}^{2}-1\u27e9.$$

4.31

where *w _{i}* are positive weights, associated with each of the six structural tensors, that balance the respective structural contributions of each fibre direction and

$${\int}_{1}^{\mathrm{\infty}}D(x){\textstyle \text{d}}x.$$

4.32

This means that the fraction of fibres that are fully taut (i.e. straight) at a stretch *λ* is given by the following integral, where *D* is an undulation probability distribution:

$${\int}_{1}^{\lambda}D(x){\textstyle \text{d}}x.$$

4.33

The stiffness of any of the six collagen bundles is:

$$\frac{{\textstyle \text{d}}}{{\textstyle \text{d}}\lambda}\left(\frac{\lambda {\rho}_{0}}{J}\frac{{\textstyle \text{d}}{\psi}_{\mathrm{c}}(\lambda )}{{\textstyle \text{d}}\lambda}\right)={E}_{\mathrm{c}}{\int}_{1}^{\lambda}D(x){\textstyle \text{d}}x.$$

4.34

Assuming that all the collagen fibres are slack and are stress-free for stretch lower or equal to unity, and further assuming that deformations are isochoric:

$$\frac{{\textstyle \text{d}}{\psi}_{\mathrm{c}}(\lambda )}{{\textstyle \text{d}}\lambda}=\frac{1}{\lambda {\rho}_{0}}{\int}_{1}^{\lambda}\left[{E}_{\mathrm{c}}{\int}_{1}^{\lambda}D(x){\textstyle \text{d}}x\right]{\textstyle \text{d}}y.$$

4.35

Flynn *et al.* [110] considered two undulation probability distributions: step and triangle distributions. They offer the advantage of being simple enough so that a closed-form analytical expression can be obtained for *ψ*_{c}. The example of a using unit step distribution is provided in electronic supplementary material, §S1.5. This strain energy function was used by Flynn *et al.* [110] to model the biaxial tensile response of rabbit skin [54] (8.7% error) and uniaxial response of porcine skin [111] (7.6% error). The constitutive model was later extended by Flynn & Rubin [112] to address shortcomings linked to the relation between fibre weight and anisotropic response, namely, the fact that for equal weights *w _{i}* the mechanical response is not necessarily isotropic. These authors introduced a generalized strain invariant

$$\gamma \mathbf{\text{= (}}\mathbf{\text{C}}+{\mathbf{\text{C}}}^{-1}\mathbf{\text{)}}\cdot \sum _{i=1}^{6}{w}_{i}{\mathfrak{B}}_{0}^{i}-2,$$

4.36

where *w _{i}* are positive weights, associated with each of the six structural tensors ${\mathfrak{B}}_{0}^{i}={\mathit{n}}_{0}^{i}\otimes {\mathit{n}}_{0}^{i}$, that balance the respective structural contributions of each fibre direction.

Unlike for the earliest model of Flynn *et al.* [110], here, equal weights can only lead to isotropy of mechanical properties. Finally, Flynn & Rubin [112] proposed the following strain energy density:

$$\psi =\frac{{\sigma}_{0}}{2{\rho}_{0}}\left[\gamma +\sum _{m=1}^{M}\frac{{\gamma}_{\mathrm{m}}}{m}\left(\frac{\gamma}{{\gamma}_{\mathrm{m}}}\right)\right],$$

4.37

where *ρ*_{0} is the mass density in the reference configuration, *σ*_{0} is a stress-like material parameter, *γ*_{m} are dimensionless positive material parameters and *M* is the order of the polynomial expansion. Using the same experimental data as in their 2011 paper [54,111], Flynn & Rubin [112] identified constitutive parameters of their new strain energy function (equation (4.37)). The fit was not as good as in their previous study [110] as the fitting errors for rabbit and pig skin data were, respectively, 12% and 17%. However, in that case, the weights were a pure measure of anisotropy.

Most of the microstructurally based constitutive models of soft tissues assume an additive decomposition of the strain energy function into decoupled matrix and fibre elastic potentials. This means that fibre–fibre and matrix–fibre interactions are not *explicitly* captured in the constitutive formulation. The network models based on BAG's formulation captures only *implicitly* and *globally* these interactions. To address this first shortcoming, Limbert [113] developed a novel invariant-based multi-scale constitutive framework to characterize the transversely isotropic and orthotropic elastic responses of biological soft tissues. The constitutive equations were particularized to model skin. The model was not only capable to accurately reproduce the experimental multi-axial behaviour of rabbit skin, as in [114] but could also *a posteriori* predict stiffness values of individual tropocollagen molecules in agreement with physical and molecular-dynamics-based computational experiments [113]. Another key motivational aspect of the formulation is that the constitutive parameters can be directly extracted from physical measurements by segregating the orthogonal deformation modes of its constituents. Another desirable feature of the constitutive equations is that the response is based on physical geometrical/structural parameters that can be measured experimentally or determined *ab initio* from molecular dynamic simulations.

Limbert's formulation is based on the constitutive framework of Lu & Zhang [115] for transversely isotropic materials which make use of four invariants that characterize decoupled deformations modes solely related to purely volumetric (*J*), deviatoric stretch in the fibre direction $(\overline{\lambda})$, cross-fibre shear $({\alpha}_{1}^{i})$ and fibre-to-fibre/matrix-to-fibre shear $({\alpha}_{2}^{i})$ stress responses. Orthotropic symmetry is accounted for by introducing a second family of fibres. The index *i*=1,2 identifies each fibre family:

$$J=\sqrt{{I}_{3}};\phantom{\rule{1em}{0ex}}{\overline{\lambda}}_{i}={{I}_{3}}^{-1/6}\sqrt{{I}_{4}^{i}};\phantom{\rule{1em}{0ex}}{\alpha}_{1}^{i}=\frac{{I}_{1}{I}_{4}^{i}-{I}_{5}^{i}}{\sqrt{{I}_{3}{I}_{4}^{i}}}\phantom{\rule{1em}{0ex}}{\textstyle \text{and}}\phantom{\rule{1em}{0ex}}{\alpha}_{2}^{i}=\frac{{I}_{5}^{i}}{{({I}_{4}^{i})}^{2}}.$$

4.38

Limbert proposed the following strain energy function:

$$\psi ={\psi}^{\mathrm{v}}(J)+\sum _{i=1}^{2}[{\psi}_{i}^{\overline{\lambda}}({\overline{\lambda}}_{i})+{\hat{\psi}}_{i}^{1}({\alpha}_{1}^{i})+{\stackrel{~}{\psi}}_{i}^{2}({\alpha}_{2}^{i},{\overline{\lambda}}_{i})],$$

4.39

${\psi}^{\mathrm{v}}$, ${\psi}_{i}^{\overline{\lambda}}$, ${\hat{\psi}}_{i}^{1}$ and ${\stackrel{~}{\psi}}_{i}^{2}$ are, respectively, the volumetric, deviatoric fibre, cross-fibre shear and fibre-to-fibre/fibre-to-matrix shear energies. The functional forms of the energies and the constitutive parameters are detailed below.

$${\hat{\psi}}_{i}^{1}({\alpha}_{1}^{i})={\textstyle \frac{1}{2}}{\mu}_{1}^{i}({\alpha}_{1}^{i}-2)\phantom{\rule{1em}{0ex}}{\textstyle \text{and}}\phantom{\rule{1em}{0ex}}{\stackrel{~}{\psi}}_{i}^{2}({\alpha}_{2}^{i},{\overline{\lambda}}_{i})={\textstyle \frac{1}{2}}{\mu}_{2}^{i}{({\alpha}_{2}^{i}-1)}^{2}\underset{\mathrm{Sigmoid}\hspace{0.17em}\mathrm{coupling}\hspace{0.17em}\mathrm{function}}{\underset{\u23df}{\frac{1}{1+{a}_{i}\hspace{0.17em}{{\textstyle \text{e}}}^{-{b}_{i}({\overline{\lambda}}_{i}-{\overline{\lambda}}_{i}^{c})}}}}.$$

4.40

A notable feature of this constitutive approach is that collagen fibres and matrix are allowed to interact via explicit decoupled shear interactions while collagen fibres behave like a worm-like chain model [104,116].

$${\psi}_{i}^{\overline{\lambda}}({\overline{\lambda}}_{\hspace{0.17em}i})=\{\begin{array}{ll}{\mathrm{\aleph}}_{i}{\mu}_{0}^{i}\left({\overline{\lambda}}_{i}^{2}+\frac{2}{{\overline{\lambda}}_{i}}-3\right)+{\xi}_{0}^{i}{\mathrm{\aleph}}_{i}\mathrm{ln}({\overline{\lambda}}_{i}^{{r}_{0i}^{2}})& {\textstyle \text{if}}{\overline{\lambda}}_{i}\le 1,\\ {\mathrm{\aleph}}_{i}\mathfrak{K}\theta \frac{{L}_{i}}{4{L}_{\mathrm{p}}^{i}}\left(2\frac{{\overline{\lambda}}_{i}^{2}{r}_{0\beta}^{2}}{{L}_{i}^{2}}+{\left(1-\frac{{\overline{\lambda}}_{i}{r}_{0i}}{{L}_{i}}\right)}^{-1}-\frac{{\overline{\lambda}}_{i}{r}_{0i}}{{L}_{i}}\right)+{\xi}_{0}^{i}{\mathrm{\aleph}}_{i}\mathrm{ln}({\overline{\lambda}}_{i}^{{r}_{0\beta}^{2}})& {\textstyle \text{if}}{\overline{\lambda}}_{i}1.\end{array}$$

4.41

The constitutive model resulted in a set of twenty-three constitutive parameters $\mathbf{\text{p}}={\mathbf{\text{p}}}^{1}\cup {\mathbf{\text{p}}}^{2}={\{\mathfrak{K},\theta ,\kappa ,{\mu}_{0}^{i},{\mu}_{1}^{i},{\mu}_{2}^{i},{\mathrm{\aleph}}_{i},{L}_{i},{L}_{p}^{i},{r}_{i},{a}_{i},{b}_{i},{\overline{\lambda}}_{i}^{\mathrm{c}}\}}_{i=1,2}$. Moreover, all the constitutive parameters of this multi-scale formulation have a direct physical interpretation. Limbert [113] determined the parameter set **p**^{2} using a numerical global optimization algorithm while the parameter set **p**^{1} was assumed *a priori* based on existing data [104] and data obtained via visual inspection of the biaxial stress–strain curves [114]. The physical meaning of these parameters is provided in electronic supplementary material, table S1. Although the model slightly under-predicts the response of rabbit skin along the head–tail direction at low stretch (less than 1.35), an excellent fit was obtained (electronic supplementary material, figure S9). Limbert implemented the nonlinear constitutive model as tri-linear hexahedral finite element using an enhanced strain formulation [117] which has been proven to be superior to a standard displacement-based formulation, particularly for shear-dominated problems and nearly incompressible materials. Moreover, analytically exact direct sensitivity analyses subroutines were also implemented to assess the sensitivity of the shear response of the model to its constitutive parameters during a simulated indentation test [113]. From the persistence lengths of tropocollagen molecules determined from the optimization procedure (22 and 65nm for fibre families 1 and 2), Limbert [113] calculated equivalent Young's moduli of, respectively, 293 and 865MPa. These values lie within one order of magnitude less than what has been determined experimentally and obtained through computational modelling studies [118]. However, using the same equation used by Limbert to calculate Young's modulus, Sun *et al.* [119] estimated Young's modulus of collagen molecules to range between 350MPa and 12GPa.

Groves *et al.* [120] used Weiss's model [121] to model the anisotropic mechanics of skin by identifying constitutive parameters from a series of uniaxial tensile tests. The original model proposed by Weiss *et al.* [121] was formulated to capture the transversely isotropic hyperelastic response of ligaments of the knee joint but without including fibre dispersion which is not significant in these soft tissues.

The strain energy function used by Groves *et al.* [120] was defined as the sum of a Veronda–Westmann (VW) potential *ψ*^{VW} [56] and a piece-wise anisotropic fibre function ${\psi}_{i}^{\mathrm{fibre}}$ [121], to model the isotropic and anisotropic responses, respectively:

$$\psi ({I}_{1},{I}_{2},{I}_{4}^{1},{I}_{4}^{2},{I}_{4}^{3})={\psi}^{\mathrm{VW}}({I}_{1},{I}_{2})+\sum _{i=1}^{3}{\psi}_{i}^{\mathrm{fibre}}({I}_{4}^{i}={\lambda}^{i}),$$

4.42

where

$${\psi}^{\mathrm{VW}}={c}_{1}\{\mathrm{exp}[{c}_{2}({I}_{1}-3)]-1\}-\frac{{c}_{1}{c}_{2}}{2}({I}_{2}-3).$$

4.43

Here *c*_{1} is a shear-like modulus while *c*_{2} is a dimensionless parameter that scales the response connected to the second invariant and

$$\frac{\mathrm{\partial}{\psi}_{i}^{\mathrm{fibre}}({\lambda}^{i})}{\mathrm{\partial}{\lambda}^{i}}=\frac{1}{{\lambda}^{i}}\{\begin{array}{ll}0& {\textstyle \text{if}}{\lambda}^{i}\le 1,\\ {c}_{3}\{\mathrm{exp}[{c}_{4}({\lambda}^{i}-1)]-1\}& {\textstyle \text{if}}1{\lambda}^{i}{\overline{\lambda}}^{i},\\ {c}_{5}{\lambda}^{i}+{c}_{6}& {\textstyle \text{if}}{\lambda}^{i}\ge {\overline{\lambda}}^{i}.\end{array}$$

4.44

Here ${\overline{\lambda}}^{i}=\sqrt{{\overline{I}}_{4}^{i}}$ represents the stretch at which the (collagen) fibres are assumed to be fully taut. *c*_{3} is an elastic modulus-like parameter scaling the exponential response, *c*_{4} controls the rate of un-crimping of the fibres, *c*_{5} is the elastic modulus of the taut fibre while *c*_{6} is a correcting factor to ensure continuity of the stiffness response at ${\lambda}^{i}={\overline{\lambda}}^{i}$.

In their experimental procedure, Groves *et al.* [120] used eight human discoid skin samples from two female donors (aged 56 and 68 years) following mastectomy and 14 post-mortem murine skin samples obtained from eight mice (aged 18–24 months). For each sample, tensile tests were simultaneously conducted along three axes (0, 45 and 90°). Using an inverse analysis based on finite-element simulations of the testing procedures constitutive parameters were determined. For each sample, three sets of four parameters (*c*_{3}, *c*_{4}, $\sqrt{{\overline{I}}_{4}^{i}}$ and an additional parameter characterizing the deviation from the assumed fibre orientation) for each fibre energy function were obtained in addition to the two parameters of the VW function. Although Groves *et al.* [120] demonstrated a good fit between their three-fibre family model and their experimental data, they also acknowledged the limitations of using a Simplex optimization algorithm for their inverse analysis which could only capture *local* minima of their cost function.

The skin exhibits a wide range of viscoelastic effects including stress relaxation, creep, rate dependency of stress and hysteresis [8]. A large number of researchers have characterized the viscoelastic properties of skin under various testing conditions [55,122–129]. Observed macroscopic viscoelastic effects can be attributed to two main mechanisms that operate in concert. Viscoelasticity of the skin can originate from the intrinsic viscoelastic characteristics of its nano- and microstructural constitutents (e.g. proteoglycans, elastic fibres) and also from the time-dependent rearrangement of its microstructure under macroscopic loads which takes a finite amount of time to occur (flow of interstitial fluid across the porous structures, progressive sliding of collagen fibrils past each other). Proteoglycan macromolecules such as decorin bound on the collagen fibrils and play an important role in these interactions which are mediated by the side chains of glycosaminoglycans. Because the forces acting between fibrils are of a non-covalent nature, these links can break and reconnect [55,130,131]. One should also note that there are other complex physics phenomena at play. For example, the thermal motion of ions, including Na^{+} and K^{+}, displacing towards the negatively charged ends of glycosaminoglycans of fibril-associated proteoglycans, induce an attractive force between two collagen fibrils.

To date, there are only very few nonlinear viscoelastic anisotropic constitutive models of skin ([55,132–135] and references therein). This is partly due to the considerable challenges of experimentally characterizing the behaviour of skin, particularly *in vivo*, and the success of simpler theories such as quasi-linear viscoelasticity (QLV) [8,136].

The literature on general viscoelastic constitutive models is rich and, here, only a brief account of the most common approaches to model biological soft tissues (following the classification of Ehret [137]) is reported. More recent and state-of-the-art nonlinear viscoelastic models of skin are presented.

The application of QLV theory to soft tissue mechanics has been particularly popularized by the work of Fung [8,136]. The idea behind QLV is to assume that the time-dependent stress response *σ*(*t*) following uniaxial loading can be expressed as a convolution integral of the form:

$$\sigma (t)={\int}_{-\mathrm{\infty}}^{t}G(t-\tau )\frac{\mathrm{\partial}{\sigma}_{\mathrm{e}}(\tau )}{\mathrm{\partial}t}\hspace{0.17em}{\textstyle \text{d}}\tau ,$$

4.45

where *σ*_{e} is the instantaneous elastic stress and *G* the reduced relaxation function. This function controls how the current stress response is modulated by past loading history.

The notion of quasi-linearity stems from the linear relation of the integrand terms and the analogy with rate equations of linear viscoelasticity. The uniaxial relationship can be extended to a full 3D representation by introducing second-order and fourth-order tensors describing, respectively, the stress and reduced relaxation function [8]. This represents a flexible framework where anisotropic properties, and decoupling between deviatoric and volumetric responses can be captured by appropriate constitutive equations.

Bischoff [138] applied QLV at the collagen fibre level to model porcine skin by developing an anisotropic microstructurally motivated constitutive model which also includes fibre dispersion. A notable feature of the model is the ability to capture a fibre-level viscoelastic orthotropic response with only seven parameters. However, the authors recognized the need for additional experiments to fully characterize the mechanical response as numerical identification procedures are plagued by non-unicity of constitutive parameters. Additional characterization tests would likely lead to convergence towards a single set of optimal parameters.

By design, these types of model—also called viscoelastic models of the *differential* type—capture strain-rate sensitivity and *short-term* viscous effects following application of load [139]. They are based on *differential* rather than *integral* equations like for QLV, and are therefore inappropriate to capture long-term memory viscoelastic effects such as relaxation. Long-term memory effects encompass the whole deformation history of the material. The so-called *principle of fading memory* [140] states that deformations that occurred in the recent time history have greater influence on the actual stresses than those which occurred in a more distant past history.

Explicitly rate-dependent models of soft tissues are particularly well suited to high strain-rate situations such as those occurring during vehicle accidents, sport activities, impact and blast [139,141–144]. These models generally postulate the existence of a viscous potential *ψ*^{v} from which the viscous dissipative effects (denoted by scalar $\mathcal{D}$) arise by differentiation with respect to the rate of the Cauchy–Green deformation tensor $\dot{\mathbf{\text{C}}}=\mathrm{\partial}\mathbf{\text{C}}/\mathrm{\partial}t$:

$$\mathcal{D}=\frac{\mathrm{\partial}{\psi}^{\mathrm{v}}}{\mathrm{\partial}\dot{\mathbf{\text{C}}}}:\dot{\mathbf{\text{C}}}.$$

4.46

The total second Piola–Kirchhoff stress tensor which includes purely elastic effects through a potential *ψ*^{e}, and purely viscous effects through *ψ*^{v}, is obtained by differentiation of the total potential *ψ* assuming that $\dot{\mathbf{\text{C}}}$ is a parameter or internal variable:

$$\mathbf{\text{S}}=2\left(\frac{\mathrm{\partial}{\psi}^{\mathrm{e}}}{\mathrm{\partial}\mathbf{\text{C}}}+\frac{\mathrm{\partial}{\psi}^{\mathrm{v}}}{\mathrm{\partial}\dot{\mathbf{\text{C}}}}\right).$$

4.47

The framework proposed by Pioletti *et al.* [139] for modelling the rate-dependent isotropic viscohyperelastic behaviour of ligaments and tendons was later extended by Limbert & Middleton [142] to transverse isotropy using tensor formalism. This model was subsequently combined with damage equations to model the failure of human skin in response to puncturing biting loads [141].

Unidimensional (small-strain) rheological models based on combinations of spring and dashpot elements are efficient means to conceptualize particular viscoelastic behaviours [87]. The arrangement of a dashpot and a spring in series constitutes a Maxwell element with elastic modulus *E* and viscosity *η* and the total strain *ε* in the element can be decomposed into elastic and inelastic strains as $\epsilon ={\epsilon}_{\mathrm{e}}+{\epsilon}_{\mathrm{i}}$ and the total strain rate is

$$\frac{\mathrm{\partial}\epsilon}{\mathrm{\partial}t}=\frac{1}{E}\frac{\sigma}{\eta}+\frac{\mathrm{\partial}\sigma}{\mathrm{\partial}t}.$$

4.48

A very simple approach to account for the viscoelastic properties of skin would be to assume that it is a generalized standard viscoelastic solid made of *n* Maxwell elements so that it can be described by a strain energy function of the form $\psi =\mathcal{P}(t){\psi}^{\mathrm{\infty}}$, where *ψ*^{∞} is the time-independent or instantaneous strain energy which is convolved with a time-dependent kernel represented by a *n*-terms Prony series $\mathcal{P}$ defined through characteristic times and moduli. *τ _{i}* and

$$\mathcal{P}(t)=1+\sum _{i=1}^{n}{c}_{i}\left[1-\mathrm{exp}\left(-\frac{t}{{\tau}_{i}}\right)\right].$$

4.49

Borrowing the concept of linear strain decomposition and applying it to the 3D finite strain regime, it is possible to establish a multiplicative decomposition of the deformation gradient into an elastic and inelastic parts as $\mathbf{\text{F}}={\mathbf{\text{F}}}^{\mathrm{e}}\cdot {\mathbf{\text{F}}}^{\mathrm{i}}$ [145,146]. From this definition one can define the fictitious *elastic* and *inelastic* right Cauchy–Green deformation tensors as:

$${\mathbf{\text{C}}}^{\mathrm{e}}={\mathbf{\text{F}}}^{\mathrm{e}\hspace{0.17em}{\textstyle \text{T}}}\cdot {\mathbf{\text{F}}}^{\mathrm{e}}\phantom{\rule{1em}{0ex}}{\textstyle \text{and}}\phantom{\rule{1em}{0ex}}{\mathbf{\text{C}}}^{\mathrm{i}}={\mathbf{\text{F}}}^{\mathrm{i}\hspace{0.17em}\mathrm{T}}\cdot {\mathbf{\text{F}}}^{\mathrm{i}}.$$

4.50

One can then postulate the existence of a free energy $\psi ({\mathbf{\text{C}}}^{\mathrm{e}},{\mathbf{\text{C}}}^{\mathrm{i}})$ decomposed into *equilibrium* and *non-equilibrium terms* where **C**^{i} is treated as an internal variable [145]:

$$\psi (\mathbf{\text{C}})=\psi ({\mathbf{\text{C}}}^{\mathrm{e}},{\mathbf{\text{C}}}^{\mathrm{i}})={\psi}^{\mathrm{equilibrium}}({\mathbf{\text{C}}}^{\mathrm{e}})+{\psi}^{n.\mathrm{equilibrium}}({\mathbf{\text{C}}}^{\mathrm{i}}).$$

4.51

Naturally, one can generalize this concept to *n* Maxwell elements so that one introduces *n* internal (tensor) variables ${\mathbf{\text{C}}}_{k}^{\mathrm{i}}$, arguments of *n* non-equilibrium potentials ${\psi}_{k}^{n.\mathrm{equilibrium}}$:

$$\psi (\mathbf{\text{C}},{\mathbf{\text{C}}}_{k}^{\mathrm{i}})={\psi}^{\mathrm{equilibrium}}(\mathbf{\text{C}})+\sum _{k=1}^{n}{\psi}_{k}^{n.\mathrm{equilibrium}}({\mathbf{\text{C}}}_{k}^{\mathrm{i}})$$

4.52

The second Piola–Kirchhoff stress is then obtained as the sum of equilibrium and non-equilibrium terms as:

$$\mathbf{\text{S}}={\mathbf{\text{S}}}^{\mathrm{equilibrium}}+\sum _{k=1}^{n}{\mathbf{\text{S}}}_{k}^{n.\mathrm{equilibrium}}=2\left[\frac{\mathrm{\partial}{\psi}^{\mathrm{equilibrium}}({\mathbf{\text{C}}}^{\mathrm{e}})}{\mathrm{\partial}{\mathbf{\text{C}}}^{\mathrm{e}}}+\sum _{k=1}^{n}\frac{{\psi}_{k}^{n.\mathrm{equilibrium}}({\mathbf{\text{C}}}_{k}^{i})}{\mathrm{\partial}{\mathbf{\text{C}}}_{k}^{i}}\right].$$

4.53

A notable feature of the present formulation is that non-equilibrium or viscous deformations are not restricted to the small strain regime unlike the vast majority of viscoelastic models of biological soft tissues found in the literature.

Bischoff *et al.* [132] combined this type of formulation with their previous orthotropic eight-chain model [99] to model the viscoelastic behaviour of soft tissues. The model was fitted to the experimental data on rabbit skin [54]. Vassoler *et al.* [147] recently proposed a variational framework making use of the multiplicative decomposition of the deformation gradient into elastic and inelastic parts to represent the mechanics of fibre-reinforced biological composites and this would be appropriate for skin.

Limbert (private communication) extended his decoupled invariant orthotropic hyperelastic model to finite strain viscoelasticity following a similar approach to that described in this section. The particular approach for anisotropic viscoelastic effects follows the tensor formalism of Nguyen *et al*. [148] and Nedjar [149] who introduced viscosity tensors. An essential feature of the model is that the matrix and the fibre phase are treated separately, each featuring their own deformation gradient. The nonlinear creep and/or relaxation response is based on the multiplicative viscoelastic split of the deformation gradient combined with the assumption of the existence of viscoelastic potentials for each phase. The deformation gradient and its multiplicative decomposition apply to all the continua (matrix and fibre phase) linking them implicitly. Separate flow rules are specified for the matrix and fibre families.

The flow rules of the fibre families are combined to provide an anisotropic flow rule of the fibre phase. Details of the constitutive formulation and finite-element implementation of the present constitutive model for the dermis can be found in Nguyen *et al*.'s paper [148]. One should point out that, besides the *matrix* phase, the multiplicative decomposition of the deformation gradient applies to the *fibre phase*, not to the individual families of fibres (two, in our case).

$$\mathbf{\text{F}}={\mathbf{\text{F}}}_{\mathrm{matrix}\hspace{0.17em}(k)}^{\mathrm{e}}{\mathbf{\text{F}}}_{\mathrm{matrix}\hspace{0.17em}(k)}^{\mathrm{v}}={\mathbf{\text{F}}}_{\mathrm{fibre}\hspace{0.17em}(p)}^{\mathrm{e}}{\mathbf{\text{F}}}_{\mathrm{fibre}\hspace{0.17em}(p)}^{\mathrm{v}}.$$

4.54

The elastic right Cauchy–Green deformation tensors associated with the matrix and fibre phases are defined as:

$${\mathbf{\text{C}}}_{\mathrm{matrix}}^{\mathrm{e}}={({\mathbf{\text{F}}}_{\mathrm{matrix}}^{\mathrm{v}})}^{-\mathrm{T}}({\mathbf{\text{F}}}^{\mathrm{T}}\mathbf{\text{F}}){({\mathbf{\text{F}}}_{\mathrm{matrix}}^{\mathrm{v}})}^{-1}$$

4.55

and

$${\mathbf{\text{C}}}_{\mathrm{fibre}}^{\mathrm{e}}={({\mathbf{\text{F}}}_{\mathrm{fibre}}^{\mathrm{v}})}^{-\mathrm{T}}({\mathbf{\text{F}}}^{\mathrm{T}}\mathbf{\text{F}}){({\mathbf{\text{F}}}_{\mathrm{fibre}}^{\mathrm{v}})}^{-1}.$$

4.56

The general formulation of the free energy density can be expressed as the sum of equilibrium and non-equilibrium energies depending on the total and viscous deformation tensors and the fibre orientation vectors (defined by *n*_{0} and *m*_{0} if one considers two families of collagen fibres):

$$\psi (\mathbf{\text{C}},{\mathbf{\text{C}}}_{\mathrm{matrix}}^{\mathrm{v}},{\mathbf{\text{C}}}_{\mathrm{fibre}}^{\mathrm{v}},{\mathit{n}}_{0},{\mathit{m}}_{0})=\underset{\mathrm{equilibrium}\hspace{0.17em}\mathrm{part}}{\underset{\u23b5}{{\psi}^{\mathrm{\infty}}(\mathbf{\text{C}},{\mathit{n}}_{0},{\mathit{m}}_{0})}}+\underset{\mathrm{non}{\textstyle \text{-}}\mathrm{equilibrium}\hspace{0.17em}\mathrm{part}}{\underset{\u23b5}{{\psi}^{\mathrm{v}}({\mathbf{\text{C}}}_{\mathrm{matrix}}^{\mathrm{e}},{\mathbf{\text{C}}}_{\mathrm{fibre}}^{\mathrm{e}},{\mathit{n}}_{0},{\mathit{m}}_{0})}},$$

4.57

where

$${\psi}^{\mathrm{\infty}}(\mathbf{\text{C}},{\mathit{n}}_{0},{\mathit{m}}_{0})={\psi}_{M}^{\mathrm{\infty}}({I}_{1},{I}_{2},{I}_{3})+{\psi}_{F}^{\mathrm{\infty}}({I}_{4},{I}_{5},{I}_{6},{I}_{7})$$

4.58

and

$${\psi}^{\mathrm{v}}({\mathbf{\text{C}}}_{\mathrm{matrix}}^{\mathrm{e}},{\mathbf{\text{C}}}_{\mathrm{fibre}}^{\mathrm{e}},{\mathit{n}}_{0},{\mathit{m}}_{0})={\mathcal{W}}_{\mathrm{M}}^{\mathrm{v}}({I}_{1(M)}^{\mathrm{e}},{I}_{2(M)}^{\mathrm{e}},{I}_{3(M)}^{\mathrm{e}})+{\mathcal{W}}_{\mathrm{F}}^{\mathrm{v}}({I}_{4(\mathrm{F})}^{\mathrm{e}},{I}_{5(\mathrm{F})}^{\mathrm{e}},{I}_{6(\mathrm{F})}^{\mathrm{e}},{I}_{7(\mathrm{F})}^{\mathrm{e}}),$$

4.59

with

$$\begin{array}{rl}{I}_{4}& \hspace{0.17em}=\mathbf{\text{C}}:({\mathit{n}}_{0}\otimes {\mathit{n}}_{0});\phantom{\rule{1em}{0ex}}{I}_{5}={\mathbf{\text{C}}}^{2}:({\mathit{n}}_{0}\otimes {\mathit{n}}_{0});\phantom{\rule{1em}{0ex}}{I}_{6}=\mathbf{\text{C}}:({\mathit{m}}_{0}\otimes {\mathit{m}}_{0});\phantom{\rule{1em}{0ex}}{I}_{7}={\mathbf{\text{C}}}^{2}:({\mathit{m}}_{0}\otimes {\mathit{m}}_{0}),\end{array}$$

4.60

$$\begin{array}{rl}{I}_{1(k)}^{\mathrm{e}}& \hspace{0.17em}={\mathbf{\text{C}}}_{(k)}^{\mathrm{e}}:\mathbf{\text{I}};\phantom{\rule{1em}{0ex}}{I}_{2(k)}^{\mathrm{e}}={\textstyle \frac{1}{2}}({{\mathbf{\text{C}}}_{(k)}^{\mathrm{e}}}^{2}:\mathbf{\text{I}}-{I}_{1(k)}^{\mathrm{e}});\phantom{\rule{1em}{0ex}}{I}_{3(k)}^{\mathrm{e}}={\textstyle \text{determinant}}({\mathbf{\text{C}}}_{(k)}^{\mathrm{e}})\end{array}$$

4.61

$$\begin{array}{rl}{\textstyle \text{and}}\phantom{\rule{2em}{0ex}}{I}_{4(k)}^{\mathrm{e}}& \hspace{0.17em}={\mathbf{\text{C}}}_{(k)}^{\mathrm{e}}:({\mathit{n}}_{0}\otimes {\mathit{n}}_{0});\phantom{\rule{1em}{0ex}}{I}_{5(k)}^{\mathrm{e}}={{\mathbf{\text{C}}}_{(k)}^{\mathrm{e}}}^{2}:({\mathit{n}}_{0}\otimes {\mathit{n}}_{0});\phantom{\rule{1em}{0ex}}{I}_{6(k)}^{\mathrm{e}}={\mathbf{\text{C}}}_{(k)}^{\mathrm{e}}:({\mathit{m}}_{0}\otimes {\mathit{m}}_{0});\phantom{\rule{1em}{0ex}}{I}_{7(k)}^{\mathrm{e}}={\mathbf{\text{C}}}_{(k)}^{\mathrm{e}},\end{array}$$

4.62

where

$$k={\textstyle \text{matrix or fibre}}$$

4.63

Finally, the specific form of the free energy using the decoupled invariants introduced in §4.2.6 is defined as follows:

$$\psi (\mathbf{\text{C}},{\mathbf{\text{C}}}_{\mathrm{matrix}}^{\mathrm{v}},{\mathbf{\text{C}}}_{\mathrm{fibre}}^{\mathrm{v}},{\mathit{n}}_{0},{\mathit{m}}_{0})=\underset{\mathrm{equilibrium}\hspace{0.17em}\mathrm{part}}{\underset{\u23b5}{{\psi}_{\mathrm{M}}^{\mathrm{\infty}}(J)+\sum _{\beta =1}^{2}[{\psi}_{\beta}^{\mathrm{\infty}}({\overline{\lambda}}_{\beta},{\alpha}_{1}^{\beta},{\alpha}_{2}^{\beta})]}}+\underset{\mathrm{non}{\textstyle \text{-}}\mathrm{equilibrium}\hspace{0.17em}\mathrm{part}}{\underset{\u23b5}{{\mathrm{\Phi}}_{\mathrm{M}}^{\mathrm{v}}({J}^{\mathrm{e}})+\sum _{\beta =1}^{2}[{\mathrm{\Phi}}_{\beta}^{\mathrm{v}}({\overline{\lambda}}_{\beta}^{\mathrm{e}},{\alpha}_{1e}^{\beta},{\alpha}_{2e}^{\beta})]}},$$

4.64

where invariants with superscript/subscript ‘e’ refer to those associated with the *elastic* right Cauchy–Green deformation tensor **C**^{e}. The equilibrium part of the free energy is identical to the total elastic energy defined in equation (4.39). The reduced dissipation inequality is

$$\underset{{\mathbf{\text{S}}}_{\mathrm{M}}^{\mathrm{v}}}{\underset{\u23df}{-2\frac{\mathrm{\partial}{\psi}_{\mathrm{M}}^{\mathrm{v}}}{\mathrm{\partial}{\mathbf{\text{C}}}_{\mathrm{matrix}}^{\mathrm{v}}}}}:\frac{1}{2}{\dot{\mathbf{\text{C}}}}_{\mathrm{matrix}}^{\mathrm{v}}\underset{{\mathbf{\text{S}}}_{\mathrm{F}}^{\mathrm{v}}}{\underset{\u23df}{-2\frac{\mathrm{\partial}{\psi}_{\mathrm{F}}^{\mathrm{v}}}{\mathrm{\partial}{\mathbf{\text{C}}}_{\mathrm{fibre}}^{\mathrm{v}}}}}:\frac{1}{2}{\dot{\mathbf{\text{C}}}}_{\mathrm{fibre}}^{\mathrm{v}}\ge 0,$$

4.65

where ${\mathbf{\text{S}}}_{\mathrm{M}}^{\mathrm{v}}$ and ${\mathbf{\text{S}}}_{\mathrm{F}}^{\mathrm{v}}$ are the stresses driving the viscous relaxation of the matrix and fibre phases. To satisfy the positive dissipation criterion for the matrix/fibre phase, the following evolution equations are proposed [145,148]:

$$\frac{1}{2}}{\dot{\mathbf{\text{C}}}}_{\mathrm{matrix}}^{\mathrm{v}}={\mathbb{V}}_{\mathrm{M}}^{-1}:{\mathbf{\text{S}}}_{\mathrm{M}}^{\mathrm{v}}\phantom{\rule{1em}{0ex}}{\textstyle \text{and}}\phantom{\rule{1em}{0ex}}{\textstyle \frac{1}{2}}{\dot{\mathbf{\text{C}}}}_{\mathrm{matrix}}^{\mathrm{v}}={\mathbb{V}}_{\mathrm{M}}^{-1}:{\mathbf{\text{S}}}_{\mathrm{M}}^{\mathrm{v}},$$

4.66

where ${\mathbb{V}}_{\mathrm{M}}^{-1}$ and ${\mathbb{V}}_{\mathrm{F}}^{-1}$ are the inverse of positive definite fourth-order isotropic/anisotropic viscosity tensors. These are defined in electronic supplementary material, §S1.7, together with particular non-equilibrium potentials for anisotropic viscoelasticity.

Flynn & Rubin [150] extended the original Flynn & Rubin's [112] and Flynn *et al.*'s [110] formulations to phenomenologically model dissipative effects in soft tissues. They followed the theoretical and numerical approach of Hollenstein *et al.* [151] which can capture both rate-independent and rate-dependent anelastic response including stress relaxation effects. The model was applied to stress relaxation data of rabbit skin [54] with various degrees of success depending on the magnitude of stretch [150].

Another popular approach for describing the viscous nonlinear elasticity of biological soft tissues is to consider a set of non-measurable (i.e. internal variable) strain-like internal variables **E*** _{k}* in the reference configuration, the free energy of a viscoelastic material can be defined as follows [87]:

$$\psi (\mathbf{\text{C}},{\mathbf{\text{E}}}_{k})={\psi}^{\mathrm{equilibrium}}(\mathbf{\text{C}})+\sum _{k=1}^{n}{\psi}_{k}^{n.\mathrm{equilibrium}}(\mathbf{\text{C}},{\mathbf{\text{E}}}_{k}),$$

4.67

where, in analogy with equation (4.52), the free energy (and stress) can be split into equilibrium and non-equilibrium contributions associated, respectively, with elastic and viscous deformation mechanisms:

$$\mathbf{\text{S}}={\mathbf{\text{S}}}^{\mathrm{equilibrium}}+\sum _{k=1}^{n}{\mathbf{\text{Q}}}_{k}^{n.\mathrm{equilibrium}}.$$

4.68

The over-stress ${\mathbf{\text{Q}}}_{k}^{n.\mathrm{equilibrium}}$ are work-conjugate to **E*** _{k}*. Drawing a parallel with rheological elements of linear viscoelasticity, and instead of formulating evolution equations for the strain-like variables

Over the course of a life time, human skin undergoes a series of chemo-mechanobiological alterations (i.e. intrinsic ageing) which occur in combination with the effects of external environmental factors such as solar radiations or chemical pollution (i.e. extrinsic ageing). These microstructural and material biophysical changes typically translate into the formation of wrinkles which become more pronounced as intrinsic and extrinsic ageing progress [158]. Skin wrinkles have been the subject of much attention because of their societal importance in relation to age and beauty but also, in the case of expression wrinkles, because of their role in conveying emotions, and thus, as a vehicle for communication. Besides these important aspects of human life, unveiling the underlying mechanical principles that condition the morphologies and patterns of wrinkles are essential in evaluating, and ultimately, predicting, how an aged skin interacts with its environment. Many products that interact with, support or protect human skin actually cause discomfort for the user and even irritation or damage to the skin. Examples include incontinence products, medical devices and sports equipment. Through adhesive and deformation-induced friction, skin wrinkles play a central role in these contact interactions, [21,68,159] and any alteration of the skin surface (e.g. wrinkles) is likely to modulate these effects. Moreover, the growing market of stretchable and wearable epidermal electronics [160] put wrinkles in the spotlight as these systems, do not only interact with the skin and its wrinkles, but are also often multi-layer systems themselves which are prone to wrinkling behaviour [161].

Here, one should provide a slightly more nuanced definition of skin wrinkles. *Skin wrinkles* (or skin folds) [162] are generally associated with ageing and manifest as amplification of the natural skin microrelief (see electronic supplementary material, §S1.7 and figure S10) which is believed to be the results of material and structural alterations of the skin internal structure [158] in combination with changes in the mechanical environment of the body (e.g. relaxation of tension lines and resorption of hypodermal tissue).

*Expression (or temporary) wrinkles* (also called *expression lines*) are associated with skin and facial movement. They can originate from facial muscular activation (i.e. smiling) or by simple mechanical actions on the skin surface such as twisting, shear or compression. *Gravitational folds* on the face arise from the combined effect of constant gravitational loads applied to a microstructurally ageing-altered epidermis, dermis and hypodermis.

Wrinkling, ubiquitous in nature, is a multi-scale spatial phenomenon which can span over eight orders of magnitude [163]. From a physical point of view, wrinkles are the result of a complex interplay between material and structural properties, boundary and loading conditions, the exact nature of which remains to be elucidated [164], particularly for nonlinear materials assembled into geometrically complex structures such as the skin. Although the study of surface instabilities dates back to the seminal work of Biot [165,166], Cerda and Mahadevan revisited the Föppl–von Kármán equations for thin elastic plates using asymptotic and scaling arguments and established a simple theory of wrinkling, valid far from the onset of surface instability. They showed that the wavelength of wrinkles *λ* is proportional to *K*^{−0.25} and *A* where *K* is an equivalent stiffness arising from an ‘elastic substrate’ effect and *A* is the amplitude of wrinkles. Wrinkles can be induced by a variety of ways including residual strains in bi-layer structures [167], differential growth in multi-layer biological structures [168] ( figure 1, the convoluted geometry of the basement membrane separating the viable epidermis from the papillary dermis can be explained by differential growth), compressive and tensile [169] loads combined with geometric constraints (i.e. inextensibility, geometrical coupling with a thicker foundation) in very thin structures (i.e. membranes and shells) or any rebalancing of material and structural properties in these single or multi-layer systems. A simple conceptual model to explain the formation of wrinkles consists of a bi-layer plate structure featuring a thin and stiff layer laying on a much thicker and softer foundation [164]. Upon in-plane compression of both layers, due to its higher bending stiffness, the top (thin) layer will favour large wavelength bending deformations while the bottom (soft) layer will penalize these wavelengths and rather promote smaller wavelengths. As a result of this energy-minimizing structural response to competition between bending energies of the top and bottom layers, intermediate wavelengths, and so, length scales, will emerge under the form of wrinkles. If compression continues after the formation of wrinkles, the latter will evolve into folds and, finally, creases when self-contact of the surface occurs. Since the paper of Cerda & Mahadevan [164] a large body of work on surface instabilities has been being produced, particularly in the physics community, as understanding and harnessing surface instabilities opens up a broad range of industrial applications from nano-metrology to programmable surfaces, whereby generation of self-organized wrinkling patterns can be controlled by judicious structural/material arrangements and tuned localized loading conditions [163,170]. Holland *et al.* [171] recently provided an in-depth analysis of instabilities of soft films on compliant substrates and showed the high sensitivity of these soft bi-layer systems to stiffness ratio, boundary conditions and modes of compression. They also introduced novel metrics, the *effective strain*, *effective stiffness* and *effective wavelength* and recommend their standardized use in the analysis of instabilities.

In the last decade, there has also been a growing interest in studying surface instabilities (e.g. wrinkles, folds and creases) in biological systems as they play a central role in many physiological functions and processes including morphogenesis/growth of normal and cancerous tissues as well as ageing [26,168,172–177]. Of particular relevance to skin and ageing, skin wrinkles constitute a rich test-bed for many theories of surface instability phenomena. However, the majority of resulting mathematical models are often restricted to idealized geometries (e.g. perfectly flat structures, i.e. plates), simplified boundary and loading conditions or simplistic constitutive models.

Because the skin can sustain finite deformations and strains *in vivo*, it is sensible to model it using appropriate corresponding theories. Although precluding material anisotropic effects, a neo-Hookean constitutive formulation is probably the simplest form of materially stable hyperelasticity valid for finite kinematics. In the physics of wrinkling, the ratio of stiffness and that of thickness between the different layers making up a structure are central in conditioning the onset of wrinkle formation as well as wrinkle characteristics [164,167].

The *stratum corneum* is the skin layer most sensitive to external environmental conditions (see §§2.1 and 2.2) and its Young's modulus can span several orders of magnitude as a result of alterations in relative humidity [25]. Considering the skin as an assembly of a thin layer of thickness *h* (i.e. the *stratum corneum*) bonded to an underlying infinitely deep half-space, both made of a neo-Hookean material with distinct mechanical properties (shear moduli *μ*_{f} and *μ*_{s} for, respectively, the thin film and thick substrate) is a reasonable assumption to study the characteristics of wrinkles as a function of layer stiffness ratio. Cao & Hutchinson [167] proposed such a model which was based on an extension of Biot's exact finite strain bifurcation analysis to bi-layer systems. Assuming isochoric behaviour and application of plane-strain compression to both layers, the compression bifurcation strain is

$${\epsilon}_{0}=\frac{1}{4}\sqrt[3]{{\left(\frac{3{\mu}_{\mathrm{s}}}{{\mu}_{\mathrm{f}}}\right)}^{2}}.$$

5.1

They showed that, for both layers in a stress-free reference state, the compressive bifurcation strain at the onset of wrinkling matches very well that determined by Allen [178] who considered the film to be linear elastic provided the modulus ratio *μ*_{f}/*μ*_{s} is large enough (greater than 10). The critical wave number is then determined as:

$${k}_{0}=\frac{1}{h}\sqrt[3]{\frac{3{\mu}_{\mathrm{s}}}{{\mu}_{\mathrm{f}}},}$$

5.2

leading to a wrinkle wavelength *λ*_{0}:

$${\lambda}_{0}=\frac{2\pi}{{k}_{0}}=2\pi h\sqrt[3]{\frac{{\mu}_{\mathrm{f}}}{3{\mu}_{\mathrm{s}}}}.$$

5.3

When the material properties of both layers are identical (*μ*_{f}/*μ*_{s}=1), *ε*_{0}=0.456 which is the solution established by Biot [165] for a homogeneous neo-Hookean half-space. However, based on another study [179], Cao and Hutchinson recognize that, given the high sensitivity of wrinkling to imperfections, the likelihood of the film reaching a strain value *ε*_{0}=0.456 before onset of wrinkling is very low.

Lejeune *et al.* [180] recently developed a comprehensive quad-layer model to shed light on observed geometric instabilities in thin film-substrate systems, thereby accounting for interfacial layers absent from theoretical bi-layer models which do not account for material or structural imperfections from manufacturing processes. Their modelling framework would be applicable to skin provided strains are small and materials could be considered as isotropic linear elastic. Motivated by epidermal electronics applications, Lejeune *et al.* [161] later generalized their previous approach by extending it with an arbitrary number of layers featuring arbitrary thicknesses and by using an algorithmic approach which was validated by testing it against finite-element computations. Through analytical and numerical calculations, Goriely *et al.* [177] studied the role of structural (e.g. thin or thick shell assumption) and constitutive parameters (e.g. type of hyperelastic constitutive formulation) on compression- and extension-induced structural instabilities in elastomers and in biological soft tissues. For certain cases, the typical strain-hardening behaviour of soft tissues was found not to be a stabilizing factor delaying the onset of instabilities. Besides, for a compressed half-space, a soft tissue material is more stable than an elastomeric material while an opposite behaviour is observed when considering a compressed spherical thick shell. Dervaux & Ben Amar [168] proposed a mechanobiological model of constrained growth in bi-layered structures that could explain the wrinkled nature of the basement membrane separating the viable epidermis from the papillary dermis. However, they recognized that the Föppl–von Kármán approximation for the thin top layer is violated when the stiffness of the two layers are of the same order of magnitude.

Ciarletta *et al.* [181] developed an analytical model of shear-induced wrinkling instabilities in skin by assuming the skin to be a mono-layer neo-Hookean solid endowed with a surface energy and subjected to a residual stretch corresponding to a cleavage line (i.e. Langer line). By conducting a perturbation analysis these authors corroborated experimental observations showing that shearing the skin surface along cleavage lines requires a greater magnitude of shear than shearing along the direction of compression. It was also demonstrated that surface energy delays the onset of instability. The model was then extended to transverse isotropy by introducing a single family of fibres (i.e. collagen and elastin) and studied using a Stroh formalism. The mathematical analysis showed that the shear threshold at which surface instabilities occur is lower in the presence of fibres, and further decreases as the apparent stiffness of fibres increase (like that which is often reported in the context of ageing). Motivated by metrology applications, Gower [182] later devised an analytical framework to study wrinkles in soft fibre-reinforced solids in order to establish a link between wrinkle characteristics and material constitutive parameters. He considered a transversely isotropic hyperelastic formulation featuring two distinct invariants, one for representing fibre extension and the other to represent fibre compression. The model showed how the angle between the fibres and the surface wrinkle orientation indicates the mechanical status of fibres: fibres resisting only extension, only compression, or a combination of both. Moreover, at the onset of surface wrinkling, and for the constitutive parameters considered, this angle was shown to be discrete by taking only three or four possible values. Against experimental evidence and physical intuition, Carfagna *et al.* [183] demonstrated that oblique wrinkles (i.e. not aligned with a principal direction of deformation) can arise in soft isotropic solids. Understanding the mechanics of wrinkle formation is an essential step towards the development of applications aiming at reducing or removing wrinkles [184] by straightening of wrinkles [185]. For example, Yu *et al.* [184] recently developed a wearable and physico-chemically tuneable polysiloxane-based material that can be topically applied to the skin surface. After curing, the cross-linked polymer layer bonded to the skin effectively straightens wrinkles. Beyond cosmetics, applications of this technology extend to restoration of the skin barrier function, pharmaceutical delivery and wound dressings.

Although analytical models of skin wrinkles can be extremely insightful in unravelling physical mechanisms that explain experimental observations, they quickly become too restrictive if research objectives are to investigate naturally complex biological structures, non-uniform loading conditions or to account for more sophisticated nonlinear materials. In these cases, numerical techniques such as the ubiquitous finite-element method are the only way forward.

Flynn & McCormack [101] developed an experimentally-based computational model of compression-induced macroscopic skin wrinkles. The model provided a valuable insight into the influence of the number of layers and pre-stress on wrinkle formation. In the *in vitro* and *in silico* experiments, it was observed that the depth of wrinkles significantly increases under 20–30% strain. This model was subsequently used to simulate ageing effects on compressive wrinkle formation but still did not take into account real skin microrelief [103]. Based on the modelling framework of their multi-stage buckling theory which suggested that wrinkle morphology suddenly changes from *stratum corneum* wrinkling to epidermis wrinkling leading, respectively, to shallow fine furrows and deep coarse wrinkles. Kuwazuru *et al.* [186] devised an experimental apparatus to address this question. It was demonstrated that the skin wrinkling rate experiences a step increase at the age of 33 translating into a sudden change in the morphology of skin wrinkles. What these authors called ‘buckling mode switch’ is the result of what had already been explained by Cerda & Mahadevan [164], namely the competition between bending energies of individual layers which favour particular wavelengths during wrinkle formation. Any change in the ratio of bending energies (through changes in layer thickness and/or material properties) affect the morphology of wrinkles. Recently, Shiihara *et al.* [187] conducted a series of finite-element analyses to study the effect of skin microrelief on large wrinkles. Their simplified two-dimensional (2D) analyses were based on non-anatomical skin microrelief and dynamic implicit-based buckling-post-buckling procedures [188]. Shiihara *et al.* [187] suggested that fine microrelief at the surface inhibit large wrinkles (i.e. deeper than 0.35mm). Until the recent work of Limbert & Kuhl [189], no physics-based study had examined the role of (anatomically realistic) skin microrelief on the characteristics of surface micro-wrinkles.

Using an anatomically-based bilayer finite-element model featuring neo-Hookean materials, these authors highlighted the critical role of layer stiffness ratio on wrinkle characteristics and surface strains in simulated in-plane compression (figure 4, see also electronic supplementary material, §S1.13). As predicted by theoretical and computational models of growth-induced surface instabilities of flat bi-layer systems [190], period-doubling effects were also observed in the finite-element analyses. The highly nonlinear nature of surface instabilities make their numerical realization very challenging, and standard Newton–Raphson incremental solving procedures might be insufficient [188], thus requiring better path-following procedures such as arc-length methods [189,191] and/or dynamic regularization techniques [189]. Cubic finite-element shape functions continuous across element boundaries are desirable to capture steep curvature gradients such as those occurring in surface instabilities. Isogeometric finite-element formulations [192–194] are therefore superior to those based on Lagrange polynomial interpolation. Identifying critical conditions that trigger surface instabilities is also a challenging endeavour but progress in this area is steadily moving forward [195].

From this review which offered a relatively large and eclectic selection of constitutive models of the skin, it is clear that these models have now reached a high level of sophistication, and can capture a wide range of relevant biophysical features. Nevertheless, as pointed out by Jor *et al.* [6] in their review paper on the computational and experimental characterization of skin, a significant limiting factor in the development and adoption of advanced constitutive theories is the scarcity and *relevance* of captured experimental data. There are several fundamental aspects, or more precisely, reasons and questions associated with this observation:

- (1)Large intra- and inter-individual variability of biophysical properties;
- (2)The
*in vivo*and*ex vivo*biophysical environments of the skin are fundamentally different; - (3)Extreme sensitivity of the skin to environmental conditions; and
- (4)How to best integrate multi-modality and multi-scale imaging, characterization and modelling techniques?

Intra-individually, the skin is a complex heterogeneous adaptive structure which varies according to body location, health status and history, diet, age, lifestyle, *external* environmental conditions (e.g. temperature, humidity, pollution level, water quality, sun exposure, contact with external surfaces) and *internal* environmental conditions (e.g. hormones, pregnancy, water and glucose levels, tension lines). Besides intra-individual variability, there is a strong biophysical variability, particularly that associated with sex and ethnicity. There are, therefore, formidable challenges in *representatively* characterizing the skin ultrastructure as well as its biophysical properties. Should we describe the skin as an individual-specific system, or rather, as a statistical system describing particular populations? Or, even more daringly, as a combination of both? This issue is, of course, not unique to the skin, but pervasive throughout many complex biological systems. Data mining and machine learning techniques [196] are likely to play an increasing role in the future to make sense of large and complex heterogeneous datasets, whether they originate from physical or computer experiments, expert knowledge (e.g. anatomists, clinicians, nurses, vets) or from any other means (e.g. patient's observations, shamanic knowledge).

Perhaps, another related buzz word of modern times to throw in here, is *artificial intelligence* (AI)! AI is simply the processing ability of a (non-human) computer to solve problems without having been taught how. Only when the techniques discussed above will have assisted us in establishing a more fundamental and mechanistic understanding of the multi-factorial nature of skin biophysics, will we be able to seamlessly integrate these huge datasets into, first, *descriptive*, and then, *predictive*, mathematical and computational models. Multi-variate and multi-scale data-based and/or physics-based statistical models of the skin built from the results of machine learning (i.e. meta-models) could then replace computationally expensive physics-based finite-element models, and be used to predict a variety of scenarios and outcomes. For example, one could ask: What should be the optimum stiffness/thickness of a silicon layer on a seating surface so that pressure on the skin does not reach a critically-damaging threshold? How would this change if relative humidity increases by 30%? And if the patient is between 50 and 65 years of age instead of 25? Of course, it depends! And it depends on a lot of factors. The answers to these questions are likely to be statistical distributions rather than single deterministic values. Perhaps, in the same way that statistical information is currently being integrated into constitutive models (e.g. fibre dispersion, degree of crimping of collagen fibres), additional stochastic elements could be added into constitutive equations effectively developing what could be termed *stochastic partial differential equations* [197]. This would be another approach to encompass variability, directly into the constitutive equations, rather than at the input data level when conducting probabilistic simulations.

The reliance of structural and structurally based phenomenological models on accurate descriptions of the skin ultrastructure and the associated mechanical properties of its sub-components, means that imaging and characterization techniques needs to span multiple spatial and temporal scales if one wants to capture more reliably the complex biophysics of the skin. Additionally, and as a corollary to that, multi-scale modelling techniques needs to be developed in concert with these physical experimental protocols. Computational multi-scale methods combining atomistic, molecular and continuum techniques are likely to play an increasingly important role in the modelling of skin biophysics in the forthcoming decades [57,198,199]. Moreover, integration of imaging and physical testing (i.e. *in situ* imaging) [200] is a very promising approach that should be followed, as much can be learnt about microstructural mechanisms under real loading conditions. For a detailed discussion of challenges/opportunities of tissue structure imaging and mechanical characterization, see the review by Jor *et al*. [6]. These authors also discuss the critically important point concerning the differences between *in vivo* and *ex vivo* testing of skin. They rightly point out that, to better translate research into the clinical setting, major attention should be directed towards the development of *in vivo* quantitative characterization techniques. The advent of cheap, powerful personal and cloud computing is rapidly transforming the face of research in the life sciences. In the very near future, on-demand physics-based computing accessible from any personal electronic device such as tablets and smart phones, and from any location on the planet, is likely to become an everyday commodity, like electricity. It is anticipated that patients will be able to collect *in vivo* data, process and transmit it so that health practitioners will be able to use these tools, remotely, or in clinics, to assist their decision process while engineers and scientists will continue to embrace, develop and push forward this trend. In the meantime, biophysical modellers are likely to keep busy!

Click here to view.^{(875K, pdf)}

The idea for this article arose following many fruitful discussions with Prof. Michel Destrade (NUI Galway, Ireland), Prof. T. Christian Gasser (KTH Stockholm, Sweden) and Dr Artur L. Gower (University of Manchester, UK) initiated during the EPSRC-sponsored workshop ‘Constitutive behaviour of soft tissues: connecting experimental and modelling perspectives’, 31 August–2 September 2016, University of Manchester organized by Prof. William Parnell and Dr Tom Shearer (University of Manchester, UK). G.L. thanks them for their invitation to participate in the workshop. G.L. would also like to thank Dr Maria-Fabiola Leyva-Mendivil (University of Southampton) for providing the histological section used in figure 1.

This article has no additional data.

I declare I have no competing interests.

G.L. acknowledges the financial support he has received over the last few years from the Royal Society, The Royal Academy of Engineering, The British High Commission in South Africa, EPSRC, Procter and Gamble, L'Oréal and the US Air Force.

1. Burns T, Breathnach S, Cox N, Griffiths C.
2004.
Rook's textbook of dermatology, 7th edn
Oxford, UK: Blackwell Science.

2. Silver FH, Siperko LM, Seehra GP
2003.
Mechanobiology of force transduction in dermal tissue. Skin Res. Technol.
9, 3–23. (doi:10.1034/j.1600-0846.2003.00358.x) [PubMed]

3. Jablonski NG.
2006.
Skin: a natural history. Berkeley, CA: University of California Press.

4. Belytschko T, Liu WK, Moran B
2000.
Nonlinear finite elements for continua and structures. Oxford, UK: Wiley.

5. Forrester AIJ.
2010.
Black-box calibration for complex systems simulation. Phil. Trans. R. Soc. A
368, 3567–3579. (doi:10.1098/rsta.2010.0051) [PubMed]

6. Jor JWY, Parker MD, Taberner AJ, Nash MP, Nielsen PMF
2013.
Computational and experimental characterization of skin mechanics: identifying current challenges and future directions. Wiley Interdiscip. Rev. Syst. Biol. Med.
5, 539–556. (doi:10.1002/wsbm.1228) [PubMed]

7. Li W.
2015.
Modelling methods for *in vitro* biomechanical properties of the skin: a review. Biomed. Eng. Lett.
5, 241–250. (doi:10.1007/s13534-015-0201-3)

8. Fung YC.
1981.
Biomechanics: mechanical properties of living tissues. New York, NY: Springer.

9. Humphrey JD.
2003.
Continuum biomechanics of soft biological tissues. Proc. R. Soc. Lond. A
459, 3–46. (doi:10.1098/rspa.2002.1060)

10. Lanir Y.
2016.
Multi-scale structural modeling of soft tissues mechanics and mechanobiology. J. Elast.
1–42. (doi:10.1007/s10659-016-9607-0)

11. Shimizu H.
2007.
Shimizu's textbook of dermatology. Hokkaido, Japan: Hokkaido University Press - Nakayama Shoten Publishers.

12. Buganza Tepole A, Kuhl E
2016.
Computational modeling of chemo-bio-mechanical coupling: a systems-biology approach toward wound healing. Comput. Methods Biomech. Biomed. Eng.
19, 13–30. (doi:10.1080/10255842.2014.980821) [PubMed]

13. Kvistedal YA, Nielsen PMF
2009.
Estimating material parameters of human skin *in vivo*. Biomech. Model Mechanobiol.
8, 1–8. (doi:10.1007/s10237-007-0112-z) [PubMed]

14. Lanir Y.
1987.
Skin mechanics. In Handbook of bioengineering (eds Skalak R, Chien S). New York, NY: McGraw-Hill.

15. Vierkötter A, Krutmann J
2012.
Environmental influences on skin aging and ethnic-specific manifestations. Dermato-endocrinology
4, 227–231. (doi:10.4161/derm.19858) [PMC free article] [PubMed]

16. Silver FH, Freeman JW, DeVore D
2001.
Viscoelastic properties of human skin and processed dermis. Skin Res. Technol.
7, 18–23. (doi:10.1034/j.1600-0846.2001.007001018.x) [PubMed]

17. Limbert G.
2014.
Chapter 4: State-of-the-art constitutive models of skin biomechanics. In Computational biophysics of the skin (ed. Querleux B, editor. ), pp. 95–131. Singapore: Pan Stanford Publishing Pte. Ltd.

18. Marieb EN, Hoehn K
2010.
Human anatomy and physiology, 8th edn
San Francisco, CA: Pearson International Edition.

19. Chan LS.
1997.
Human skin basement membrane in health and autoimmune diseases. Front. Biosci.
2, 343–352. (doi:10.2741/A196) [PubMed]

20. Leyva-Mendivil MF, Page A, Bressloff NW, Limbert G
2015.
A mechanistic insight into the mechanical role of the stratum corneum during stretching and compression of the skin. J. Mech. Behav. Biomed. Mater.
49, 197–219. (doi:10.1016/j.jmbbm.2015.05.010) [PubMed]

21. Leyva-Mendivil MF, Lengiewicz J, Page A, Bressloff NW, Limbert G
2017.
Skin microstructure is a key contributor to its friction behaviour. Tribology Lett.
65, 12 (doi:10.1007/s11249-016-0794-4)

22. Dandekar K, Raju BI, Srinivasan MA
2003.
3-D finite-element models of human and monkey fingertips to investigate the mechanics of tactile sense. J. Biomech. Eng.
125, 682–691. (doi:10.1115/1.1613673) [PubMed]

23. Xu F, Lu T
2011.
Introduction to skin biothermomechanics and thermal pain. Heidelberg, Germany: Springer.

24. Biniek K, Levi K, Dauskardt RH
2012.
Solar UV radiation reduces the barrier function of human skin. Proc. Natl Acad. Sci USA
109, 17111–17116. (doi:10.1073/pnas.1206851109) [PubMed]

25. Wu KS, van Osdol WW, Dauskardt RH
2006.
Mechanical properties of human stratum corneum: effects of temperature, hydration, and chemical treatment. Biomaterials
27, 785–795. (doi:10.1016/j.biomaterials.2005.06.019) [PubMed]

26. Ciarletta P, Ben Amar M
2012.
Papillary networks in the dermal-epidermal junction of skin: a biomechanical model. Mech. Res. Commun.
42, 68–76. (doi:10.1016/j.mechrescom.2011.12.001)

27. Burgeson RE, Christiano AM
1997.
The dermal-epidermal junction. Curr. Opin. Cell Biol.
9, 651–658. (doi:10.1016/S0955-0674(97)80118-4) [PubMed]

28. Ribeiro JF, dos Anjos EHM, Mello MLS, de Campos Vidal B
2013.
Skin collagen fiber molecular order: a pattern of distributional fiber orientation as assessed by optical anisotropy and image analysis. PLoS ONE
8, e54724 (doi:10.1371/journal.pone.0054724) [PMC free article] [PubMed]

29. Gosline J, Lillie M, Carrington E, Guerette P, Ortlepp C, Savage K
2002.
Elastic proteins: biological roles and mechanical properties. Phil. Trans. R. Soc. Lond. B
357, 121–132. (doi:10.1098/rstb.2001.1022) [PMC free article] [PubMed]

30. Sherratt MJ.
2013.
Age-related tissue stiffening: cause and effect. Adv. Wound Care
2, 11–17. (doi:10.1089/wound.2011.0328) [PMC free article] [PubMed]

31. Langer K.
1861.
Zur Anatomie und Physiologie der Haut. Über die Spaltbarkeit der Cutis. Sitzungsbericht der Mathematisch-naturwissenschaftlichen Classe der Wiener Kaiserlichen Academie der Wissenschaften Abt 44.

32. Langer K.
1978.
On the anatomy and physiology of the skin: I. The cleavability of the cutis. Br. J. Plastic Surg.
31, 3–8. (doi:10.1016/0007-1226(78)90003-6) [PubMed]

33. Langer K.
1978.
On the anatomy and physiology of the skin: II. Skin tension (with 1 figure). Br. J. Plastic Surg.
31, 93–106. (doi:10.1016/S0007-1226(78)90109-1) [PubMed]

34. Ní Annaidh A, Bruyère K, Destrade M, Gilchrist MD, Otténio M
2011.
Characterization of the anisotropic mechanical properties of excised human skin. J. Mech. Behav. Biomed. Mater.
5, 139–148. (doi:10.1016/j.jmbbm.2011.08.016) [PubMed]

35. Alexander H, Cook TH
1977.
Accounting for natural tension in the mechanical testing of human skin. J. Invest. Dermatol.
69, 310–314. (doi:10.1111/1523-1747.ep12507731) [PubMed]

36. Flynn C, Taberner AJ, Nielsen PMF, Fels S
2013.
Simulating the three-dimensional deformation of *in vivo* facial skin. J. Mech. Behav. Biomed. Mater.
28, 484–494. (doi:10.1016/j.jmbbm.2013.03.004) [PubMed]

37. Flynn C, Stavness I, Lloyd J, Fels S
2013.
A finite element model of the face including an orthotropic skin model under *in vivo* tension. Comput. Methods Biomech. Biomed. Eng.
18, 571–582. (doi:10.1080/10255842.2013.820720) [PubMed]

38. Deroy C, Destrade M, Mc Alinden A, Ní Annaidh A
2016.
Non-invasive evaluation of skin tension lines with elastic waves. Skin Res. Technol. 23, 326–335 (doi:10.1111/srt.12339) [PubMed]

39. Rosado C, Antunes F, Barbosa R, Fernando R, Estudante M, Silva HN, Rodrigues LM
2016.
About the *in vivo* quantitation of skin anisotropy. Skin Res. Technol. 23, 429–436 (doi:10.1111/srt.12353) [PubMed]

40. Wan Abas WAB.
1994.
Biaxial tension test of human skin *in vivo*. Biomed. Mater. Eng.
4, 473–486. [PubMed]

41. Ní Annaidh A, Bruyère K, Destrade M, Gilchrist MD, Maurini C, Otténio M, Saccomandi G
2012.
Automated estimation of collagen fibre dispersion in the dermis and its contribution to the anisotropic behaviour of skin. Ann. Biomed. Eng.
40, 1666–1678. (doi:10.1007/s10439-012-0542-3) [PubMed]

42. Ottenio M, Tran D, Annaidh AN, Gilchrist MD, Bruyère K
2015.
Strain rate and anisotropy effects on the tensile failure characteristics of human skin. J. Mech. Behav. Biomed. Mater.
41, 241–250. (doi:10.1016/j.jmbbm.2014.10.006) [PubMed]

43. Kvistedal YA, Nielsen PMF
2004.
Investigating stress-strain properties of *in-vivo* human skin using multiaxial loading experiments and finite element modeling. In *Proc. of the 26th Annu. Int. Conf. of the IEEE Engineering in Medicine and Biology Society, San Francisco, CA, 1--4 September*, vol. 7, pp. 5096–5099. Piscataway, NJ: IEEE.

44. Batisse D, Bazin R, Baldeweck T, Querleux B, Lévêque JL
2002.
Influence of age on the wrinkling capacities of skin. Skin Res. Technol.
8, 148–154. (doi:10.1034/j.1600-0846.2002.10308.x) [PubMed]

45. Delalleau A, Josse G, Lagarde JM, Zahouani H, Bergheau JM
2006.
Characterization of the mechanical properties of skin by inverse analysis combined with the indentation test. J. Biomech.
39, 1603–1610. (doi:10.1016/j.jbiomech.2005.05.001) [PubMed]

46. Diridollou S, Patat F, Gens F, Vaillant L, Black D, Lagarde JM, Gall Y, Berson M
2000.
*In vivo* model of the mechanical properties of the human skin under suction. Skin Res. Technol.
6, 214–221. (doi:10.1034/j.1600-0846.2000.006004214.x) [PubMed]

47. Dobrev HQ.
2000.
Use of Cutometer to assess epidermal hydration. Skin Res. Technol.
6, 239–244. (doi:10.1034/j.1600-0846.2000.006004239.x) [PubMed]

48. Hendriks FM, Brokken D, van Eemeren J, Oomens CWJ, Baaijens FPT, Horsten J
2003.
A numerical-experimental method to characterize the non-linear mechanical behaviour of human skin. Skin Res. Technol.
9, 274–283. (doi:10.1034/j.1600-0846.2003.00019.x) [PubMed]

49. Weickenmeier J, Jabareen M, Mazza E
2015.
Suction based mechanical characterization of superficial facial soft tissues. J. Biomech.
48, 4279–4286. (doi:10.1016/j.jbiomech.2015.10.039) [PubMed]

50. Tonge TK, Atlan LS, Voo LM, Nguyen TD
2013.
Full-field bulge test for planar anisotropic tissues: part I—experimental methods applied to human skin tissue. Acta Biomater.
9, 5913–5925. (doi:10.1016/j.actbio.2012.11.035) [PubMed]

51. Geerligs M, Oomens CWJ, Ackermans PAJ, Baaijens FPT, Peters GWM
2011.
Linear shear response of the upper skin layers. Biorheology
48, 229–245. (doi:10.3233/BIR-2011-0590) [PubMed]

52. Geerligs M, van Breemen LCA, Peters GWM, Ackermans PAJ, Baaijens FPT, Oomens CWJ
2011.
*In vitro* indentation to determine the mechanical properties of epidermis. J. Biomech.
44, 1176–1181. (doi:10.1016/j.jbiomech.2011.01.015) [PubMed]

53. Lamers E, van Kempen THS, Baaijens FPT, Peters GWM, Oomens CWJ
2013.
Large amplitude oscillatory shear properties of human skin. J. Mech. Behav. Biomed. Mater.
28, 462–470. (doi:10.1016/j.jmbbm.2013.01.024) [PubMed]

54. Lanir Y, Fung YC
1974.
Two-dimensional mechanical properties of rabbit skin—II: experimental results. J. Biomech.
7, 171–182. (doi:10.1016/0021-9290(74)90058-X) [PubMed]

55. Wong WLE, Joyce TJ, Goh KL
2016.
Resolving the viscoelasticity and anisotropy dependence of the mechanical properties of skin from a porcine model. Biomech. Model. Mechanobiol.
15, 433–446. (doi:10.1007/s10237-015-0700-2) [PubMed]

56. Veronda DR, Westmann R
1970.
Mechanical characterization of skin—finite deformations. J. Biomech.
3, 111–124. (doi:10.1016/0021-9290(70)90055-2) [PubMed]

57. Marino M.
2016.
Molecular and intermolecular effects in collagen fibril mechanics: a multiscale analytical model compared with atomistic and experimental studies. Biomech. Model. Mechanobiol.
15, 133–154. (doi:10.1007/s10237-015-0707-8) [PubMed]

58. Spencer AJM.
1984.
Constitutive theory for strongly anisotropic solids. In Continuum theory of the mechanics of fibre-reinforced composites (ed. Spencer AJM, editor. ), pp. 1–32. Vienna, Austria: Springer.

59. Šolinc U, Korelc J
2015.
A simple way to improved formulation of FE^{2} analysis. Comput. Mech.
56, 905–915. (doi:10.1007/s00466-015-1208-4)

60. Saeb S, Steinmann P, Javili A
2016.
Aspects of computational homogenization at finite deformations: a unifying review from Reuss' to Voigt's bound. Appl. Mech. Rev.
68, 050801 (doi:10.1115/1.4034024)

61. Buganza Tepole A, Gart M, Gosain AK, Kuhl E
2014.
Characterization of living skin using multi-view stereo and isogeometric analysis. Acta Biomater.
10, 4822–4831. (doi:10.1016/j.actbio.2014.06.037) [PMC free article] [PubMed]

62. Buganza Tepole A, Gosain AK, Kuhl E
2012.
Stretching skin: the physiological limit and beyond. Int. J. Non-Linear Mech.
47, 938–949. (doi:10.1016/j.ijnonlinmec.2011.07.006) [PMC free article] [PubMed]

63. Buganza Tepole A, Gosain AK, Kuhl E
2014.
Computational modeling of skin: Using stress profiles as predictor for tissue necrosis in reconstructive surgery. Comput. Struct.
143, 32–39. (doi:10.1016/j.compstruc.2014.07.004) [PMC free article] [PubMed]

64. Buganza Tepole A, Ploch CJ, Wong J, Gosain AK, Kuhl E
2011.
Growing skin: a computational model for skin expansion in reconstructive surgery. J. Mech. Phys. Solids
59, 2177–2190. (doi:10.1016/j.jmps.2011.05.004) [PMC free article] [PubMed]

65. Zöllner A, Buganza Tepole A, Gosain AK, Kuhl E
2012.
Growing skin: tissue expansion in pediatric forehead reconstruction. Biomech. Model Mechanobiol.
11, 855–867. (doi:10.1007/s10237-011-0357-4) [PMC free article] [PubMed]

66. Zöllner AM, Buganza Tepole A, Kuhl E
2012.
On the biomechanics and mechanobiology of growing skin. J. Theor. Biol.
297, 166–175. (doi:10.1016/j.jtbi.2011.12.022) [PMC free article] [PubMed]

67. Zöllner AM, Holland MA, Honda KS, Gosain AK, Kuhl E
2013.
Growth on demand: reviewing the mechanobiology of stretched skin. J. Mech. Behav. Biomed. Mater.
28, 495–509. (doi:10.1016/j.jmbbm.2013.03.018) [PMC free article] [PubMed]

68. Leyva-Mendivil MF, Lengiewicz J, Page A, Bressloff NW, Limbert G
2017.
Implications of multi-asperity contact for shear stress distribution in the viable epidermis—an image-based finite element study. Biotribology (doi:10.1016/j.biotri.2017.04.001)

69. Young PG, Beresford-West TBH, Coward SRL, Notarberardino B, Walker B, Abdul-Aziz A
2008.
An efficient approach to converting three-dimensional image data into highly accurate computational models. Phil. Trans. R. Soc. A
366, 3155–3173. (doi:10.1098/rsta.2008.0090) [PubMed]

70. Limbert G, van Lierde C, Muraru OL, Walboomers XF, Frank M, Hansson S, Middleton J, Jaecques S
2010.
Trabecular bone strains around a dental implant and associated micromotions--a micro-CT-based three-dimensional finite element study. J. Biomech.
43, 1251–1261. (doi:10.1016/j.jbiomech.2010.01.003) [PubMed]

71. Linder-Ganz E, Shabshin N, Itzchak Y, Gefen A
2007.
Assessment of mechanical conditions in sub-dermal tissues during sitting: A combined experimental-MRI and finite element approach. J. Biomech.
40, 1443–1454. (doi:10.1016/j.jbiomech.2006.06.020) [PubMed]

72. Limbert G, Bryan R, Cotton R, Young P, Hall-Stoodley L, Kathju S, Stoodley P
2013.
On the mechanics of bacterial biofilms on non-dissolvable surgical sutures: a laser scanning confocal microscopy-based finite element study. Acta Biomater.
9, 6641–6652. (doi:10.1016/j.actbio.2013.01.017) [PubMed]

73. Gerhardt LC, Strässle V, Lenz A, Spencer ND, Derler S
2008.
Influence of epidermal hydration on the friction of human skin against textiles. J. R. Soc. Interface
5, 1317–1328. (doi:10.1098/rsif.2008.0034) [PMC free article] [PubMed]

74. Adams MJ, Briscoe BJ, Johnson SA
2007.
Friction and lubrication of human skin. Tribology Lett.
26, 239–253. (doi:10.1007/s11249-007-9206-0)

75. Derler S, Gerhardt LC, Lenz A, Bertaux E, Hadad M
2009.
Friction of human skin against smooth and rough glass as a function of the contact pressure. Tribology Int.
42, 1565–1574. (doi:10.1016/j.triboint.2008.11.009)

76. Kwiatkowska M, Franklin SE, Hendriks CP, Kwiatkowski K
2009.
Friction and deformation behaviour of human skin. Wear
267, 1264–1273. (doi:10.1016/j.wear.2008.12.030)

77. Wolfram LJ.
1983.
Friction of skin. J. Soc. Cosmet. Chem.
34, 465–476.

78. Stupkiewicz S, Lewandowski MJ, Lengiewicz J
2014.
Micromechanical analysis of friction anisotropy in rough elastic contacts. Int. J. Solids Struct.
51, 3931–3943. (doi:10.1016/j.ijsolstr.2014.07.013)

79. Goldstein B, Sanders J
1998.
Skin response to repetitive mechanical stress: a new experimental model in pig. Arch. Phys. Med. Rehabil.
79, 265–272. (doi:10.1016/S0003-9993(98)90005-3) [PubMed]

80. Weickenmeier J, Jabareen M
2014.
Elastic–viscoplastic modeling of soft biological tissues using a mixed finite element formulation based on the relative deformation gradient. Int. J. Numer. Methods Biomed. Eng.
30, 1238–1262. (doi:10.1002/cnm.2654) [PubMed]

81. Li W, Luo XY
2016.
An invariant-based damage model for human and animal skins. Ann. Biomed. Eng.
44, 3109–3122. (doi:10.1007/s10439-016-1603-9) [PMC free article] [PubMed]

82. McBride A, Bargmann S, Pond D, Limbert G
2016.
Thermoelastic modelling of the skin at finite deformations. J. Therm. Biol
62, 201–209. (doi:10.1016/j.jtherbio.2016.06.017) [PubMed]

83. Vermolen FJ, Gefen A, Dunlop JWC
2012.
*In vitro* ‘wound’ healing: experimentally based phenomenological modeling. Adv. Eng. Mater.
14, 76–88. (doi:10.1002/adem.201180080)

84. Sherratt JA, Dallon JC
2002.
Theoretical models of wound healing: past successes and future challenges. C. R. Biol.
325, 557–564. (doi:10.1016/S1631-0691(02)01464-6) [PubMed]

85. Buganza Tepole A.
2017.
Computational systems mechanobiology of wound healing. Comput. Methods Appl. Mech. Eng.
314, 46–70. (doi:10.1016/j.cma.2016.04.034)

86. Marsden JE, Hughes TJR
1994.
Mathematical foundations of elasticity. New York, NY: Dover.

87. Holzapfel GA.
2000.
Nonlinear solid mechanics. A continuum approach for engineering. Chichester, UK: John Wiley & Sons.

88. Boehler
1978.
Lois de comportement anisotrope des milieux continus. J. Mécanique
17, 153–190.

89. Limbert G, Taylor M
2002.
On the constitutive modeling of biological soft connective tissues. A general theoretical framework and tensors of elasticity for strongly anisotropic fiber-reinforced composites at finite strain. Int. J. Solids Struct.
39, 2343–2358. (doi:10.1016/S0020-7683(02)00084-7)

90. Spencer AJM.
1992.
Continuum theory of the mechanics of fibre-reinforced composites. New York, NY: Springer.

91. Criscione JC, Humphrey JD, Douglas AS, Hunter WC
2000.
An invariant basis for natural strain which yields orthogonal stress response terms in isotropic hyperelasticity. J. Mech. Phys. Solids
48, 2445–2465. (doi:10.1016/S0022-5096(00)00023-5)

92. Criscione JC, Douglas AS, Hunter WC
2001.
Physically based strain invariant set for materials exhibiting transversely isotropic behavior. J. Mech. Phys. Solids
49, 871–897. (doi:10.1016/S0022-5096(00)00047-8)

93. Lanir Y.
1983.
Constitutive equations for fibrous connective tissues. J. Biomech.
16, 1–22. (doi:10.1016/0021-9290(83)90041-6) [PubMed]

94. Meijer R, Douven LFA, Oomens CWJ
1999.
Characterisation of anisotropic and non-linear behaviour of human skin *in vivo*. Comput. Methods Biomech. Biomed. Eng.
2, 13–27. (doi:10.1080/10255849908907975) [PubMed]

95. Jor JWY, Nash MP, Nielsen PMF, Hunter PJ
2011.
Estimating material parameters of a structurally based constitutive relation for skin mechanics. Biomech. Model. Mechanobiol.
10, 767–778. (doi:10.1007/s10237-010-0272-0) [PubMed]

96. Holzapfel GA, Ogden RW
2016.
On fiber dispersion models: exclusion of compressed fibers and spurious model comparisons. J. Elast., 1–20. (doi:10.1007/s10659-016-9605-2)

97. Gasser TC, Ogden RW, Holzapfel GA
2006.
Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. J. R. Soc. Interface.
3, 15–35. (doi:10.1098/rsif.2005.0073) [PMC free article] [PubMed]

98. Arruda EM, Boyce MC
1993.
A three-dimensional constitutive model for the large stretch behavior of rubber elastic-materials. J. Mech. Phys. Solids
41, 389–412. (doi:10.1016/0022-5096(93)90013-6)

99. Bischoff JE, Arruda EA, Grosh K
2002.
A microstructurally based orthotropic hyperelastic constitutive law. J. Appl. Mech.
69, 570–579. (doi:10.1115/1.1485754)

100. Bischoff JE, Arruda EM, Grosh K
2000.
Finite element modeling of human skin using an isotropic, nonlinear elastic constitutive model. J. Biomech.
33, 645–652. (doi:10.1016/S0021-9290(00)00018-X) [PubMed]

101. Flynn C, McCormack BAO
2008.
Finite element modelling of forearm skin wrinkling. Skin Res. Technol.
14, 261–269. (doi:10.1111/j.1600-0846.2008.00289.x) [PubMed]

102. Flynn C, McCormack BAO
2008.
A simplified model of scar contraction. J. Biomech.
41, 1582–1589. (doi:10.1016/j.jbiomech.2008.02.024) [PubMed]

103. Flynn CO, McCormack BAO
2010.
Simulating the wrinkling and aging of skin with a multi-layer finite element model. J. Biomech.
43, 442–448. (doi:10.1016/j.jbiomech.2009.10.007) [PubMed]

104. Kuhl E, Garikipati K, Arruda E, Grosh K
2005.
Remodeling of biological tissue: Mechanically induced reorientation of a transversely isotropic chain network. J. Mech. Phys. Solids
53, 1552–1573. (doi:10.1016/j.jmps.2005.03.002)

105. Limbert G, Middleton J
2005.
A polyconvex anisotropic strain energy function. Application to soft tissue mechanics. In *ASME Summer Bioengineering Conf., Vail, CO, 22--26 June*. Piscataway, NJ: IEEE.

106. Itskov M, Ehret AE, Mavrilas D
2006.
A polyconvex anisotropic strain-energy function for soft collagenous tissues. Biomech. Model Mechanobiol.
5, 17–26. (doi:10.1007/s10237-005-0006-x) [PubMed]

107. Itskov M, Aksel N
2004.
A class of orthotropic and transversely isotropic hyperelastic constitutive models based on a polyconvex strain energy function. Int. J. Solids Struct.
41, 3833–3848. (doi:10.1016/j.ijsolstr.2004.02.027)

108. Holzapfel GA, Gasser TC, Ogden RW
2000.
A new constitutive framework for arterial wall mechanics and a comparative study of material models. J. Elast.
61, 1–48. (doi:10.1023/A:1010835316564)

109. Tonge TK, Voo LM, Nguyen TD
2013.
Full-field bulge test for planar anisotropic tissues: Part II – A thin shell method for determining material parameters and comparison of two distributed fiber modeling approaches. Acta Biomater.
9, 5926–5942. (doi:10.1016/j.actbio.2012.11.034) [PubMed]

110. Flynn C, Rubin MB, Nielsen P
2011.
A model for the anisotropic response of fibrous soft tissues using six discrete fibre bundles. Int. J. Numer. Methods Biomed. Eng.
27, 1793–1811. (doi:10.1002/cnm.1440)

111. Ankersen J, Birkbeck AE, Thomson RD, Vanezis P
1999.
Puncture resistance and tensile strength of skin simulants. Proc. Inst. Mech. Eng. H J. Eng. Med.
213, 493–501. (doi:10.1243/0954411991535103) [PubMed]

112. Flynn C, Rubin MB
2012.
An anisotropic discrete fibre model based on a generalised strain invariant with application to soft biological tissues. Int. J. Eng. Sci.
60, 66–76. (doi:10.1016/j.ijengsci.2012.04.006)

113. Limbert G.
2011.
A mesostructurally-based anisotropic continuum model for biological soft tissues--Decoupled invariant formulation. J. Mech. Behav. Biomed. Mater.
4, 1637–1657. (doi:10.1016/j.jmbbm.2011.07.016) [PubMed]

114. Bischoff JE, Arruda EM, Grosh K
2002.
Finite element simulations of orthotropic hyperelasticity. Finite Elem. Anal. Des.
38, 983–998. (doi:10.1016/S0168-874X(02)00089-6)

115. Lu J, Zhang L
2005.
Physically motivated invariant formulation for transversely isotropic hyperelasticity. Int. J. Solids Struct.
42, 6015–6031. (doi:10.1016/j.ijsolstr.2005.04.014)

116. Kratky O, Porod G
1949.
Röntgenuntersuchungen gelöster Fadenmoleküle. Rec. Trav. Chim. Pays-Bas. Belg.
68, 1106–1122. (doi:10.1002/recl.19490681203)

117. Korelc J, Šolinc U, Wriggers P
2010.
An improved EAS brick element for finite deformation. Comput. Mech.
46, 641–659. (doi:10.1007/s00466-010-0506-0)

118. Gautieri A, Vesentini S, Redaelli A, Buehler MJ
2011.
Hierarchical structure and nanomechanics of collagen microfibrils from the atomic scale up. Nano Lett.
11, 757–766. (doi:10.1021/nl103943u) [PubMed]

119. Sun YL, Luo ZP, Fertala A, An KN
2002.
Direct quantification of the flexibility of type I collagen monomer. Biomech. Biophys. Res. Commun.
295, 382–386. (doi:10.1016/S0006-291X(02)00685-X) [PubMed]

120. Groves RB, Coulman SA, Birchall JC, Evans SL
2013.
An anisotropic, hyperelastic model for skin: Experimental measurements, finite element modelling and identification of parameters for human and murine skin. J. Mech. Behav. Biomed. Mater.
18, 167–180. (doi:10.1016/j.jmbbm.2012.10.021) [PubMed]

121. Weiss JA, Maker BN, Govindjee S
1996.
Finite element implementation of incompressible transversely isotropic hyperelasticity. Comput. Methods Appl. Mech. Eng.
135, 107–128. (doi:10.1016/0045-7825(96)01035-3)

122. Barbenel JC, Evans JH
1973.
The time-dependent mechanical properties of skin. J. Invest. Dermatol.
69, 165–172. [PubMed]

123. Pereira JM, Mansour JM, Davis BR
1990.
Analysis of shear-wave propagation in skin; application to an experimental procedure. J. Biomech.
23, 745–751. (doi:10.1016/0021-9290(90)90021-T) [PubMed]

124. Pereira JM, Mansour JM, Davis BR
1991.
Dynamic measurement of the viscoelastic properties of skin. J. Biomech.
24, 157–162. (doi:10.1016/0021-9290(91)90360-Y) [PubMed]

125. Lanir Y.
1979.
The rheological behavior of the skin: experimental results and a structural model. Biorheology
16, 191–202. [PubMed]

126. Wu JZ, Cutlip RG, Welcome D, Dong RG
2006.
Estimation of the viscous properties of skin and subcutaneous tissue in uniaxial stress relaxation tests. Bio-Med. Mater. Eng.
16, 53–66. [PubMed]

127. Khatyr F, Imberdis C, Vescovo P, Varchon D, Lagarde JM
2004.
Model of the viscoelastic behaviour of skin *in vivo* and study of anisotropy. Skin Res. Technol.
10, 96–103. (doi:10.1111/j.1600-0846.2004.00057.x) [PubMed]

128. Boyer G, Laquieze L, Le Bot A, Laquieze S, Zahouani H
2009.
Dynamic indentation on human skin in vivo: ageing effects. Skin Res. Technol.
15, 55–67. (doi:10.1111/j.1600-0846.2008.00324.x) [PubMed]

129. Boyer G, Zahouani H, Le BA, Laquieze L
2007.
*In vivo* characterization of viscoelastic properties of human skin using dynamic micro-indentation. In *2007 Annu. Int. Conf. of the IEEE Engineering in Medicine and Biology Society, Lyon, France, 22--26 August*, vols 1–16, pp. 4584–4587. Piscataway, NJ: IEEE.

130. Goh KL, Listrat A, Béchet D
2014.
Hierarchical mechanics of connective tissues: integrating insights from nano to macroscopic studies. J. Biomed. Nanotechnol.
10, 2464–2507. (doi:10.1166/jbn.2014.1960) [PubMed]

131. Redaelli A, Vesentini S, Soncini M, Vena P, Mantero S, Montevecchi FM
2003.
Possible role of decorin glycosaminoglycans in fibril to fibril force transfer in relative mature tendons—a computational study from molecular to microstructural level. J. Biomech.
36, 1555–1569. (doi:10.1016/S0021-9290(03)00133-7) [PubMed]

132. Bischoff JE, Arruda EM, Grosh K
2004.
A rheological network model for the continuum anisotropic and viscoelastic behavior of soft tissue. Biomech. Model Mechanobiol.
3, 56–65. (doi:10.1007/s10237-004-0049-4) [PubMed]

133. Kearney SP, Khan A, Dai Z, Royston TJ
2015.
Dynamic viscoelastic models of human skin using optical elastography. Phys. Med. Biol.
60, 6975–6990. (doi:10.1088/0031-9155/60/17/6975) [PMC free article] [PubMed]

134. Lokshin O, Lanir Y
2009.
Viscoelasticity and preconditioning of rat skin under uniaxial stretch: microstructural constitutive characterization. J. Biomech. Eng.
131, 031009 (doi:10.1115/1.3049479) [PubMed]

135. Lokshin O, Lanir Y
2009.
Micro and macro rheology of planar tissues. Biomaterials
30, 3118–3127. (doi:10.1016/j.biomaterials.2009.02.039) [PubMed]

136. Fung YC.
1973.
Biorheology of soft tissues. Biorheology
10, 139–155. [PubMed]

137. Ehret A.
2011.
Generalised concepts for constitutive modelling of soft biological tissues. PhD thesis, RWTH Aachen University, pp. 1–230.

138. Bischoff J.
2006.
Reduced parameter formulation for incorporating fiber level viscoelasticity into tissue level biomechanical models. Ann. Biomed. Eng.
34, 1164–1172. (doi:10.1007/s10439-006-9124-6) [PubMed]

139. Pioletti DP, Rakotomanana LR, Benvenuti JF, Leyvraz PF
1998.
Viscoelastic constitutive law in large deformations: application to human knee ligaments and tendons. J. Biomech.
31, 753–757. (doi:10.1016/S0021-9290(98)00077-3) [PubMed]

140. Coleman BD, Noll W
1961.
Foundations of linear viscoelasticity. Rev. Mod. Phys.
3, 239–249. (doi:10.1103/RevModPhys.33.239)

141. Limbert G.
2004.
Development of an advanced computational model for the simulation of damage to human skin, pp. 1–95. Cardiff, UK: Welsh Development Agency (Technology and Innovation Division) - FIRST Numerics Ltd.

142. Limbert G, Middleton J
2004.
A transversely isotropic viscohyperelastic material: application to the modelling of biological soft connective tissues. Int. J. Solids Struct.
41, 4237–4260. (doi:10.1016/j.ijsolstr.2004.02.057)

143. Limbert G, Middleton J
2005.
An anisotropic viscohyperelastic constitutive model of the posterior cruciate ligament suitable for high loading-rate situations. In *IUTAM Symp. on Impact Biomechanics: from Fundamental Insights to Applications, Dublin, Ireland*, pp. 1--8. Dordrecht, The Netherlands: Springer.

144. Limbert G, Middleton J
2006.
A constitutive model of the posterior cruciate ligament. Med. Eng. Phys.
28, 99–113. (doi:10.1016/j.medengphy.2005.03.003) [PubMed]

145. Reese S, Govindjee S
1998.
A theory of finite viscoelasticity and numerical aspects. Int. J. Solids Struct.
35, 3455–3482. (doi:10.1016/S0020-7683(97)00217-5)

146. Lubarda VA.
2004.
Constitutive theories based on the multiplicative decomposition of deformation gradient: thermoelasticity, elastoplasticity and biomechanics. Appl. Mech. Rev.
57, 95–108. (doi:10.1115/1.1591000)

147. Vassoler JM, Reips L, Fancello EA
2012.
A variational framework for fiber-reinforced viscoelastic soft tissues. Int. J. Numer. Methods Eng.
89, 1691–1706. (doi:10.1002/nme.3308)

148. Nguyen TD, Jones RE, Boyce BL
2007.
Modeling the anisotropic finite-deformation viscoelastic behavior of soft fiber-reinforced composites. Int. J. Solids Struct.
44, 8366–8389. (doi:10.1016/j.ijsolstr.2007.06.020)

149. Nedjar B.
2007.
An anisotropic viscoelastic fibre–matrix model at finite strains: continuum formulation and computational aspects. Comput. Methods Appl. Mech. Eng.
196, 1745–1756. (doi:10.1016/j.cma.2006.09.009)

150. Flynn C, Rubin MB
2014.
An anisotropic discrete fiber model with dissipation for soft biological tissues. Mech. Mater.
68, 217–227. (doi:10.1016/j.mechmat.2013.07.009)

151. Hollenstein M, Jabareen M, Rubin MB
2013.
Modeling a smooth elastic–inelastic transition with a strongly objective numerical integrator needing no iteration. Comput. Mech.
52, 649–667. (doi:10.1007/s00466-013-0838-7)

152. Holzapfel GA, Gasser TC
2001.
A viscoelastic model for fiber-reinforced composites at finite strains: continuum basis, computational aspects and applications. Comput. Methods Appl. Mech. Eng.
190, 4379–4403. (doi:10.1016/S0045-7825(00)00323-6)

153. Pena E, Calvo B, Martinez MA, Doblare M
2007.
An anisotropic visco-hyperelastic model for ligaments at finite strains. Formulation and computational aspects. Int. J. Solids Struct.
44, 760–778. (doi:10.1016/j.ijsolstr.2006.05.018)

154. Pena E, Calvo B, Martinez MA, Doblare M
2008.
On finite-strain damage of viscoelastic-fibred materials. Application to soft biological tissues. Int. J. Numer. Methods Eng.
74, 1198–1218. (doi:10.1002/nme.2212)

155. Ehret AE, Itskov M, Weinhold GW
2009.
A micromechanically motivated model for the viscoelastic behaviour of soft biological tissues at large strains. Nuovo Cimento C
32, 73–80. (doi:10.1393/ncc/i2009-10334-7)

156. Gasser TC, Forsell C
2011.
The numerical implementation of invariant-based viscoelastic formulations at finite strains. An anisotropic model for the passive myocardium. Comput. Methods Appl. Mech. Eng.
200, 3637–3645. (doi:10.1016/j.cma.2011.08.022)

157. Simo JC.
1987.
On a fully three-dimensional finite-strain viscoelastic damage model: formulation and computational aspects. Comput. Methods Appl. Mech. Eng.
60, 153–173. (doi:10.1016/0045-7825(87)90107-1)

158. Kligman AM, Zheng P, Lavker RM
1985.
The anatomy and pathogenesis of wrinkles. Br. J. Dermatol.
113, 37–42. (doi:10.1111/j.1365-2133.1985.tb02042.x) [PubMed]

159. Veijgen NK, Masen MA, van der Heide E
2013.
Variables influencing the frictional behaviour of *in vivo* human skin. J. Mech. Behav. Biomed. Mater.
28, 448–461. (doi:10.1016/j.jmbbm.2013.02.009) [PubMed]

160. Wang S, Li M, Wu J, Kim D-H, Lu N, Su Y, Kang Z, Huang Y, Rogers JA
2012.
Mechanics of epidermal electronics. J. Appl. Mech.
79, 031022 (doi:10.1115/1.4005963)

161. Lejeune E, Javili A, Linder C
2016.
An algorithmic approach to multi-layer wrinkling. Extreme Mech. Lett.
7, 10–17. (doi:10.1016/j.eml.2016.02.008)

162. Piérard GE, Uhoda I, Piérard-Franchimont C
2004.
From skin microrelief to wrinkles. An area ripe for investigation. J. Cosm. Derm.
2, 21–28. (doi:10.1111/j.1473-2130.2003.00012.x) [PubMed]

163. Genzer J, Groenewold J
2006.
Soft matter with hard skin: from skin wrinkles to templating and material characterization. Soft Matter.
2, 310–323. (doi:10.1039/b516741h)

164. Cerda E, Mahadevan L
2003.
Geometry and physics of wrinkling. Phys. Rev. Lett.
90, 074302 (doi:10.1103/PhysRevLett.90.074302) [PubMed]

165. Biot MA.
1963.
Surface instability of rubber in compression. Appl. Sci. Res. A
12, 168–182. (doi:10.1007/BF03184638)

166. Biot MA.
1963.
Incremental elastic coefficients of an isotropic medium in finite strain. Appl. Sci. Res. A
12, 151–167.

167. Cao Y, Hutchinson JW
2012.
Wrinkling phenomena in neo-Hookean film/substrate bilayers. J. Appl. Mech.
79, 031019 (doi:10.1115/1.4005960)

168. Dervaux J, Ben Amar M
2010.
Localized growth of layered tissues. IMA J. Appl. Math.
75, 571–580. (doi:10.1093/imamat/hxq023)

169. Taylor M, Bertoldi K, Steigmann DJ
2014.
Spatial resolution of wrinkle patterns in thin elastic sheets at finite strain. J. Mech. Phys. Solids
62, 163–180. (doi:10.1016/j.jmps.2013.09.024)

170. Efimenko K, Rackaitis M, Manias E, Vaziri A, Mahadevan L, Genzer J
2005.
Nested self-similar wrinkling patterns in skins. Nat. Mater
4, 293–297. (doi:10.1038/nmat1342) [PubMed]

171. Holland MA, Li B, Feng XQ, Kuhl E
2017.
Instabilities of soft films on compliant substrates. J. Mech. Phys. Solids
98, 350–365. (doi:10.1016/j.jmps.2016.09.012)

172. Ben Amar M, Goriely A
2005.
Growth and instability in elastic tissues. J. Mech. Phys. Solids
53, 2284–2319. (doi:10.1016/j.jmps.2005.04.008)

173. Ben Amar MBA, Chatelain C, Ciarletta P
2011.
Contour instabilities in early tumor growth models. Phys. Rev. Lett.
106, 148101 (doi:10.1103/PhysRevLett.106.148101). [PubMed]

174. Dervaux J, Ben Amar M
2008.
Morphogenesis of growing soft tissues. Phys. Rev. Lett.
101, 068101 (doi:10.1103/PhysRevLett.101.068101) [PubMed]

175. Dervaux J, Ciarletta P, Ben Amar M
2009.
Morphogenesis of thin hyperelastic plates: a constitutive theory of biological growth in the Föppl-von Kármán limit. J. Mech. Phys. Solids
57, 458–471. (doi:10.1016/j.jmps.2008.11.011)

176. Goriely A, Ben Amar M
2005.
Differential growth and instability in elastic shells. Phys. Rev. Lett.
94, 198103 (doi:10.1103/PhysRevLett.94.198103) [PubMed]

177. Goriely A, Destrade M, Ben Amar M
2006.
Instabilities in elastomers and in soft tissues. Q. J. Mech. Appl. Math.
59, 615–630. (doi:10.1093/qjmam/hbl017)

178. Allen HG.
1969.
Chapter 8: Wrinkling and other forms of local instability. In Analysis and design of sandwich panels, pp. 156–189. New York, NY: Pergamon.

179. Cao Y, Hutchinson JW
2011.
From wrinkles to creases in elastomers: the instability and imperfection-sensitivity of wrinkling. Proc. R. Soc. A
468, 94–115. (doi:10.1098/rspa.2011.0384)

180. Lejeune E, Javili A, Linder C
2016.
Understanding geometric instabilities in thin films via a multi-layer model. Soft Matter
12, 806–816. (doi:10.1039/C5SM02082D) [PubMed]

181. Ciarletta P, Destrade M, Gower AL
2013.
Shear instability in skin tissue. Q. J. Mech. Appl. Math. 66, 273–288. (doi:10.1093/qimam/hbt007)

182. Gower AL.
2015.
Connecting the material parameters of soft fibre-reinforced solids with the formation of surface wrinkles. J. Eng. Math.
95, 217–229. (doi:10.1007/s10665-014-9736-z)

183. Carfagna M, Destrade M, Gower AL, Grillo A
2017.
Oblique wrinkles. Phil. Trans. R. Soc. A
375, 20160158 (doi:10.1098/rsta.2016.0158) [PubMed]

184. Yu B.
et al.
2016.
An elastic second skin. Nat. Mater.
15, 911–918. (doi:10.1038/nmat4635) [PubMed]

185. Destrade M, Ogden RW, Sgura I, Vergori L
2014.
Straightening wrinkles. J. Mech. Phys. Solids
65, 1–11. (doi:10.1016/j.jmps.2014.01.001)

186. Kuwazuru O, Miyamoto K, Yoshikawa N, Imayama S
2012.
Skin wrinkling morphology changes suddenly in the early 30s. Skin Res. Technol.
18, 495–503. (doi:10.1111/j.1600-0846.2011.00598.x) [PubMed]

187. Shiihara Y, Sato M, Hara Y, Iwai I, Yoshikawa N
2014.
Microrelief suppresses large wrinkling appearance: an in silico study. Skin Res. Technol.
21, 184–191. (doi:10.1111/srt.12175) [PubMed]

188. Wong YW, Pellegrino S
2006.
Wrinkled membranes part III: numerical simulations. J. Mech. Mater. Struct.
1, 63–95. (doi:10.2140/jomms.2006.1.63)

189. Limbert G, Kuhl E
In preparation
On skin microrelief and the emergence of expression micro-wrinkles. Soft Matter.

190. Budday S, Kuhl E, Hutchinson JW
2015.
Period-doubling and period-tripling in growing bilayered systems. Philos. Mag.
95, 3208–3224. (doi:10.1080/14786435.2015.1014443) [PMC free article] [PubMed]

191. de Souza Neto EA, Feng YT
1999.
On the determination of the path direction for arc-length methods in the presence of bifurcations and ‘snap-backs’. Comput. Methods Appl. Mech. Eng.
179, 81–89. (doi:10.1016/S0045-7825(99)00042-0)

192. Dortdivanlioglu B, Javili A, Linder C
2016.
Computational aspects of morphological instabilities using isogeometric analysis. Comput. Methods Appl. Mech. Eng.
316, 261–279. (doi:10.1016/j.cma.2016.06.028)

193. Chen L, Nguyen-Thanh N, Nguyen-Xuan H, Rabczuk T, Bordas SPA, Limbert G
2014.
Explicit finite deformation analysis of isogeometric membranes. Comput. Methods Appl. Mech. Eng.
277, 104–130. (doi:10.1016/j.cma.2014.04.015)

194. Buganza Tepole A, Kabaria H, Bletzinger K-U, Kuhl E
2015.
Isogeometric Kirchhoff–Love shell formulations for biological membranes. Comput. Methods Appl. Mech. Eng.
293, 328–347. (doi:10.1016/j.cma.2015.05.006) [PMC free article] [PubMed]

195. Javili A, Steinmann P, Kuhl E
2014.
A novel strategy to identify the critical conditions for growth-induced instabilities. J. Mech. Behav. Biomed. Mater.
29, 20–32. (doi:10.1016/j.jmbbm.2013.08.017) [PubMed]

196. Barber D.
2012.
Bayesian reasoning and machine learning. Cambridge, UK: Cambridge University Press.

197. Liu W, Röckner M
2015.
Stochastic partial differential equations: an introduction. 1st edn
Berlin, Germany: Springer.

198. Buehler MJ.
2006.
Large-scale hierarchical molecular modeling of nanostructured biological materials. J. Comput. Theor. Nanos
3, 603–623. (doi:10.1166/jctn.2006.002)

199. Rim JE, Pinsky PM, van Osdol WW
2009.
Multiscale modeling framework of transdermal drug delivery. Ann. Biomed. Eng.
37, 1217–1229. (doi:10.1007/s10439-009-9678-1) [PubMed]

200. Bancelin S, Lynch B, Bonod-Bidaud C, Ducourthial G, Psilodimitrakopoulos S, Dokládal P, Allain J-M, Schanne-Klein M-C, Ruggiero F
2015.
*Ex vivo* multiscale quantitation of skin biomechanics in wild-type and genetically-modified mice using multiphoton microscopy. Sci. Rep. UK
5, 17635 (doi:10.1038/srep17635) [PMC free article] [PubMed]

Articles from Proceedings. Mathematical, Physical, and Engineering Sciences are provided here courtesy of **The Royal Society**