Three-dimensional simulation of scalp soft tissue expansion using finite element method.

Scalp soft tissue expansion is one of the key medical techniques to generate new skin tissue for correcting various abnormalities and defects of skin in plastic surgery. Therefore, it is very important to work out the appropriate approach to evaluate the increase of expanded scalp area and to predict the shape, size, number, and placement of the expander. A novel method using finite element model is proposed to solve large deformation of scalp expansion in this paper. And the procedure to implement the scalp tissue expansion with finite element method is also described in detail. The three-dimensional simulation results show that the proposed method is effective, and the analysis of simulation experiment shows that the volume and area of the expansion scalp can be accurately calculated and the quantity, location, and size of the expander can also be predicted successfully with the proposed model.

1. Introduction

3D imaging has been widely applied in the current plastic and reconstructive surgery. Noninvasive computed tomography (CT), magnetic resonance imaging (MRI), and three-dimensional laser scanning are increasingly used to generate tissue structural views for 3D anatomical model [1] to assess facial growth [2, 3], facial expressions [4, 5], facial asymmetry of cleft lip and palate patients [6], and facial reconstruction [7]. This kind of techniques is also expected to model the reconstruction accurately and to make plastic surgery planning as a truly interactive procedure. Many literatures have also shown that various 3D skeleton models have been made in maxillofacial surgery. During a planning process, the simulation of osteotomy can be checked by the computer-aided design (CAD), where the actual lines of osteotomy can be clearly identified with the help of CAD.

However, the assessment of soft tissue in the plastic surgery is more difficult than bones due to the biomechanical properties of soft tissue, such as nonhomogeneous, quasi-incompressible, and nonlinear plastic-viscoelastic material properties. Therefore, computer-aided based soft tissue expansion technique is introduced into three-dimensional imaging to assess the structure of soft tissue. Skin soft tissue expansion technique (as shown in Figure 1) is also one of the key medical techniques in surgery planning, which is used to generate new skin flaps for correcting various skin abnormalities and defects. Therefore, it is widely applied in many fields such as plastic and reconstructive surgery, cosmetic surgery, and reparative and reconstructive surgery. It is also very useful and practical in the treatment of soft tissue defects in the head and facial area. Meanwhile, it is critical for a successful surgery planning to effectively select the shape, size, number, and buried location of the expander according to the practical state of the defect skin. Because inadequate expanded flap size may result in failure of covering the skin defect without enough tension and overexpanded flap size means wasting tissue, Ji et al. [7] obtained the data of the scar excision in a child with burned injuries and the expanded cervicofacial flap by using a 3D digital scanner. The proposed result shows that the scar area planned for excision matches the area of the face and anterior neck with tissue expansion well. But the result is of theoretical value, and the influence of biomechanical properties of soft tissue has not been considered during the whole procedure, such as wound retraction and flap shrinkage.

Figure 1

Skin soft tissue expansion surgery.

So far, the accuracy of operation still mainly relies on the surgeon's clinical experience accumulated from long-term practices. It is quite subjective and unstable. To solve the problems well, simulation of tissue deformation during its expansion process can provide an additional modality to improve the mission success rate in tissue expansion surgery, so it will be discussed in detail in this paper. The main purpose of this study is to find out the appropriate approach to obtain more accurate data considering wound retraction and flap shrinkage after removing the tissue expanders in surgery, which is based on the images generated from CT scans. The CAD-based mathematical model is then to be constructed to simulate the process of the surgery, which consisted of selecting the proper expanders, removing the tissue expanders, and covering the soft tissue defect with the expanded flaps.

In Figure 1(a), there is a practical example of two liquid expanders implanted into the patient's head; what is shown in Figure 1(b) is the healed wound after removing expanders.

The rest of this paper is organized as follows. Section 2 provides the methods employed to simulate scalp soft tissue expansion. Section 3 focuses on the experiments and the discussion of the method and experiment results. Section 4 presents a conclusion and the future work.

2. Simulation of Scalp Soft Tissue Expansion with Finite Element Method

Geometric model is the foundation of deformable simulation. And three-dimensional model is widely used in the simulation of medical tissue deformation. The scalp of the model usually has a certain thickness. According to the structural characteristics of the scalp, the thickness is set to 3.4 mm [8]. Besides, the model of tetrahedral grids [9] is obtained by free division method [10] using software Abaqus. The result is shown in Figure 2.

Figure 2

Scalp model represented by tetrahedral grids.

Scalp expansion is a kind of deformation in large range and nonlinear problem that contains geometric and material issues [11]. In this paper, a dynamic finite element model is introduced to solve the nonlinearity of the problem [12]. The finite element theory based geometric equations and equilibrium conditions for small deformation are no longer suitable due to the geometric nonlinearity. In general, there are two main ways to describe large deformation problem, material description (Lagrange description), and space description (Euler description) [13]. At the initial time t 0 = 0, the coordinate of number of i points is X (i = 1,2, 3) and becomes x after motion at any time. The motion of the point can be represented by the following equation:

