Computational mechanical modeling of human skin for the simulation of reconstructive surgery procedures

Skin is the most extended organ in human body representing 16% of the total body weight with a surface extension up to 2 m2. From a mechanical viewpoint, skin can be described by an hyperelastic membrane, particularly when computational modeling for insilico testing of reconstructive surgery procedures is needed. These procedures often involves complex topological manipulations of the skin tissue in order to minimize post-operative scarring. In this paper, the simulation of reconstructive surgery procedures is described by FE membrane models developed within the framework of finite strain elasticity (an hyperelastic incompressible model for skin is adopted). An algorihm is presented to generally describe complex topologies of cutting and removing of material, while suturing is enforced by suitable multi-point constraints along wound boundaries. The archetypal reconstructive surgery of the Z-plasty is here considered, where a rotational transposition of resulting triangular flaps is involved, leading to severe stress/strain localization and displacement discontinuities. The results are discussed in terms of key deformation parameters commonly used to guide surgical decisions during reconstructive procedures. Apart from the direct applications to surgery of human skin, the computational tool proposed can be used with reference to artifical materials (like for instance polymeric hydrogels produced with advanced 3D printing technologies), whose mechanical behaviour resambles that of the natural skin tissue. © 2021 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/) Peer-review under responsibility of the scientific committee of the IGF ExCo.


Introduction
Skin is a soft tissue membrane covering the whole body, with different biological and mechanical functions, which has a noticeable self-repair capacity. After damages, reconstructive processes on skin immediately activate in order to restore the original integrity. Depending on the extension of the damage, the result may generate some defects. For instance, deep cut injuries or extended burn scars may lead to a contracted integument, which could hinder proper joint movements. In both cases, surgical correction is suggested, and modern procedures can achieve high quality results.
Among the various procedures available today, Z-plasty, invented by Horner (1837) and improved by Berger (1904), is the oldest one technique and it is highly effective with almost immediate results. It consists in a Z-shaped incision creating two opposite triangular flaps which, after being delaminated from the underneath tissues, are transposed and sutured in place. This procedure allows to elongate the central limb of about 70%, relaxing the contraction along its direction significantly.
The operation is planned and performed directly on patients drawing the best fitting Z-shaped incision. The procedure must take into account possible side effects due to final high stresses, such as blood flow reduction which can lead to flaps necrosis (Gibson and Kenedi, 1967;Larrabee et al., 1984;Raposio et al., 2000). Despite doctors experience, prognosis is not always clear and it is strongly dependent on the specific mechanical response of skin. Pioneering experiments on this problem were conducted on dogs by Furnas and Fischer (1971), performing a series of Z-plasties with different dimensions and configurations, and giving interesting relationships between wound closure tension and Z limbs angles. However, beside the ethical issues of animal testing, the reliability of such results is weak as skin properties may be different among species.
To overcome the above mentioned problems, in-silico testing of surgical procedures through the Finite Element Method (FEM) is deemed to be a strong and reliable tool. Similar studies have been conducted to investigate the mechanical behavior of simple skin wounds with different shapes of the boundary (Lott-Crumpler and Chaudhry, 2001;Flynn, 2010).
In this paper, a series of Z-plasty operations has been analyzed by the commercial FEM code ABAQUS, taking into account the highly non-linear response of skin (Gibson and Kenedi, 1967) and its initial stress state (Borges and Alexander, 1962), through the Ogden (1972) hyperelastic model. The process of flaps transposition and suturing has been simulated applying a series of implicit kinematic constraints along the boundaries of the incision. Due to the geometrical complexity of the operations, a preprocessing tool has been developed (Alberini et al., 2021) within Matlab environment, which includes an automeshing (Persson and Strang, 2004) tool, and an automatic algorithm for kinematic constraint application.

Modeling and FEM implementation
Skin is a soft tissue membrane composed by three layers, namely, epidermis, dermis, and hypodermis, each one with different mechanical characteristics and biological purposes. Epidermis is the outermost and thinnest layer, protecting the rest of the tissue from bacteria, chemicals and radiations (Martini and Nath, 2010). Just below, dermis gives the skin its mechanical strength, being able to undergo large strains without damages. This layer is mainly constituted by a matrix of collagen fibers, which uncrimples under tension, displaying a highly non linear response (Gibson and Kenedi, 1967;Brown, 1973;Yang et al., 2015). Finally, hypodermis connects the upper layers with the underneath tissues, such as muscles and bones, allowing large displacements under low tensions (Oxlund et al., 1988;Groves et al., 2012). The overall mechanical behavior of skin is characterized by the dermis and epidermis, which can be treated as a unique membrane.
In the following, finite strain continuum mechanics is briefly introduced as the fundamental tool to analyze bodies undergoing large displacements. Then, the implementation of skin corrective surgeries into FE framework is discussed.