The above equation is the Lagrange description used in this paper that examines the movement and deformation using the motion of a specific point. The configuration of the object before deformation is already known and is also as a reference model. The configuration of the object after deformation is computed by finite element method. Currently, the pair of second class Piola-kirchhoff stress tensor and Lagrant-Green strain tensor is widely used to express the energy item using finite element method to solve nonlinear problem. This pair tensor takes initial configuration as reference configuration. The solving steps of scalp soft tissue expansion with finite element method are as follows.

2.1. Discrete and Equal Parameter Unit Interpolation

For the initial configuration, its geometry of internal units is interpolated with the coordinates of the unit points. Besides, unit displacement is also obtained by the same interpolation function as follows: where X is the coordinate of point K before deformation in the direction i, u is the displacement of the point k in the direction i, and m is the number of the unit points. They can be represented as vectors, X = NX and U = Na , N = [N 1 I, N 2 I,…, N I], where I is a 3 × 3 unit matrix and N is a 3 × 3m shape function matrix. Vector X is the coordinate vector of the initial unit point and a is the displacement vector of the unit point. Both of them are 3m in dimension.

2.2. Derivation of Strain Matrix B

Strain matrix is represented by the strain tensor of Green in this paper. At first, strain tensor of Green E is represented by two parts, linear part E and nonlinear part E . So, there are E = E + E and B = B + B . B is the transformational matrix between E and a ; that is, E = B a . B is the transformational matrix between E and a , and E = A∗θ/2. The detailed solving process of B and B is described as follows: where θ is a displacement gradient matrix, ∂U/∂X is a 3 × 1 matrix, and 0 represents a 3 × 1 zero matrix. And then B can be calculated as in the following equation. That is B = AG, where G is obtained by

2.3. System Balance Equation

For a triangle composed of three points, second class Piola-kirchhoff stress tensor can be described as . Physical force load and surface force load are, respectively, represented as , . System balance equation using virtual work equation is defined as where matrix c combines displacement vector a of element point with the total displacement vector a of finite element system via the equation of δa = c δa. The balance equation of the whole system can be defined as (7) and obtained by adding the total unit balance equations given as (6)

To solve (7), scalp constitutive model is needed to be introduced and defined as follows [14]: where w represents strain energy density, I 1 represents the first strain invariant, λ 1, λ 2, and λ 3 are the stretch variables, respectively, in directions x, y, z, and a, and b represents parameter s of the expander. Because scalp constitutive model is a superelastic material model and W = W(I 1, I 2, I 3), the scalp is defined as superelastic material in this paper. Considering that external force has nothing to do with the deformation path of the material, I 1 = E , I 2 = [(E )2 − E × E ]/2, and E = E = E . For incompressible material, I 3 = 1. The constitutive equation with rate form can be obtained by derivation of deformation variables

where S ′ is the derivation of time by second class Piola-kirchhoff stress tensor, E ′ is the derivation of time by Green strain tensor, and D is tangent modulus tensor. For a three-dimensional model, object region is V 0 and boundary is A 0. The equivalent nodal load R of FEM system can be gotten by combining D with balance equation as described in the following equation:

2.4. System Tangent Stiffness Matrix K

In general, it is difficult to obtain accurate stiffness matrix for nonlinear material, so stiffness matrix can be replaced by tangent stiffness matrix which is the curve tangent of stress-strain defined as follows: where K = K + K and K = K + K . K is a tangent stiffness matrix associated with constitutive matrix, K is a general small displacement stiffness matrix, K only including a linear or quadratic term is caused by a large displacement, and K is a tangent stiffness matrix by stress D in formula (12) contacts unit stress tensor and strain tensor, representing dS = D dE. K can be calculated by (13) where is the transformation matrix of displacement gradient vector and unit nodal displacement vector.

2.5. Solving Equations by Newton Method

Step 1 .

Solving linear elastic problem K a − R = 0, first approximate solution a 1 is obtained. a is the total displacement vector of FEM system.

Step 2 .

Compute matrix A from a 1 based on (4). Displacement gradient vector θ = Ga . Compute E through δAθ = A(δθ) based on the definitions of A and θ. Then, compute S 1 and B 1 based on scalp constitutive model D and S = DE

Step 3 .

Get unbalanced force φ via the left part of (9a) minus the right part of (9b). Get φ 1 through S 1 and B 1.

Step 4 .

Solve tangent stiffness matrix K 1 with (11) to (13).

Step 5 .

Compute the correction of displacement via Δa 1 = −(K 1)−1 φ 1. Get second approximate solution a 2 = a 1 + Δa 1.

Step 6 .

Iteration on a 2 from Step 2 to Step 5 till φ is small enough.

2.6. Conditions of Loading and Boundary

The expander usually is set near the defect area of patient. The stress of the scalp packaging expander is caused by liquid in expander and the direction of stress is along the normal expander's surface. Therefore, the stress of scalp contact parts is the same as the stress of expander, which is set as 30 N/cm2 in this paper. Besides, the expander is set as constant regular ellipsoid whose triaxial proportion is 1.53 : 1 : 1. Suppose the volume of the virtual ellipsoid is zero and the center point of the virtual ellipsoid is at the center of specified area at the beginning. During the scalp deformation, the ellipsoid expands and its volume increases. The distance between the center point and scalp surface is computed to decide the virtual ellipsoid and scalp surface whether the virtual ellipsoid and the scalp surface is contacted or not.