Continuum mechanics of membranes
Following the notation of Holzapfel (2000), let B be a material body in R 3 moving from the space region Ω at time t 0 = 0 to ω at time t. A material point belonging to B moves from the position X to x according to the biunivocal function x = χ(X, t) = X + u(X, t), where u is the displacement field. During the motion χ the body may undergo deformations which can be described by the fundamental strain measure F = ∂x ∂X , also known as deformation gradient. The polar decomposition leads to F = RU, where R is an orthogonal tensor representing a pure rotation, and U is a pure stretch tensor, which admits the spectral decomposition in terms of principal stretches λ a and principal basis vectorsN a , U = λ aNa ⊗N a . The volume ratio is given by the determinant of F, that is J = det(F) = λ 1 λ 2 λ 3 . As a further strain measure, the left Cauchy-Green is defined as C = F T F = U 2 , and admits the spectral decomposition as well, C = λ 2 aNa ⊗N a . Assuming skin as an hyperelastic, isotropic and homogeneous membrane, the constitutive relationship can be found from a scalar valued strain-energy density function ψ = ψ(F) = ψ(C) in terms of second Piola-Kirchhoff stress S, with The Cauchy true stress tensor σ is obtained by means of the push forward operation on S, σ = J −1 FSF T . Both stress measures can be also written as S = S aNa ⊗N a , σ = σ ana ⊗n a , withn a = RN a , where the principal stresses are expressed as Skin behavior is here described by the Ogden (1972) model, assuming the material as incompressible. The strainenergy density function is defined as where ψ p = (J − 1) has been introduced to account for the incompressibility constrains, in which p is an unknown hydrostatic pressure. This model is commonly used to describe the actual skin behavior due to its simplicity and accuracy with just one series (N = 1) (Ogden et al., 2004;Shergold and Fleck, 2004;Groves et al., 2012). In plane stress conditions the pressure p can be obtained explicitly by substituting Eq. (3) in Eq.

FE modeling
In silico testing of skin corrective surgeries are performed in three steps. First, the geometry of the cut is drawn directly on the skin depending on the extension of the defect, and the amount of neighboring tissue. Then, the incision is performed and skin is delaminated from subcutaneous tissues in a circular region underneath the cut, allowing flaps to move and the tissue to rearrange under the induced deformations. Finally, flaps are transposed and the wound is sutured. This procedure is rather complex and a precise simulation must include as much physical aspects as possible.
The domain considered includes the circular region of delaminated skin, which behaves as a free membrane, while the surrounding skin, still connected to the underneath tissues, has been replaced with proper constraints around the circular boundary ∂Ω E . The domain is taken sufficiently large, so that the operation has negligible effects on the external boundary, which can therefore be considered fixed.
The internal boundary ∂Ω I , represented by the edges of the cut, has unknown boundary conditions. Skin flaps undergo large displacements, and the consequent geometrical non-linearity makes the final configuration of the edges, as well as the relative traction forces, an unknown of the problem. Thus, in order to simulate the process of wound closure, the constraint must be given as an implicit kinematic equation, which progressively reduce to a vanishing distance the material points along two approaching edges ∂Ω I− and ∂Ω I+ of the inner boundary, (∂Ω I− ∪ ∂Ω I+ = ∂Ω I , ∂Ω I− ∩ ∂Ω I+ = ∅). Therefore, for every couple of points X − ⊂ ∂Ω I− , X + ⊂ ∂Ω I+ , in the reference configuration Ω, the final positions x − ⊂ ∂ω I− and x + ⊂ ∂ω I+ , respectively, must coincide ( Figure 1). The implicit kinematic equation Φ can be written as and it is applied continuously along the edges ∂Ω I− , ∂Ω I+ . Accordingly, in order to preserve the integrity of the material, every coupling (X − , X + ) must be defined by a linear biunivocal function Γ : ∂Ω I+ → ∂Ω I− , such that dX − = c dX + , where c is the ratio between the total lengths of ∂Ω I− and ∂Ω I+ . Due to the discrete nature of FEM framework, Eq. (5) is applied on the nodes of the internal boundary, which must be placed according to the function Γ. Moreover, the boundary of the wound may include several geometrical discontinuities that require smaller elements for the discretization in their influence area. Therefore, nodes along the inner boundary must be placed also according to the refinement rule for the elements mesh. In order to respect both requirements, a code for automatic generation of FEM models has been developed in Matlab environment, and has been integrated with the open source automeshing code DistMesh (Persson and Strang, 2004). The code uses 3-node plane stress elements, and it is amenable to handle efficiently multiple geometrical parameters to generate the optimal discretization. Furthermore, the code is capable to generate models for fairly different geometries of the internal boundary (Alberini et al., 2021).

Z-plasty optimization
The classical Z-plasty is made by three incision of equal length, in which the lateral limbs are slanting 60 • with respect to the central one. Surgeons came up with this configuration during years of practice, but no mechanical analysis demonstrated that this is the best configuration possible. Moreover, depending on the angles chosen for the lateral limbs, even unsymmetrically, different final configuration can be achieved. Therefore, in order to identify an optimal configuration of the operation, a series of Z-plasties is analyzed within the FEM code ABAQUS, adopting different combination of slanting angles of the lateral limbs. The round region of delaminated skin is characterized by a radius R = 200 mm, and a thickness t = 2 mm. The central incision limb, and the lateral ones as well, is made with a length l = 5 mm, leaving a small gap s = 0.05 mm between the edges of the cut, for numerical reasons (Figure 2a). The domain has been automatically discretized using an element dimension of h max = l/5, which reduces to h min = l/150 near the geometrical discontinuities.
Generally, the angles used in practical surgery for a symmetric Z-plasty are α = β (Figure 2a) ranging from 30 • to 60 • . For the optimization, this interval has been extended to 20 • ≤ α, β ≤ 90 • , taking combinations of α and β with step of 10 • . For symmetric reasons, only 36 combinations of 64 have been analyzed to save computational effort (Figure 2b).
The material properties used are taken from an in-vivo measurement which provides µ 1 = 110 Pa and α 1 = 26 (Mahmud et al., 2013). To take into account the contracted skin, an anisotropic tension field is uniformly applied to the domain considering the principal stretches λ 0,1 = 1.05 and λ 0,2 = 1.1 applied along the vertical and horizontal directions, respectively, (Gambarotta et al., 2005). Accordingly, the initial stress tensor σ 0 has been computed using the principal stresses featuring in Eq. (2).
The performance of the operation is evaluated through the elongation ratio, defined by the ratio λ AB between the distance of the two cut extremes after (points A and B , respectively) and before (points A and B, respectively) the operation, namely λ AB = A B AB . Results of Figure 4(a) show that λ AB always increases with increasing angles α and β. Despite this interesting property, higher values for α and β are not necessarily an appropriate choice. As shown in Figure 3(a) and (b), wider angles generate increased deformations, and this could lead to the failure of skin flaps for necrosis, due to the reduced blood flow (Gibson and Kenedi, 1967;Stell, 1980). The same problem may occur for low angles α and β, since blood flow also depends on the width of the base of the flap, along which it remains connected to the rest of the skin. Thus, in order to assess the elongation performances of each combination of angles, an efficiency parameter η AB is ruled out by penalizing the ratio λ AB with the highest effective strain ε e f f = 2 3 (ε 2 1 + ε 2 2 ) within the domain, and the sine of the lower angle between α and β (note that the area of each triangular flap turns out to be proportional to this sine function), that is η AB = λ AB ε e f f,max sin(min(α, β)).
The results, reported in Figure 4(b), show that for every fixed value of β the peek efficiency is achieved with α = β. This feature is significantly helpful as the choice of the angles can be reduced to a single variable. In particular, the combination α = β = 60 • gives the highest elongation efficiency, and this justify its wide use in the common surgical practice.

Conclusion
Z-plasty is an old surgical procedure performed to elongate contracted skin. This technique is still used today thanks to its immediate relaxation effects, and speed of execution. The geometry of the cut is simple and can be performed with different angles of the lateral Z-shaped limbs. The classical choice with equal angles of 60 • has become a standard, but this came up from years of experience instead of a mechanical-based justification. In the present work, a geometrical optimization has been addressed, simulating the process of suturing on a series of Z-plasties with different combinations of α and β angles, between 20 • and 90 • , by means of the FE models. For this purpose, a Matlab code has been developed in order to efficiently generate the models. The code is capable to discretize complex domains with localized mesh refinement, and automatically apply the suturing constraint along the inner boundary of the wound. The analysis shows that the elongation ratio λ AB , the key parameter to evaluate the elongation performance of the operation, always increases with increasing angles α and β, but the optimal configurations, which minimize the risk of necrosis, is achieved with α = β. In particular the combination α = β = 60 • reaches the highest elongation efficiency, confirming it as the best choice in common surgical practice. Using the numerical model recently proposed by the first three authors (Alberini et al., 2021), further investigations are under way by considering the simulation of Z-plasty in the skin membrane undergoing wrinkling deformations (e.g. see the problem of wrinkling in thin metal foils in Bolzon et al. (2017a,b)).