3. Experiments Result and Analysis

3.1. Experimental Result

The simulation is based on the large deformation of scalp expansion and conditions of loading and boundary above. Eight steps are set in software Abaqus and the time of each step is 0.04 s. Initial state and main process results of the 8 steps are shown in Figure 3.

Figure 3

Results of scalp expansion simulation.

3.2. Stress-Strain Analysis

The color in Figure 3 represents the average stress distribution: red means the maximum and blue means minimum. It is obvious that the stress is relatively large in confined areas and the stress decreases from the top of scalp to the bottom of scalp. Five tetrahedrons are selected and the average strain-time relationship of four points of each tetrahedron is shown in Figure 4.

Figure 4

Stress-strain analysis diagram.

As shown in Figure 4, the stress of tetrahedron 1, 3, or 4 at the edge of area is larger than that of tetrahedron 5 from 0 s to 0.016 s. When the expander is implanted into head, the downward stress from scalp is large, the strain of top tetrahedron 5 is restricted, and tetrahedrons 1, 2, 3, and 4 around expand outward. Then, when the expansion strain of the edge regional is growing and the volume of expander is increasing, the expansion to the top of scalp becomes obvious. Therefore, the strain of tetrahedron 5 increases faster than that of tetrahedrons 1, 2, 3, and 4. With the restrictions of the fixed scalp around, the strain of tetrahedrons 1 and 2 is quite small. The strain of tetrahedron 3 is relatively large for the corner position. The strain of tetrahedron 4 is neither large nor small for the middle layer position.

3.3. Area Analysis

Before the expansion, coordinates of scalp surface model are obtained at random. Then, these points are connected in delta in real time using Delaunay triangulation method [15] and the area of each triangle is computed. After accumulation, the area of defect is 4266.04 mm2 at last. The process is shown in Figures 5 and 6.

Figure 5

Random sampling inside and outside the model.

Figure 6

Getting area by interactive sampling in real time.

With the process of calculation result by software Abaqus, the area of scalp after deformation is shown in Table 1.

Table 1

The area of scalp after deformation.

Step	S (mm2)	ΔS (mm2)	∑ΔS (mm2)
0	117733.11	0	0
1	118504.79	771.68	771.68
2	119319.77	814.98	1586.66
3	120392.57	1072.8	2659.46
4	122278.13	1855.56	4545.02
5	123252.50	974.37	5519.39
6	123638.54	386.04	5905.43
7	124240.53	601.99	6507.42
8	124522.65	282.12	6798.54

After removing the expander, skin shrinkage phenomenon will occur. The part of shrinkage will be 30% based on medical experience, so about 1.43 times more new skin is needed. After the end of the experiment, the area is increased to 6789.54 mm2 which is 1.59 times of the defective area (4266.04 mm2). Considering the utilizable efficiency of new skin, the surgery can be satisfied.

3.4. Volumetric Analysis

In order to estimate the volume needed in the operation of scalp tissue expansion, each step of deformation is processed by software and the volume is shown in Table 2.

Table 2

The change of volume.

Step	V (mm3)	ΔV (mm3)	∑ΔV (mm3)
0	3107136.50	0	0
1	3166015.68	58879.18	58879.18
2	3220374.85	54359.17	113238.35
3	3289162.78	68787.93	182026.28
4	3404764.16	115601.38	297627.66
5	3467026.40	62262.24	359889.9
6	3491925.98	24899.58	384789.48
7	3530934.93	39008.95	423798.43
8	3549148.62	18213.69	442012.12

As Table 2 shows, the volume of head increases totally 442012.12 mm3 (442.01212 mL). It is almost equal to the volume of expander indeed. So, a 450 mL ellipsoid expander is needed to be implanted. The results show that the large deformation method proposed in this paper is effective.

4. Conclusion

Tissue expansion is a good option for covering the soft tissue defect. Successful reconstruction is depended on the precise judgement on the amount of tissue provided by expansion to cover the defect. The three-dimensional anatomy alters the situation that elasticity and contractility of the expanded flap make it extremely difficult to accurately predict the proper size implants and the size of skin flaps required to cover the defects.

Based on small deformation and linear elastic problems with finite element method, a novel solution to large deformation of scalp expansion is put forward in this paper. Then, the concrete steps to implement the scalp tissue expansion process with finite element method are also given in detail. The scalp tissue is simulated as a shell with certain thickness and is split into tetrahedral meshes. The deformation results prove that the solution for large deformation is effective. Then, the stress during the deformation process is also analyzed, and the volume and the area of the scalp are accurately calculated. With the proposed method, the quantity, placement location, quantity, and size are predicted successfully.

However, the proposed model used to approximate flap shrinkage is rough, and it does not take into consideration other variables such as thickness of flaps and the length of expansion time which will make the rate of shrinkage different. Further studies are required to make the model more accurate.