A FINITE ELEMENT STUDY OF HUMAN THIGH AREA IN SEATED POSTURE FOR PRESSURE ULCER PREDICTION AND PREVENTION By Sheng Chen A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirement s for the degree of Mechanical Engineering Doctor of Philosophy 2019 ABSTRACT A FINI TE ELEMENT STUDY OF HUMAN THIGH AREA IN SEATED POSTURE FOR PRESSURE ULCER PREDICTION AND PREVENTIO N By Sheng Chen Pressure ulcers (PUs), also known as pressure sores, are localized damage to the skin and underlying tissues, usually occurring over a bony prominence and caused by sitting or lying in one position for long time. PUs are a detriment to the well - being of people who lose their mobility either permanently or temporarily, and high morbidity and mortality are associated with P U s . Although the initiating mechanism of PUs is still unclear, it is commonly accepted that internal normal and shear stresses, due to the presence of unrelieved external loads, play a central role in the formation and development of these wounds . Despite the significance of inter nal stresses in PUs formation, interfacial pressures, which are a surface measure of stress, are the indicators commonly used to develop practices and protocols to minimize loading on the soft tissue . However, no direct correlation exists between interfacial pressure and internal stresses of soft tissue . Therefore, tools and methods that can show strains as a response to external loading are needed . The ability of finite element (FE) models to accurately represent the anatomical structure of the leg and buttocks area and to estimate the localized stress/strain fields within highly deformable media, makes them powerful tools to investigate soft tissue response to external loadings . Despite the significant advancement previous studies have achieved, there are still important aspects in human thigh - b uttock soft tissue modeling area that need to be improved. Two challenges are identified in this dissertation: 1) Microstructurally motivated skin modeling for an individual skin layer in finite element model. 2) Parameters estimation associated with large deformations. To address the first challenge , a micro structurally based constitutive model is proposed to describe the mechanical behavior of skin. The constitutive model incorporated the distribution of collagen fiber bundle orientations and relative col lagen content measured from histology, and shows good agreement with the tensile test data. To address the second challenge , an optimization procedure that is able to match nonlinear behaviors between FE simulation and in vivo experimental data is developed. The difference between 3D and semi - 3D model is quantified, and the accuracy of four commonly used constitutive model representing soft tissue nonlinear mechanical behavior is compared. Finally , a thigh FE model that has det ailed anatomical representation of different soft tissue types, i.e., skin, fat, and muscle, is developed. The subject - specific in vivo experimental data are used to inform the optimization procedure to obtain best - fit constitutive parameters for different soft tissue types. The research in this dissertation provides an approach to describe the in vivo mechanical behavior of soft tissues in thigh - buttock area accurately through FE modeling. The constitutive parameters informed by in vivo data in this disser tation are valuable to facilitate future FE modeling studies to achieve accurate internal stress/strain distribution of soft tissues in thigh - buttock area. iv To my parents and my wife. v ACKNOWLEDGEMENTS I would like to thank my advisor, Dr. Sara Roccabianca, who has been a great mentor and a dear friend to me while pursuing my PhD. Her patient guidance, rich knowledge, positive energy and hard - working spirit motivated me along this journ ey. I would also like to thank my dissertation committee, Dr. Tamara Reid Bush, Dr. Seungik Baek, and Dr. Weiyi Lu for their guidance and insights on my research. Last but not least, I would like to thank my parents and wife, for their endless love and sup port. vi TABLE OF CONTENT S LIST OF TABLES ................................ ................................ ................................ ..................... viii LIST OF FIGURES ................................ ................................ ................................ ..................... ix INTRODUCTION AN D LITERATURE REVIEW ................................ .......................... 1 1.1. Overview ................................ ................................ ................................ ......................... 1 1.2. Skin mechanical model ................................ ................................ ................................ .. 2 1.3. FE models of leg and buttock area ................................ ................................ ............. 14 1.4. Objectives of this dissertation ................................ ................................ ..................... 23 CHARACTERIZATION OF THE MECHANICAL PROPERTIES OF SKIN USING RAT AND A MICROSTRUCTRUALLY BASED CONSTITUVIE MODEL .............. 24 2.1. Overview ................................ ................................ ................................ ....................... 24 2.2. Methods ................................ ................................ ................................ ......................... 24 2.2.1. Mechanical testing ................................ ................................ ................................ 24 2.2.2. Histology and collagen network analysis ................................ ............................ 29 2.2.3. Constitutive model ................................ ................................ ................................ 31 2.2.4. Parameter estimation ................................ ................................ ............................ 34 2.3. R esults ................................ ................................ ................................ ........................... 35 2.3.1. Storage condition effect ................................ ................................ ........................ 35 2.3.2. Location, orientation and sex effect ................................ ................................ .... 37 2.3.3. Fibers orientation, straightness and relative collagen content ......................... 38 2.3.4. Constitutive model parameters ................................ ................................ ............ 40 2.4. Discussion and conclusions ................................ ................................ .......................... 42 A SUBJECT SPECIFIC MODEL THAT ADDRESSES NONLINEAR MECHANICAL BEHAVIOR OF HUMAN THIGH INFORMED BY IN VIVO EXPERIMENTAL DATA ................................ ................................ ................................ ................................ .... 48 3.1. Overview ................................ ................................ ................................ ....................... 48 3.2. Methods ................................ ................................ ................................ ......................... 48 3.2.1. Mechanical testing ................................ ................................ ................................ 48 3.2.2. Geometry generation ................................ ................................ ............................ 49 3.2.3. Finite element model ................................ ................................ ............................. 49 3.2.4. Constitutive models and parameter optimization ................................ .............. 52 3.2.5. Statistical analysis ................................ ................................ ................................ . 55 3.3. Results ................................ ................................ ................................ ........................... 55 3.3.1. Effects of the choice of constitutive model and of modeling region dimension 55 3.3.2. Effect of sex and location ................................ ................................ ...................... 58 3.4. Discussion and conclusions ................................ ................................ .......................... 60 A NONLINEAR HUMAN THIGH MODEL WITH HIGH ANATOMICAL AND MECHANICAL FIDELITY ................................ ................................ ............................... 69 vii 4.1. Overview ................................ ................................ ................................ ....................... 69 4.2. Methods ................................ ................................ ................................ ......................... 69 4.2.1. Mechanical behavior of soft tissues in vivo ................................ ......................... 69 4.2.2. Geometry generation ................................ ................................ ............................ 70 4.2.3. Finite element model ................................ ................................ ............................. 72 4.2.4. Constitutive models and parameter estimation ................................ ................. 74 4.2.5. Parameter estimation ................................ ................................ ............................ 76 4.3. Results ................................ ................................ ................................ ........................... 79 4.3.1. Best - fit constitutive parameters for fat and muscle ................................ ........... 79 4.3.2. Skin constitutive parameters in vivo ................................ ................................ ... 80 4.3.3. Effect of in vitro and in vivo skin constitutive parameters on thigh mechanical behavior ................................ ................................ ................................ ............................... 80 4.4. Discussion and conclusions ................................ ................................ .......................... 82 CONCLUSIONS AND FUTURE WORK ................................ ................................ ......... 86 BIBLIOGRAPHY ................................ ................................ ................................ ....................... 89 viii LIST OF TABLES Table 1 . The NRMSD between each storage protocol and the control group for cyclic loading stress - stretch curves and rupture loading stress - stretch curves. ................................ .................... 36 Table 2 . Means and standard deviations of best - fit constitutive parameters , , , and , grouped based on location, orientation, and sex. ................................ ................................ .... 41 Table 3 . Initial guess constitutive parameters for each constitutive model. ......................... 55 Table 4 . Best - fit constitutive parameters at proximal, middle, and distal thigh location for each constitutive model (i.e., neo - Hookean, Money - Rivlin, Ogden, Fung orthotropic). Experimental data collected for representative subject F6 have been used to in form the 3D and semi - 3D FE model. The percentage difference of best - fit constitutive parameters between the two indicate semi - 3D model have higher parameter values, and negative values indicate 3D model have higher parameter values). ................................ ................................ ................................ ..... 57 Table 5 . P - - test comparing sex difference on Ogden and Fung orthotropic constitutive parameters at different thigh locations. The gray areas highlight the statistical significance, i.e., p < 0.05. ................................ ................................ ................................ ............ 59 Table 6 . P - values of Holm - Sidak post - hoc pairwise tests comparing model parameter difference between different thigh locations, in male and female subjects group, for both Ogden and Fung constitutive models, respectiv ely. The gray areas highlight the statistical significance, a lighter gray corresponds to statistical significance ( p < 0.05), and a darker gray corresponds to strong statistical significance ( p < 0.001). ................................ ................................ .................... 59 Table 7 . Best - fit constitutive parameters for fat and muscle tissues that describe the in vivo indentation test data collected in quadruped and prone postures. ................................ ................. 80 Table 8 . The constitutive parameters of the thigh FE model with in vitro skin constitutive model and the thigh FE model with in v ivo skin constitutive model. The NRMSD values show the normalized root - mean - square deviations between the FE simulated force - deflection curve and the in vivo indentation data collected in quadruped posture. ................................ .............................. 82 ix LIST OF FIGURES Figure 1 . ( a ) Schematic of skin structure. ( b ) A typical load - strain graph obtained in an in vitro uniaxial test of human abdominal skin, three phases of extension are indicated in the graph. ( c - f ) Collagen network, parallel section through mid - dermis, under the scanning electron microscope at different strain levels. ( c ) Unstressed skin, collagen fibe rs appear wavy and multidirectional. ( d ) Strain level corresponding to the end of phase 1, fibers have reoriented towards load axis. ( e ) Strain level corresponding to the end of phase 2, fibers appear aligned in the direction of the load and almost straigh t. ( f ) Strain level corresponding to mid - phase 3, fibers appear aligned, straight, and compacted. Pictures ( b - f ) from [31] , numbering changed from original. ................................ ................................ ................................ ................................ ........... 4 Figure 2 . Computer - generated images ( a, c ) and the power plots of Fourier analysis ( b, d ) previously published in [35] . Specifically, ( a ) and ( c ) show two examples of the computer - generated images, which were used to validate the accuracy of Fourier analysis generated power plots in representing fiber orientation preference. The images consist of ellipses with a width/length ratio of 1:10. ( b ) and ( d ) are the power plot generated by Fourier analysis, of ( a ) and ( c ), respectivel y. The network represented in ( b ) corresponds to an orientation index of 1, and the network represented in ( d ) corresponds to an orientation index of 0. ( e - h ) Laser scatter images and first - order power plots previously published in [36] . ( e ) The laser scat ter image of lesional sclerodermic skin (LS) produces an elongated scatter plot, indicating a higher orientation of collagen bundles when compared to control skin (CS), shown in ( f ), which produces a circular scatter plot. ( g - h ) First - order power plots fro m the fast Fourier transform of the images shown in ( e - f ). ( g ) LS creates a power plot with a preferential distribution with two maxima, and CS, shown in ( h ) is characterized by a power plot with a circular pattern. Numbering changed from original studies. ................................ ................................ ................................ ................................ ............ 5 Figure 3 . ( a ) Geometry of a model collagen fiber as described in [46] . ( b - c ) Schematics of the collagen and elastin network struc ture in flat tissue from [42], [43] . ( b ) Tissues with high density of crosslinks and elastin induced collagen undulation. ( c ) Tissues with low density of crosslinks and inherent collagen undulation. Thick lines represent collagen; thin lines represent elastin, and dotted lines represent overall collagen fiber direction. ................................ ................ 8 Figure 4 . Left : Maximum intensity projection for a 60 transverse section showing collagen fibers in bright red. Right : A mixture of two von Mises distribution fitted to fiber orientations presented in the form of histogram, under various symmetry assumptions. Pictures are from [37] . ................................ ................................ ................................ ................................ ........ 9 Figure 5 . Images output from automated algorithm. ( a ) Original histology slide. Scale bar is 1 mm. ( b ) Binarized image after automated thresholding. ( c ) Bi narized image after erosion step. ( d ) All identified collagen bundles outlined in green. ( e ) Remaining bundles which meet area and eccentricity criteria. ( f ) Best fit ellipse about each fiber that meets the specified criteria. ( g ) Histogram of collagen orientations. The two distinct peaks correspond to the preferred orientation of the two fiber families. All images are from [38] . ................................ ................................ ..... 11 x Figure 6. a posterior trunk and scalp. ( b c leg. ( d ) Wrinkle lines on the poste rior trunk and scalp. ( e ) Wrinkle lines on the anterior leg. ( f ) Wrinkle lines on the posterior leg. ................................ ................................ ................................ 13 Figure 7 . FE model of the p elvis and thighs published in [21] . ( a ) Bony structures, ( b ) bony structures and human soft tissues, and ( c ) complete model (i.e., bony structures, human tissues, and skin). ( d - e ) Pressure distributions (Pa) predicted by the FE pelvis model (shown in d ) and measure experimentally for a human male (shown in e ) for a rigid surface boundary condition. ( f - g ) Pressure distributions (Pa) predicted by the FE pelvis model (shown in f ) and measure experimentally for a human male (shown in g ) for a soft cushion b oundary condition. All results previously published in [21] . ................................ ................................ ................................ ........ 15 Figure 8 . Left : Examples of coronal MRI images of the buttock from a healthy 29 - year - old female at a non - weight - bearing posture (top frame) and during weight - bearing sitting (bottom frame). Right : ( a ) Mesh of the subject - specific FE model (for subject #4). ( b ) An example of the fit between MRI measured (gray areas) and FE - calculated (meshed areas) gluteal muscle contours during weight - bearing sitting. The right and left frames show data for the gluteus muscles at the right and left body sides, respectively. Pictures are from [19] . ................................ .................... 16 Figure 9 . ( a ) 3D thigh - buttock FE model and the included muscle groups and bony structure. ( b - c ) Representation of the 30 regions of interest that were identified on the thigh - buttock surface and definition of the local reference system. ( d - e ) Sagittal view of the small area right beneath the ischial tuberosity that was selected to calculate the deformation shift of the gluteus muscle. ( d - e ) Representation of the selected area in the MRI and corresponding FE simulation, respectively. All figures are from [20] . ( f ) Results are from [63] . ................................ ................................ ........... 17 Figure 10 . ( a - e ) Three - dimensional CAD model of the left buttock and thigh. The model was created from MRI scan segmentation. This included: manual tracing over contours of muscles, fat and bones in transverse MRI scans using, and then creating a lofted geometry that connects the segmentation contours. ( a ) An isometric view of the 3D subject - specific CAD model. ( b ) A side view of the 3D subject - specific CAD model overlaid on to the sagittal non - weight - bearing MRI scan of the upper thigh. ( c ) The right view of the mesh ed muscles and bones, ( d ) the left view of the mesh for the fat layer and the sacrum, and ( e ) the left view of the meshed muscles and bones. ( f - h ) An illustration of the detailed FE model ( f ) and the simplified FE model ( h ). ( g ) The anatomical structures i ncluded in the segmentation of the simplified model. ( h ) The boundary and loading conditions applied to the simplified model. All figures are from [28] . ........................... 19 Figure 11 . ( a ) The buttock model geometry from [64] . ( b ) Cross - section of the FE model of seat - interface and buttock soft tissue from [27] . ( c ) The finite element model of buttock and seat surface from [60] . ( d - e ) Finite eleme nt model [61] during sitting position: ( d ) before loading. ( e ) after loading on stiff target. Plane stress was assumed on the boundary plane on this 2D model. 22 Figure 12. ( a ) Sample location and orientation on animals in Group I (left in the picture) and Group II (right in the picture). ( b ) Dimensions of dog bone shape skin sample. ( c ) A skin sample xi lution (HBSS) after it was cut out of the skin flap. ( d ) A mounted skin sample during the tensile test. The position of the four markers placed on the central part of the sample were tracked by a CCD camera to calculate its stretch ratio during the test. T he gauge length (i.e., the distance between the highest and lowest fiducial marker on the skin sample) was , as measured before the test. ( e ) Loading scheme of the uniaxial tensile test. RC: reference configuration. ................................ ................................ ................................ ......... 25 Figure 13. Mechanical parameters on a representative rupture stress - stretch curve of a rat skin sample. ................................ ................................ ................................ ................................ .......... 29 Figure 14. a ), were segmented into images with different components based on color s, i.e., ( b ) red (cytoplasm, muscle, erythrocyte), ( c ) blue (collagen), ( d ) black (nuclei). Collagen fiber bundles in a PSR stained rat back skin histology slide observed under polarized light, as shown in ( e ), were grouped based on colors, i.e., ( f ) red, ( g ) yellow, and ( h ) green, which indicated thick, medium, and thin fibers, respectively. Straightness of a single collagen fiber, as shown in ( i ), was measured in the following way as shown in ( j ). The green curve following the shape of the fiber bundle re presented the original length , while the white straight line connecting the two ends of the fiber bundle represented the compressed length . Straightness was calculated as . ..................... 31 Figure 15. Influence of storage condition on initial slope, maximum slope, UTS, rupture strain, and toughness. Values shown in the bar graph are means and stan dard deviations. Values shown below each bar graph are p - values of one - way ANOVA comparing the mechanical parameter difference between four different storage conditions. ................................ ................................ .. 36 Figure 16 ( a ) Averaged cyclic loading curves and ( b ) averaged rupture loading curves of each storage protocol. Horizontal bars on each curve represent standard errors. ................................ . 36 Figure 17. Influence of location, orientation and sex on ( a ) initial slope, ( b ) maximum slope, ( c ) rupture strain, ( d ) UTS, and ( e ) toughness. Means and standard deviations are shown in the bar graphs. P - values of the three - way ANOVA are shown below the bar graph for comparison within each factor, i.e., upper vs. lower, axial vs. transverse, male vs. female. Statistical significance ( and strong statistical significance ( are marked by single and double asterisks, respectively. ................................ ................................ ................................ ................................ .. 38 Figure 18. ( a ) left location on rat back. ( b ) Means and standard deviations of collagen fiber straightness at four back locations measured in PSR stained histology images. 50 fibers/5samples were measured for each locations. UL upper left, UR upper right, LL lower left, and LR lower right. ( c ) Means and standard deviations of relative collagen content at upper and lower back measured by collagen tained histology images. Statistical significance ( ) is marked using asterisk. ................................ ................................ ................................ ............... 39 Figure 19. Representative experimental data and constitutive model with corresponding best - fit constitutive parameters at upper and lower rat back locations. A xial experimental data are from xii animal 5, while transverse experimental data are from animal 10. Best - fit constitutive parameters of the simulated stress - stretch curves shown in the figure are as follows. Upper back sample in axial direction (animal 5) , , , ; and transverse direction (animal 10) , , , . Lower back sample in axial direction (animal 5) , , , ; and transverse direction (animal 10) , , , . ................................ ................................ ................................ ....................... 41 Figure 20. Influence of location, orientation, and sex on best - fit constitutive parameters, ( a ) , ( b ) , ( c ) , ( d ) . Means and standard deviations are shown in the bar graphs. P - values of a three - way ANOVA are shown below the bar graph for comparison within each factor, i.e., upper vs. lower, axial vs. transverse, male vs. female. Statistical significance ( and strong statistical significance ( are marked by single and double asterisks, respectively. ................................ ................................ ................................ ................................ ....................... 42 Figure 21 . ( a ) Schematic of the seating positioning for a participant during 3D scanning. ( b ) Thigh CAD model obtained from 3D scanning of the thigh of a representative participant. ( c ) Representative 3D FE model. ( d ) Dimensions and positioning of testing areas for the 3D FE model, on top, and the semi - 3D FE model, on the bottom. From left to right, proximal, middle, and distal thigh testing locations are s hown in red. ( e ) Flow chart of the constitutive parameter optimization protocol. ................................ ................................ ................................ ................................ ........ 51 Figure 22 . Compression force - deflection curves from experimental measures (symbols) and model predictions (lines) for each location and each constitutive model considered. ( a - d ) The same experimental data are shown in every section, representing the mechanical behavior for representative participant F6 at the proximal (triangles), middle (circles), and distal (asterisks) locations. Also shown are the model predictions after parameter optimization for each location, namely proximal (solid line), middle (dashed line), and distal (dotted line). The models considered are ( a ) neo - Hookean, ( b ) Mooney - Rivlin, ( c ) Ogden, and ( d ) Fung orthotropic. ........................ 56 Figure 23 . Mechanical behavior comparison between 3D and semi - 3D FE thigh models of participant F6 at ( a ) proximal thigh, ( b ) middle thigh, and ( c ) distal thigh. Horizontal axis represents forces applied to the soft tissue of thigh at each loading step, while the vertical axis represents deflecti on difference between semi - 3D FE model and 3D FE model at each loading step, i.e., . ................................ ................................ ................................ .................... 58 Figure 24 . Best - fit Ogden model and Fung orthotropic model parameters, i.e., ( a ) , ( b ) , ( c ) , ( d ) , ( e ) , and ( f ) for 20 tested subjects, at proximal, middle, and distal thigh locations, respectively. Parameter values are grouped based on location (marked at the top of each plot) and sex (marked at the bottom of each plot). The central mark within the box indicates the group median, and the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively. The whiskers extend to the most extreme data points not considered outliers, and the outliers are plotted individually using the '+' symbol. ................................ ................................ .. 60 Figure 25 . ( a ) A cylinder device that applies suction pressure on skin and measures the displacement of skin by an ultrasound transducer. . ( b ) The ascendant part of pressure - xiii displacement data of skin under suction pressure. The solid line represents experimental data, while the dotted line represents theoretical data. ( c ) FEM geometry to reproduce the suction test. The thickness of t he skin is 0.95 mm. The diameter of the outer circle is 30 mm, corresponding to the cylinder suction device shown in ( a) ( ). The diameter of the inner circle is 6 mm, corresponding to the aperture size shown in ( a) ( ) . ( d ) Mesh of the F EM of the skin. ( a ) and ( b ) are from [121] . ................................ ................................ ................................ ............ 70 Figure 26 . ( a ) A cross - section MRI image of the right thigh of the subject. ( b ) Isolated pixels of the muscle tissues from the cross - section MRI image. Femur geometry will be subtracted from the muscle geometry when the 3D models of femur and muscle are formed. ( c ) Boundary of the muscle tissues recognized from the cross - section muscl e tissues. ( d ) Boundary lines of the muscle tissues stacked in the 3D space. ( e ) 3D geometry of the muscle tissues generated from the cross - section boundary lines. ( f ) 3D geometry of the thigh with fat, muscle, and femur. ( g ) Mesh of the thigh 3D geometry. ................................ ................................ ................................ ....................... 72 Figure 27 . indentation test areas at ( a ) proximal thigh, ( b ) middle thigh, and ( c ) distal thigh. ( d ) The centr al part of skin that is under suction pressure. ........................... 73 Figure 28 . ( a ) Structure of two crossing collagen fiber families, where and represent the preferred collagen fiber directions of collagen family 1 and 2, respectively. The bisector of the direction and direction is defined as the preferred fiber orientation of the skin, and is the angle between the direction an d the preferred fiber orientation of skin. ( b ) The wrinkle line of the human leg that indicates the mean fiber direction. At the anterior side of the thigh, wrinkle lines run ~37° with respect to the circumferential direction. At the posterior side of the thigh, wrinkle lines run in transverse direction. Around the knee area, wrinkle lines run in transverse direction. ( b ) is from [57] . ................................ ................................ ................................ ............ 75 Figure 29 . Preferred fiber orientation vector of each skin shell element at ( a ) anterior side of the thigh, ( b ) posterior side of the thigh, ( c ) knee area. ................................ ................................ 76 Figure 30 . ( a ) Force - deflection behavior of thigh soft tissue from in vivo indentation data in quadruped posture and the FE simulation. ( b ) Force - deflection be havior of thigh soft tissue from in vivo indentation data in prone posture and the FE simulation. ................................ ................. 79 Figure 31 . ( a ) Displacement con tour map of the skin suction model with best - fit constitutive parameters with maximum pressure load. The displacement unit in the figure is millimeter. ( b ) Pressure - displacement curves of experimental data collected in vivo and FE simulation data with best - fit constitutive parameters. ................................ ................................ ................................ .... 80 Figure 32 . Force - deflection behavior of thigh soft tissue collected from in vivo indentation test and FE simulation. The in vivo force - deflection curve in ( a ) and ( b ) are the same data from the quadruped posture. The FE simulated force - deflection curves in ( a ) are from thigh FE model with in vitro skin constitutive parameters. The FE simulated force - deflection curves in ( b ) are from thigh FE model with in vivo skin constitutive parameters. ................................ ........................... 81 1 INTRODUCTION AND LITERATURE REVIEW 1.1. Overview Pressure ulcers (PUs), also known as pressure sores, are localized damage to the skin and underlying tissues, usually occurring over a bony prominence and caused by sitting or lying in one position for long time. People with reduced mobility are critically vulnerable to PUs, including elders, people with spinal cord injuries, and people temporarily disabled due to diseases or medical treatments. PUs are a detriment to the well - being of people who lose their mobility either permanently or temporari ly, and high morbidity and mortality are associated with P U s [1] . Each year in the US an estimated 2.5 million people develop PUs [2] . In the year 2013, almost 30,000 dea ths were caused by complications associated with this condition globally [3] . Further, the healing time of PUs can be long, ranging from several weeks to several months, depending on the stage of the ulcer and factors like age and other medical complications. Even when healing occurs, it is often not complete, leading to a high reoccurrence rate [4] . Several hypotheses have been made on the biomechanical mechanism involved with the formation of PUs. These include cellular death due to mechanical distortion [5] , tissue decay due to reduced i nterstitial flow and lymphatic drainage [6] or reduced blood perfusion [7] , and localized is chemia [8], [9] . Although the initiating mechanism of PUs is still unclear, it is commonly accepted that in ternal normal and shear stresses, due to the presence of unrelieved external loads, play a central role in the formation and development of these wounds [9] [12] . Therefore, it is crucial to evaluate the internal stresses caused by the physiological loading associated with sitting or lying down to assess the risk of tissue injury. Despite the significance of i nternal stresses in PUs formation, interfacial pressures, which are a surface measure of stress, are the indicators commonly used to develop practices and protocols to minimize loading on the soft 2 tissue, for example for timing between posture changing o r for the selection of materials for hospital mattresses [13], [14] . However, previous published studies, that focused both on animal [15] and numerical models [16], [17] , have pointed out that no direct correlati on exists between interfacial pressure and internal stresses of soft tissue. Therefore, tools and methods that can show needed. The ability of finite element (FE) models to accurately represent the anatomical structure of the leg and buttocks area and to estimate the localized stress/strain fields within highly deformable media, makes them powerful tools to investigate soft tissue response to external loadings [1 8] . Indeed, FE models were employed in previous studies to identify and investigate external and internal factors that may affect stress/strain distributions within the soft tissue during sitting [19] [21] , a position in which both disabled and normal people sp end majority of time during the day [22] [25] . Despite the significant advancement previous studies have achieved, there are still important aspects in human thigh - buttock soft tissue modeling area that need to be improved. Two challenges are identified in this dissertation: 1) Microstructurally moti vated skin modeling for an individual skin layer in finite element model. 2) Parameters estimation associated with large deformations. 1.2. Skin mechanical model In the past 20 years, there has been a growing effort on developing more anatomically accurate FE m odels, however, one component that has still been significantly overlooked in this effort is the skin. This might be due to the fact that FE models have been mostly used as a tool to assess the stresses and strains in deep tissue, because deep tissue injur ies are harder to detect compared to superficially formed PUs. However, superficially formed PUs, which are more involved with 3 friction and shear forces are also a major form of PUs [26] . Previous FE studies either model the bulk properties of all soft tissues ( i.e., skin, fat, muscle) together [21], [27] , or model skin and fat as one component [19], [28] . The only few studies that modeled the skin separately used simple phenomenological constitutive models [16], [20] . A separate, accurate representation of the skin in a FE model is important to investigate PUs formation due to the following reaso ns: (1) Different mechanical properties of tissues . Skin is a very thin layer compared to muscle and fat. Collagen fibers in skin, which are the main weight - bearing structure, have preferred directions parallel to the surface [29] . Hence, mechanical response from the skin layer may come from extension rather than compression during sitting contact. (2) Differ ent vulnerability to pressure . Due to the differences in composition, structure, and metabolic rate, different tissues ( i.e., muscle, fat and skin) may have different pressure tolerance [30] . Thus, str ess or strain threshold criteria for PUs formation should be considered for different types of tissues separately, both in experimental testing and in numerical simulation. Skin microstructure. major barrier, both biologically and mechanically, against the external environment by fulfilling multiple biological purposes ( i.e., regulation of heat and water exchange with the surroundings, protection from mechanical, bacterial or viral insults). Its structure can be approximately subdivided in three principal layers: the epidermis, the dermis and the hypodermis. The outermost layer (epidermis) is a thin cellular membrane that serves as a barrier but does not contribute to the mechanical behaviors. Th e middle layer (dermis) is mostly composed by a matrix of collagen and elastin fibers , 4 Figure 1 . ( a ) Schematic of skin structure. ( b ) A typical load - strain graph obtained in an in vitro uniaxial test of human abdominal skin, three phases of extension are indicated in the graph. ( c - f ) Collagen network, parallel section through mid - dermis, under the scanning electron microscope at different strain levels. ( c ) Unstressed skin, collagen fibe rs appear wavy and multidirectional. ( d ) Strain level corresponding to the end of phase 1, fibers have reoriented towards load axis. ( e ) Strain level corresponding to the end of phase 2, fibers appear aligned in the direction of the load and almost straigh t. ( f ) Strain level corresponding to mid - phase 3, fibers appear aligned, straight, and compacted. Pictures ( b - f ) from [31] , numbering changed from original. and endows strength and compl iance to the skin. The innermost layer (hypodermis) serves as a supportive substratum and is primarily made of adipose tissue ( Figure 1 a ). Collagen is the major load - bearing constituent of the dermis and accounts for approximately 60%~80% of skin dry weight [32], [33] . Brown [31] examined the structure of human skin under scanning electron microscopy. His study showed that the nonlinear stress - strain response of skin ( Figure 1 b ) is due to the re - orientation, alignment, and straightening of undulated collagen fibers during extension ( Figure 1 c f ). A study by Finlay [34] showed that anisotropy of skin is attributed to the non - uniform dis tribution of collagen fiber orientation. These observations demonstrate the important role collagen fibers play in skin mechanical behavior. As a result, the quantification of the collagen network structure is considered a key step in the development of a microstructural model of skin 5 Figure 2 . Computer - generated images ( a, c ) and the power plots of Fourier analysis ( b, d ) previously published in [35] . Specifically, ( a ) and ( c ) show two examples of the computer - generated images , which were used to validate the accuracy of Fourier analysis generated power plots in representing fiber orientation p reference. The images consist of ellipses with a width/leng th ratio of 1:10. ( b ) and ( d ) are the power plot generated by Fourier analysis, of ( a ) and ( c ), respectively. The network represented in ( b ) corresponds to an orientation index of 1, and the network represented in ( d ) corresponds to an orientation index of 0. ( e - h ) Laser scatter images and first - order power plots previously published in [36] . ( e ) The laser scatter image of lesional sclerodermic skin (LS) produces an elongated scatter plot, indicating a higher orientation of collagen bundles when compared to control skin (CS) , shown in ( f ), which produces a circular scatter plot. ( g - h ) First - order power plots from the fast Fourier transform of the images shown in ( e - f ). ( g ) LS creates a power plot with a preferential distribution with two maxima, and CS, sh own in ( h ) is characterized by a power plot with a circular pattern. Numbering changed from original studies. mechanical behavior. For example, van Zuijlen et al. [35] developed an image analysis technique that used Fourier analysis to assess collagen orientation on image of human skin acquired by confocal laser scanning microscopy. The fiber orientation index was evaluated by calculating the width to length ratio of the power plot of histology images generated by Fourier analysis ( Figure 2 a - d ). The orientation index ranges between 0 an d 1, which represents the cases of perfectly aligned fibers and randomly distributed fibers, respectively. In another study, Vries et al. [36] used two methods, specifically the laser scatter method and the fast Fourier transform (FFT), to describe collagen orientation, then collagen structure indicators, such as orientation ratio, orientation variation, and bundle spacing, were evaluated based on the laser scatter images and FFT power 6 plots ( Figure 2 e - h ). These studies, however, provided mostly qualitative measurement that indicated the preferential directi on or level of anisotropy of the collagen network, rather than quantitative measurement. Studies that quantified collagen orientation distribution which could be implemented into structural constitutive model are still rare [37], [38] . Skin constitutive modeling. Skin is a nonlinear, anisotropic, and visc oelastic material and different models have been developed to describe several aspects of skin mechanical behavior [38] [45] . The literature review that follows discusses papers that mainly focus on two aspects: (1) the development of nonlinear and anisotropic constitutive model of skin, and (2) experimental measures of collagen microstructural organization. Skin constitutive models can be categorized into two types: phenomenological models and microstructural model. Phenomenological models focus on finding a mathematical form that describes the stress - strain experimental data, without making hypothesis on th e physiological meaning of the each parameter. For example, Alexander and Cook [39] , first, mechanically tested in vivo the skin from the upper back of healthy adult males using a suction cup device and obtained biaxial mechanical data. Second, they described these data by employing a relationship of the form , ( 1 ) where is the skin tension (stress) in testing direction , and are extension ratios in the test direction and in the perpendicular direction, respectively, and and are material constants with the dimension of a tension and dimensionless, respectively. In another study, Tong and Fung [44] proposed a strain energy function of the form , ( 2 ) 7 where , , and are constitutive parameter s , , with the dimension of a stress, while , dimensionless, and are the components of the Green strain tensor. In this study, two stress - strain experimental datasets were collected on rabbit skin in two perpendicular directions, each set was collected with one principal direction of stretch and the perpendicular direction fixed. The a uthors determined the constitutive parameter s by fitting the experimental data at two selected points, one at a low stress level and one at a high stress level. Following the idea proposed by Tong and Fung to employ an exponential functional to describe the mechanics of skin, Veronda and Westmann [45] developed a strain energy function of the form , ( 3 ) where , , are constitutive parameter s, the first two with the dimension of stresses and the last dimensionless, , and are, respectively, the first, second, and third invariant of the right Cauchy - Green strain tensor, , where is the deformation gradient. Then, the Authors determined the best - fit constitutive parameter values by optimizing the form shown in Equ ation ( 3 ) to describe uniaxial stress - stretch data collected on cat skin. Phenomenological models provide accurate description of the experimental data , however, they offer little insight on skin microstructure. In contrast, microstructural models connect macroscopic level mechanical response of skin to its microstructure by associating constitutive parameters to physiologically meaningful measurements. 8 Figure 3 . ( a ) Geometry of a model collagen fiber as described in [46] . ( b - c ) Schematics of the collagen and elas tin network structure in flat tissue from [42], [43] . ( b ) Tissues with high density of crosslinks and elastin induced collagen undulation. ( c ) Tissues with low density of crosslinks and inherent collagen undulation. Thick lines represent collagen; thin lines represent elas tin, and dotted lines represent overall collagen fiber direction. As one of the early attempts to develop microstructural model, Comninou and Yannas [46] modeled fiber undulation as long sinusoidal beams ( Figure 3 a ). The model however did not include c ollagen structure experimental data to determine the geometrical parameters of the sinusoidal beams, further the effect of fiber orientation was only briefly discussed in a qualitatively manner. Then, Lanir [42], [43] proposed a modeling framework that described skin structure as a ground matrix in which collage n and elastin fiber embedded ( Figure 3 b - c ). The assumptions made in this framework have crucial influence on later skin microstructural models [47] . The total strain energy function was assumed be to the sum of strain energy functions of each component. The fibers were m odeled to have no compressive strength and to be buckled under zero load. Further, each fiber was considered to have the capability of supporting only a uniaxial strain, on. The undulation and orientation of the fibers were described using probability density functions. 9 In another study, Groves et al. [40] developed a microstructu ral model to compare human and murine skin mechanical characteristics. The model described the skin as a composite structure made of ground matrix and three families of parallel fibers. The ground matrix was described by the Veronda - Westmann material shown in Eq uation ( 3 ) . Three fiber families were included in the model, the mechanical behavior of each fiber family was described by employing the strain energy function developed in [48] . In this study by Weiss et al. , however, there is a lack of microstructural experimental data, which resulted in the necessity of increasing the number of constitutive parameter s, 4 for each collagen family and 14 for the model overall, to ensure a good mechanical description. The redundancy in constitutive para meter offered good fitting results to experimental data, but somewhat dimi nished the physical meaning of each parameter. A later study from Karimi et al. [41] had similar issues, the model employed 4 fiber families, and although the mean orientation angle of the collagen fibers was measured experimentally, it was still not enough to characterize the collagen structure in details in the constitutive model. Fig ure 4 . Left : Maximum intensity projection for a 60 transverse section showing collagen fibers in bright red. Right : A mixture of two von Mises distribution fitted to fiber orientations presented in the form of histogram, under various symmetry assumptions. Pictures are from [37] . As discussed in the section Skin microstructure , quantifying collagen structure is a key step to develop accurate skin microstructural models that have physiologically meaningful parameters. Jor et al. [37] used a structural - tensor approach to determine the collagen spatial distri bution of 10 porcine skin. Orientation data were fitted to a mixture of two von Mises distributions ( Fig ure 4 ). In a later study, Jor et al. [49] tested porcine skin in vitro using a multi - axial tensile test. With the collagen structural data parameters, a modeling framework was proposed to simulate the experiments. The fibers were assumed to be linearly elastic and the ground matrix was assumed to behave like an isotropic neo - Hookean material. One drawback of these studies is that the collagen fiber microstructure was examined in a plane perpendicular to the skin surface, while the preferr ed orientation of collagen fibers in the dermis is parallel to skin surface [29] . More recently, Ní Annaidh et al. [38] developed an automated image analysis routine ( Figure 5 a - f ) that can calcula te the orientation distribution of collagen network using Van Gieson stained histology slides of human back skin. The image analysis routine can distinguish each single collagen fiber in the field and fit an ellipse along the fiber direction. The direction of the major axis of the fitted ellipse represented the orientation of the collagen fiber. The standardized - periodic Von Mises distribution was used to fit the probability density of collagen fiber orientation distribution ( Figure 5 g ). Finally, the resulting structural parameters were then employed in the Gasser - Ogden - Holzapfel model [50] to fit the uniaxial tensile data. The model was not only cap able of fitting the known experimental data, but also able to provide good prediction of stress - stretch response of skin sample along different direction as long as the unique structural data of that sample were provided. 11 Figure 5 . Images output from automated algorithm. ( a ) Original histology slide. Scale bar is 1 mm. ( b ) Binarized image after automated thresholding. ( c ) Binarized image after erosion step. ( d ) All identified collagen bundles outlined in green. ( e ) Remaining bund les which meet area and eccentricity criteria. ( f ) Best fit ellipse about each fiber that meets the specified criteria. ( g ) Histogram of collagen orientations. The two distinct peaks correspond to the preferred orientation of the two fiber families. All im ages are from [38] . Based on image analysis routine shown in [38] , this dissertation will further improve the automated ro ( i.e., binarization step). Collagen undulation level will also be experimentally measured and implemented into the microstructural model. Preferred collagen fiber orienta tion . The studies and Ní Annaidh group showed, through histolog ical data, the existence of a two peak orientation distribution of collagen fibers within the skin dermis. This suggests that, w hen implementing a microstructurally accurate skin layer into a FE model, we need to determine the preferred collagen fiber orientation across the human body , and specifically with in the region s of interest. In this context, we consider the angle associated with the bisector of the two peak distributi on as the preferred collagen fiber orientation angle . It has been shown before that t h e preferred collagen orientation angle is closely related to the skin tension lines [51] . - known skin tension lines mapping on human body. However, Lang n misquoted and his 12 illustrations incorrectly reproduced [52] . This fact can be shown by simply searching line on a search engine lik e Google when a variety of line mappings on human body return as results , some of them showing drastic differences between one another . Indeed, i t is of great importance to determine the proper mapping of tension lines for the implementation of collagen fibers in the skin layer FE model. In 1834, G uillaume Dupuytren [53], [54] for th e first time made the observation that elliptical wounds were produced by a round instrument on human body . Inspired by Dupuytren observation, in 1861, Langer [55] used rounded shaped instruments to stab through skin of cadavers , providi ng a thorough mapping of the long axis direction of the elliptical shaped wounds ( Figure 6 a c ) . This mapping became of great importance in the context of surgica l procedure, because surgeons were seeking to identify techniques that resulted in optimal healing and least conspicuous appearance of the surgical incision . It became apparent that the proper direction of the incision is the one that allow s the collagen f ormed under the scar to be aligned along the same direction with collagen of neighboring tissues , and that the mapping provided by Langer seemed to be related to this information . In 1935, however, Webster [56] first suggested surgeons should . that lines used by surgeons for e In 1951, Kraissl [57] further provided a complete mapping of wrinkle lines of the whole body ( Figure 6 d - f ) . Wrinkle lines can be demonstrated by existent wrinkles, gentle compression of the relaxed skin, muscles contraction or a combination of these methods. Later on, Borges and Alexander [58], [59] described the relaxed skin tension lines (RSTL s ) , which can be obtained through pinching test. The difference between RSTLs and wrinkle lines are subtle and reconcilable, 13 while [56] . Nonetheless, RSTLs a nd wrinkle lines are [53] . Figure 6 . ( a ( b c ) d ) Wrinkle lines on the posterior trunk and scalp . ( e ) Wrinkle lines on the anterior leg. ( f ) Wrinkle lines on the posterior leg. When it comes to the question that which tension line mapping should be used for collage n fiber implementation in skin for a FE model, we argue that wrinkle lines are the most appropriate for the model that is proposed in this dissertation. W rinkle lines and RSTLs have been estimated on the skin of a living person, while were o btained through study on skin cadavers . For this reason, wrinkle lines and RSTLs can take into account t he action of muscle and movement of joints , which play an important role in caus ing skin tension lines. The static state of the cadaver and the relaxed state of underlying tissue s af ter death m ight significantly change the tension line distribution compared to a living person. Furthermore , b y using x - ray diffraction and micoradiographic technique, Homstrand [51] demonstrated that the major part of the collagen is parallel with the wrinkle lines . ing from surgical incisions the preferred collagen fiber directions. Therefore, the preferred fiber orientation within the FE model skin layer will be determi ned by wrinkle lines in this dissertation. 14 1.3. FE models of leg and buttock area Previous FE studies have investigated the influence of both internal and external factors that affects stresses/strains distribution in thigh - buttock area during sitting. Dabnichk i et al. [60] showed that, when sitting on a rigid surface, if friction occurs at the body/seat interface the stresses in the tissue adjacent to the seat increased, but at the same time the maximum stress, which occurred some distance away from the contact area, decreased. As regard to cushion, it h as been concluded by previous studies that a softer cushion not only reduces the interface stress substantially [16] , but also contributes to decreasing internal stresses [61] . Furthermore, it has been shown that a thicker cushions can effectively reduce seat - interface and subcutaneous compressive stress [27] . Regan et al indicated, however, the existence of an optimal cushion thickness, the showed that increasing cushion thic kness over a certain level could compromise postural stability, leading to asymmetrical loading, elevated seat - interface pressures, and loss of function [27] . Not only the mechanical properties of the seating surface, but also the mechanical characteristic of soft tissues play an important role in the magnitude and subdermal stresses/strains distribution. By employing uniaxial tension testing and FE model on rat gracilis muscles, Linder - Ganz and Gefen concluded that injured muscles were significantly stiffer than uninjured ones as a result of exposure to prolonged and intensive loads [15] . This increased muscle stiffness in widening region could lead to elevated tissue stresses that exacerbate the potential for tissue necrosis. Nonlinear mechanical behavior data in vivo and parameters estimation . Despite all the findings about contributing factors to PUs formation from previous FE studies, one critical aspect that is still missing is the address of nonlinear mechanical behavior of FE model s. There are two aspects involved to validate the nonlinear mechanical behavior of the FE model: (1) in vivo experimental data exhibiting the full range (from being relaxed to being fully compressed) nonlinear mechanical beh avior need to be collected. (2) A parameter estimation process that 15 actively reduces the difference between the experimental data and corresponding FE output is needed to achieve the accurate constitutive parameter s for the FE model. However, both of the two aspects were largely overlook ed in the literature. To further illustrate this issue, the model development and validation process of some state - of - the - art FE studies focusing on leg and buttock area are discussed as below. Figure 7 . FE model of the pelvis an d thighs published in [21] . ( a ) Bony structures, ( b ) bony structures and human soft tissues, and ( c ) complete model ( i.e., bony structures, human tissues, and skin). ( d - e ) P ressure distributions (Pa) predicted by the FE pelvis model (shown in d ) and measure experimentally for a human male (shown in e ) for a rigid surface boundary condition. ( f - g ) Pressure distributions (Pa) predicted by the FE pelvis model (shown in f ) and me asure experimentally for a human male (shown in g ) for a soft cushion boundary condition. All results previously published in [21] . Verver et al [21] developed a 3D buttock - thigh sitting model ( Figure 7 a - c ) to investigate the interface pressure distribution during sitting and how variations in human tissue and cushion properties would affect these distributions. The buttock - thigh outer geom etry and bone geometry were collected from two different subjects. The skin was modeled using triangular shells elements, while fat and muscle were modeled together as a whole using 4 node tetrahedron elements. The skin was described employing a linear ela stic isotropic constitutive model , and bulk properties of 16 fat and muscle were described using the Mooney - Rivlin constitutive model . Finally, all the constitutive parameter s were chosen from previously published work. This model was validated by comparing i nterfacial pressure fields predicted by the FE model and measured experimentally, associated with siting on a rigid surface ( Figure 7 d - e ) and on a soft cushion ( Figure 7 f - g ). The Authors concluded that FE model showed reasonable correlation with the human data, pressure distribution differences between FE model and in vivo data we re not quantified, and the Authors optimization process implemented in the study to attempt to minimize the difference between FE model and in vivo data. Figure 8 . Left : Examples of coronal MRI images of the buttock fr om a healthy 29 - year - old female at a non - weight - bearing posture (top frame) and during weight - bearing sitting (bottom frame). Right : ( a ) Mesh of the subject - specific F E model (for subject #4). ( b ) An example of the fit between MRI measured (gray areas) and FE - calculated (meshed areas) gluteal muscle contours during weight - bearing sitting. The right and left frames show data for the gluteus muscles at the right and left body sides, respectively. Pictures ar e from [19] . In another study, Linder - Ganz et al. [19] developed subject - specific semi - 3D model, defined has having one dimension significantly smaller when compared to the others ( i.e., the mod el has a thickness of 4 mm) for 6 subjects. The model geometry was constructed based on MRI images obtained from a single plane. Muscle and fat were modeled separately and skin was modeled as 17 part of the fat layer ( Figure 8 ). The FE models were then validated using two indicators: (1) the contour of the gluteal muscle from the FE solution was compared to that measured from MRI images, and (2) the peak interface pressure under the ischial tuberosity predicted by the FE analysis was compared to that measured experimentally, for each subject. This study used quantifiable indicators ( i.e., correlation coefficient and - values of two - tailed probability) to describe the o utput difference between FE model and in vivo data. However, there was still no optimization process to minimize the error between FE model and in vivo data. This model was then employed in a later study to investigate the strains and stresses differences in buttock subdermal tissue between paraplegic participants and healthy participants [62] . Figure 9 . ( a ) 3D thigh - buttock FE model and the included muscle groups and bony structure. ( b - c ) Representation of the 30 regions of interest that were identified on the thigh - buttock surface and definition of the local reference system. ( d - e ) Sagittal view of the small area right beneath the ischial tuberosity that was selected to calculate the deformation shift of the gluteus muscle. ( d - e ) Representation of the selected area in the MRI and corresponding FE simulation, r espectively. All figures are from [20] . ( f ) Results are f rom [63] . 18 A work published by Makhsous et al. [20] focused on the development of a 3D thigh - buttock FE model with detailed muscle geometry. 5 muscle groups were identified in the model, which al so contained separate fat and skin layer. In this study, fat and muscle layers were meshed using 4 nodes tetrahedral elements, while skin was meshed using 3 nodes triangular elements. All the soft tissues mechanical behaviors were described employing Moone y - Rivlin constitutive model s, and the constitutive parameter s were obtained for each component from literature. The FE model was validated comparing the values of two quantities, specifically the gross displacement of 30 regions of interest ( Figure 9 b - c ) and shift of gluteus muscle ( Figure 9 d - e ), between FE and in vivo data during the compression introduced by an inflated bladder. This is the first study that developed a FE model that could accurately represent the soft tissues anatomical structure. The model developed here was then employed in a later study to evaluate the e ffectiveness of a new - partly transfer the loa ds from the ischial tuberosities to thigh soft tissues [63] ( Figure 9 f ) . 19 Figure 10 . ( a - e ) Three - dimensional CAD model of the left bu ttock and thigh. The model was created from MRI scan segmentation. This included: manual tracing over contours of muscles, fat and bones in transverse MRI scans using, and then creating a lofted geometry that connects the segmentation contours. ( a ) An isom etric view of the 3D subject - specific CAD model. ( b ) A side view of the 3D subject - specific CAD model overlaid on to the sagittal non - weight - bearing MRI scan of the upper thigh. ( c ) The right view of the meshed muscles and bones, ( d ) the left view of the mesh for the fat layer and the sacrum, and ( e ) the left view of the meshed muscles and bones. ( f - h ) An illustration of the detailed FE model ( f ) and the simplified FE model ( h ). ( g ) The anatomical structures included in the segmentation of the simplified model. ( h ) The boundary and loading conditions applied to the simplified model. All figures are from [28] . In a recent study, Al - Dirini et al. [28] developed a subject - specific model through MRI images ( Figure 10 a - e ). The model geometry composed of very detailed muscle structure, which contained 28 muscles, and inter - muscular fat, fascia, and subcutaneous fat were modeled together. Tetrahedron elements were used to mesh the geometry, and first order Ogden model was employed to describe the material behavior. In this study, for the first time, a pa rameter optimization process was employed. Briefly, a simplified model containing only the tissues surrounding the ischial 20 tuberosity and the lesser trochanter was isolated from the full 3D model ( Figure 10 f - h ). Initial Ogden parameters were obtained through literature. Then Genetic Algorithm was used to search for proper constitutive parameter s by minimizing the difference between FE estimates and MRI measurements on 1 - D thickness change of the Gluteus Maximus, intermuscular fat, subcutaneous fat below ischial tuberosity and the lesser trochanter. The FE model with optimized constitutive parameter s was further validated by comparing the FE calculation and MRI measurem ents for the displacement of gluteus maximus muscle. This study for the first time utilized an optimization process to estimate constitutive parameter s, however the in vivo data used for the comparison only included two loading conditions, weight - bearing and non - weight bearing conditions. Due to the nonlinear nature of soft tissue, data collected in only two loading conditions might not be enough to characterize the nonlinea r mechanical behavior of soft tissue. As discussed above, there is still a need for experimental data collected in vivo that can represent mechanical behavior of soft tissues in large deformation, and for an optimization procedure that can estimate best - fi t constitutive parameter s in a 3 - dimensional FE model. 2D, semi - 3D, and 3D FE model . Geometry of the early FE models were largely simplified: the buttocks were modeled as a hemisphere [64] or a slab [27] , and the ischial tuberosity were modeled as a cylinder with hemispherical end - point [60] , as shown in Figure 11 a - c . The FE models in these studies were usually further simplified to 2D models by assuming axisymmetry, while other studies developed 2D models from the start by assuming plane strain or plane stress condition s [60], [61], [65] , as shown in Figure 11 d - e . More recently, researchers were able to develop anatomically accurate 3D models for the thigh and buttock areas by employing medical images [20], [28], [66] . Despite increasing computational capabilities, 2D and semi - 3D FE models still remain popular among researchers [19], [67], [68] . The model from Linder - Ganz et al. [19] 21 ( Figure 8 ) is a typical semi - 3D model: it has three dimensions, but the measurement in one dimension (4 mm in thickness) is significantly smaller than the other two. All elements had direct contact to the plane strain boundary condition. Therefore, semi - 3D FE model would behave more like 2D model than 3D model. Considering its nearly incompressible nature, soft biological tissue will significantly expand/shrink in the transverse directions when subjected to the large deformations associated with seating. It has been shown by rece nt studies that soft tissues in the leg - buttock area are subject to three - dimensional deformations. Al - Dirini et al. [69] measured deform ations below the ischium and the proximal femur across 6 seated subjects through MRI scans on both transverse and sagittal planes. Their results showed that deformation of gluteus maximus muscle in the direction along the thigh was almost equivalent to the deformation below the ischium perpendicular to the sitting surface. In a typical semi - 3D model ( Figure 8 ), deformation of gluteus muscle in the direction along the th igh is lost. Sonenblum et al. [70] conducted mul ti - planar scans on the buttocks of 7 subjects. Their results also show that a single slice of MRI hinders the ability to identify tissue accurately. Therefore, 2D and semi - 3D FE models are limited in their ability to accurately represent the in vivo mechan ical behavior of the soft, due to the influence of boundary conditions. 22 Figure 11 . ( a ) The buttock model geometry from [64] . ( b ) Cross - section of the FE model of seat - interface and buttock soft tissue from [27] . ( c ) The finite element model of buttock and seat surface from [60] . ( d - e ) Finite element model [61] during sitting position: ( d ) before loading. ( e ) after loading on stiff target. Plane stress was assumed on the boundary plane on this 2D model. Constitutive model choice. Early FE studies commonly emp loyed linear elastic material description for soft tissue. For example, Chow & Odell [64] measured from synthetic gel from literature for their buttock model. Todd & Thacker [71] obtained force - deflection curve. Ragan et al [27] in their study. With the growing computation power and FE platform advancement, the constitutive model s applied in the FE models have shifted to nonlinear hyperelastic material. Fo ur hyperelastic models can commonly be found in literature, namely neo - Hookean model [19], [61], [66], [68] , Mooney - Rivlin model [20], [21], [60] , first order Ogden model [16], [28], [65] , and Fung - type model [72], [73] . However, a constitutive model being nonlinear does no t guarantee that it is able to match the nonlinear behavior of certain soft tissues. Different nonlinear models have been developed to describe specific soft tissues, but might not be suited to describe others. Due to the lack of in vivo data in large defo rmations and protocols to estimate constitutive parameter s, 23 currently no guidelines exist for the choice of constitutive model when it comes to soft tissue. Researcher have acknowledged the need for proper nonlinear constitutive model s, for example Moes & Horváth [66] pointed out that a constitutive model with increasing stiffness was required to describe the behavior of soft tissues. Also, Makhsous et al. [20] suggested that the primary factor causing the difference between their model prediction and the measurement on MRI images is th e lack of an optimal constitutive model and precise material properties in the FE analysis. The apparent nonlinearity of human soft tissue cannot be precisely described by the Mooney - Rivlin model they used in that study. Taken together, there is a pressing need for proper constitutive model choice in human soft tissue modeling. 1.4. Objectives of this dissertation Taken together all the remaining issues and challenges in the literature, we determine the objectives of this dissertation to be as follows: (1) Propose a micro structurally based skin constitutive model. Quantify collagen network structure and incorporate it into the structural model. (2) Develop an optimization procedure that is able to match nonlinear behaviors between FE simulation and in vivo experimental d ata. Quantify the different between 3D and semi - 3D FE model in an individual - specific case. Compare the accuracy of four commonly used constitutive model representing soft tissue nonlinear mechanical behavior. (3) Construct a thigh FE model that has detailed anatomical representation of different soft tissue types, i.e., skin, fat, and muscle. Identify constitutive parameters for the different soft tissue types that can describe the in vivo experimental data accur ately. 24 CHARACTERIZATION OF THE MECHANICAL PROPERTIE S OF SKIN USING RAT AND A MICRO STRUCTRUALLY BASED CONSTITUVIE MODEL 2.1. Overview The goal of this chapter is to develop a microstructural constitutive model that links collagen microstructure and skin macroscopic level mechanical behavior. Some import ant aspects of mechanical properties of skin were investigated, and a microstructural model of skin was developed in this c hapter. First, storage condition effects on skin mechanical properties were studies. We developed three different storage protocols. The mechanical parameters of skin samples stored in those conditions were evaluated and compared to those of the fresh samples. The best storage protocol among the three was determined. Effects of orientation, location, and sex on skin mechanical propertie s were also investigated. Automated image analysis routine were developed to quantify collagen structure, i.e., fibers orientation distribution and collagen relative content, using histology images . Those quantified collagen structure data were then implem ented into the microstructural constitutive model. 2.2. Methods 2.2.1. Mechanical testing Sample collection . Skin samples from sixteen Sprague Dawley rats (eight males and eight females, 10 to 12 weeks old) were collected for mechanical tensile tests. The dorsal skin of each animal was shaved after sacrifice and the entire back skin flap was excised . A custom - made dog bone shaped die was used to punch samples out of the flap (see Figure 12 b - c for sample dimensions). Four samples were collected from each animal, two from the upper back, on each side of the spine (i.e., upper left and upper right locations), and two from the lower back, on each side of the spine (i.e., lower left a nd lower right locations). The sixteen animals used for 25 mechanical tensile tests were divided into two groups of eight animals (four males and four females) each: in Group I, samples were cut in the axial direction (i.e., direction parallel to the spine) ; in Group II, samples were cut in the transverse direction (i.e., direction perpendicular to the spine, see Figure 12 a ). Finally, the thickness of each sample was measure d employing a Fowler digital Vernier caliper. Samples from male rats have thickness of and samples from female rats have thickness of . Figure 12 . ( a ) Sample location and orientation on animals in Group I (left in the picture) and Group II (right in the picture). ( b ) Dimension s of dog bone shape skin sample. ( c ) A skin sample placed in ( HBSS ) after it was cut out of the skin flap. ( d ) A mounted skin sample d uring the tensile test. The position of the four markers placed on the central part of the sample were tracked by a CCD camera to calculate its stretch ratio during the test. The gauge length (i.e., the distance between the highest and lowest fidu cial mark er on the skin sample) was , as measured before the test. ( e ) Loading scheme of the uniaxial tensile test . RC: reference configuration. Storage condition . To investigate the influence of storage condition on skin mechanical properties, th ree different storage protocols were developed as follows, Protocol I : samples were 26 Protocol II : samples were flash frozen in isopentane that was cooled down by dry ice ( - 7 8.5°C). The flash frozen samples were then stored in a - 80°C chamber for up to two weeks , and thawed at 4°C for 6 h ours before testing. Protocol III : flash - freezing and storage temperature were the same as Protocol II , but samples were thawed at 4°C for 24 h ours before testing. Control: samples were tested fresh, within 2 hours of excision . Skin samples obtained from each animal in Group I were randomized into the four groups (three storage protocols plus control), while skin samples obtained from animals i n Group II were stored using Protocol II. After the mechanical testing, an a nalysis of variance (ANOVA) was carried out to examine the influence of storage condition on skin mechanical parameters. While t he averaged stress - stretch behavior difference betw een the control group and each storage protocol was evaluated using normalized root - mean - square deviation ( N RMSD ) , which was calculated as , ( 4 ) where and represented stretch ratios of a sample from a storage protocol and a control sample, respectively, at 20 (i.e., different evenly spaced stress values along the stress - stretch curve. Tensile test . Q uasi - static uniaxial tensile tests were performed at a speed 20 mm/min on skin samples . Sandpaper w as used to increase the grip between the clamp and skin sample s . Four markers were placed on the central part of the skin sample ( Figure 12 d ) to track the variation in length using a CCD camera (Hitachi model KP - M2AN). A load cell ( Futek model LSB200 ) was used to record forces applied to the samples during the tensile test. Before each part of th e test, a preload of 10g (~ 0.1N) was applied to the sample and dimensions of the sample in different 27 reference configurations were recorded. The tensile test protocol consisted of three parts ( Figure 12 e ): (1) 10 loading - unloading preconditioning cycles at 20% strain; (2) a step - increase loading protocol, starting with two loading - unloading cycles at 20% strain; followed by a series of increasing strain levels (5% ~ 8% in crease) each performed for two loading - unloading cycles. The highest level of strain for each sample was repeated for an additional two loading - unloading cycles. The 5% ~ 8% strain range was estimated to accommodate the inter - specimen variation , so that similar peak stress es could be achieved for the last cycle. Peak stresses were determined during the te st by live readings of the loads and dimensions of the cross - section measured . (3) Rupture loading. Skin mechanical behavior has been previously estimated in a multitude of conditions, e.g., through cyclic loading - unloading test [74] [76] , through rupture test [77] [79] , and through in vivo test [80] [82] . In order to cover the mechanical response of skin under quasi - static lo ading comprehensively, we chose to investigate both preconditioned hyperelastic behavior and rupture behavior in this study. Part I of the loading protocol was to precondition the lower stress range of the skin mechanical response. Part II of the loading p rotocol was to obtain preconditioned hyperelastic mechanical behavior of skin that within the similar stress range (0 - ). Due to the fact that the uniaxial machine was strain - controlled, slight adjustment had to be made in strain increments for different samples to achieve the similar stress range for the last cycle in Part II. By the data from Part II, we were able to model the hyperelastic behavior of a sample and evaluate the corresponding constitutive parameters. Part III of the loading proto col enabled us to characterized skin stress - stretch behavior to its full capacity (i.e., until rupture) and evaluate the mechanical parameters associated with it. Due to the assumption of incompressibility, , where and are the sample leng th in the reference and current configuration, respectively; and are the sample cross - section al 28 area in the reference and current configuration, respectively. True stress (i.e., Cauchy stress) is calculated as , where is the for ce value in the current configuration measured by the load cell, and the stretch in the current configuration is measured by . The stress - stretch curve of each skin sample during the loading of the last cycle in Part II of the tensile test ( Figure 12 e ) was used for constitutive modeling purpose s due to its well preconditioned behavior, while the stress - stretch curve during the rupture loading was used to evaluate the mechanical parameters of skin sample s as describe d in the following section Mechanical parameter and statistical test . Mechanical parameters and statistical test . We evaluated five mechanical parameters of skin [83] , defined as follows and as shown in Figure 13 . I nitial slope , represent ed the stiffness of the rupture stress - stretch curve at the first 2% deformation strain. W e hypothesize d that this response was dominated by the ground matrix c ontent and exhibits approximately linear behavior. M aximum slope , describe d the stiffness of the quasi - linear part of the rupture stress - stretch curve at high stress level s , which we hypothesize d was dominated by the mechanical behavior of collagen. R upture stretch and U ltimate tensile strength (UTS) , represent ed the str etch and stress value at point of rupture , respectively. T oughness , represent ed the work input done per unit volume of a sample during rupture loading, defined as the area under the rupture stress - stretch loading curve. Upon completion of tensile test s, th e mechanical parameters of each sample , as described above , were evaluated. Three - way ANOVA was utilized to determine the influence of location, orientation, and sex on skin mechanical parameters. Significant level of statistical tests was set to be . 29 Figure 13 . Mechanical parameters on a representative rupture stress - stretch curve of a rat skin sample . 2.2.2. Histology and collagen network analysis Dog bone shape skin samples along the spine were obtained from another five male Sprague Dawley rats (10 to 12 weeks old) at the same four locations where samples for mechanical tensile testing were collected ( i.e., upper left, upper right, lower left, and lower right). Biopsies were excised parallel to the skin plane at the middle of the dermis layer. Two biopsies from each sample were stained using p Stained biopsies were examined using a Nikon mic roscope Eclipse 80i at a magnification. Nikon Digital Camera DXM 1200c was used to take digital images of the examined field. Measurement of relative collagen content. Relative collagen content was measured through rome stained histology slide [84], [85] . Images of total field area of were taken for each sample. It has been shown that collagen area percentage measured stained biopsies are highly correlated with collagen dry weight measured by a hydroxyproline assay [86] . A color - based segmentation algorithm from the MATLAB image processing toolbox was used to separate different colors from the images, i.e., 30 black (nuclei), red (cytoplasm, muscle, erythrocyte), and blue (collagen), a s shown in Figure 14 a - d . Different to other algorithms which require thresholds to define colors [87] , the LAB color space based algorithm allows the investigator to choose reference colors from the image and then categorize each pixel of the image into the closest reference color . This approach provided good agreement with visual observation s . Relative collagen content was calculated as the ratio between pixels identified as collagen and the total pixels identified as non - background. Finally, the relative collagen content of each back location was obtained by averaging measurements over animals. Automated detection of collagen fiber bundle orientation . PSR stained biopsies were examined under polarized light. Images of total field area were taken for each sample. Due to the birefringent nature of collagen fiber bundles, which was further enhanced by the PSR stain, thick, medium, and thin collagen fiber bundles exhibited red, yellow, and green colors under polarized light, respectively [87], [88] . Using the LAB color space based segmentation algorithm described above, we separated each PSR stained histology image into three ima ges ( Figure 14 e - h ) that each contained only thick fibers (red), medium fibers (yellow), and thin fibers (green). After collagen fiber bundles were grouped based on thick ness, a previously validated algorithm [38] was used to detect fibers orientation. Briefly, first, each color - segmented image was binarized by turning pixels representing collagen fiber bun dles into white pixels and pixels representing background into black pixels. Second , an erosion process was performed to detach cross - linked bundle from each other, and then a hole - filling process was performed to make each individual fiber bundle a solid recognizable piece. Then , each fiber bundle was fitted with an ellipse that satisfied the eccentricity and minimum size criterion, to exclude out - of - plane fibers from the analysis. An average of 16,372 collagen fiber bundles per sample were detected. T he orientation 31 of the major axis of each ellipse (positive in counter - clockwise direction) with respect to the rat dorsal midline (i.e., spine direction) was calculated. Finally, the total fiber orientation of each sample from the same location were averag ed over animals and the probability density of collagen fiber bundle orientation s were plotted for each location. Manual examination of collagen straightness . Collagen fiber bundles are undulated when skin is in a stress - free state [42] . Collagen undulation was m easured by introducing a straightness parameter, defined as , where represent ed the distance between two ends of each fiber, and represented the total length of each fiber ( Figure 14 i & j ) [89], [90] . 200 collagen fiber bundles in total (50 collagen fiber bundles per location) were measured across all 20 samples for straightness . Figure 14 . rat back skin histology slide, as shown in ( a ), were segmented into images with different components based on colors, i.e., ( b ) red (cytoplas m, muscle, erythrocyte), ( c ) blue (collagen), ( d ) black (nuclei). Collagen fiber bundles in a PSR stained rat back skin histology slide observed under polarized light , as shown in ( e ), were grouped based on colors, i.e., ( f ) red, ( g ) yellow, and ( h ) green, which indicated thick, medium, and thin fibers, respectively. Straightness of a single collagen fiber, as shown in ( i ), was measured in the following way as shown in ( j ). The g reen curve following the shape of the fiber bundle represented the orig inal length , while the white straight line connecting the two ends of the fiber bundle represented the compressed length . Straightness was calculated as . 2.2.3. Constitutive model We describe d the nonlinear, anisotropic behavior of s kin employing a mass - averaged strain energy function (SEF) of the form 32 , ( 5 ) where and represent ed the mass fraction and the SEF for each constituent, respectively, with for the ground matrix (elastin, glycosaminoglycan (GAG), proteoglycans and glycoproteins), and the collagen fiber family . A collagen family was defined as a collection of collagen fiber bundles oriented within a specified range, defined as , with respect to the dorsal midline. The mass fraction of a collagen family was defined as , where represent ed the number of collagen fiber bundles in the collagen family, and represented the total number of collagen fiber bundles. The total collagen mass fraction of a sample was determined by location - averaged values the ground matrix of a sample was calculated as . The relative collagen content measured by area percentage in histology images should be converted to mass fraction by a coefficient. However, the collagen area percentage in our study actually fell in the 60% - 80% range of collagen dry weight [33] . Therefore, the relative collagen content measure d by percentage area was implemented directly into the constitutive model for the mass fraction parameter . A dditionally, the skin was assumed to be an incompressible material, and plain strain conditions were assumed during uniaxial tensile test. The mechanical behavior of the g round matrix has been described by using a neo - Hookean constitutive model [50] , 33 , ( 6 ) where was a constitutive parameter with the dimension of a stress and represented the right Cauchy - Green deformation tensor for the ground matrix, defined as , where represent ed the deformation gradient for the ground matrix. The mechanical behavior of the collagen fiber network has been described by using a Fung - type exponential SEF [91] , , ( 7 ) where and were constitutive parameter s, the first with the dimension of a stress and the second being dimensionless, and describ ed the stretch in the direction of collagen fiber family defined as , where was the straightness parameter, and represent ed the direction of the collagen fiber family. The angle has been estimated by analyzing the histological images, and we considered the majority of the fibers to be contained within the plane of the skin [38] . Although in principle each family of fibers could show different mechanical behaviors, we assumed that all collagen families shared the same constitutive parameter s within each sample, namely , and , for . This simplification of the constitutive model helped t o provide a clearer interpretation of the contribution of the skin microstructure to the macroscopic level mechanical response of skin. Finally, and for each collagen family were estimated based on the location - specific collagen orientation distribution data ( Figure 18 a ) . To determine the correct number of collagen families, we conducted a parameter convergence study with the constitutive model. Results 34 showed that the value of and converged when . Therefore, 60 collagen families were adopted for the collagen fiber orientation distribution and constitutive parameter estimation. 2.2.4. Parameter estimation Constitutive model parameter estima tion was carried out using the nonlinear least - squares algorithm implemented in MATLAB (i.e., lsqnonlin function). The objective function employed in the optimization problem was , ( 8 ) where and represent ed the stress values at different stretch ratios in experimental data and constitutive model, respectively. Stress values at 30 evenly spaced stretch values were considered in the parameter fitting proc ess i.e. . Constitutive parameters , , , and were obtained through the parameter estimation process. The steps to obtain constitutive parameters for each sample were as follows. First, was estimated by fitting the stress - stretch curve at the initial stretch interval using only the neo - Hookean model, second, from the previous step was applied to the constitutive model, then , , and were obtained by fitting the entire stress - stretch curve of the sample using the complete constitutive model, see Equation ( 5 ) . Initial values for each parameters were randomly g enerated, between 0 to for , , and , and between 0.9 to 1 for . Initial values for each parameters were randomly generated, between 0 to for , , and , and between 0.9 to 1 for . 35 2.3. Results A total of 64 samples were successfully tested. Upon inspection of experimental data, the stress - was therefore excluded from the analysis of rupture data. 2.3.1. Storage condition effect The five mechanical parameters of Group I samples were evaluated and grouped based on storage protocol ( Figure 15 ). A one - way ANOVA showed that storage protocol did not have a significant effect ( on initial slope, maximum slope, UTS, rupture strain, or toughness. Additionally, averaged cyclic loading stress - stretch curves and ruptur e loading stress - stretch curves of Group I samples from each storage protocol were shown in Fig ure 16 . NRMSD between the control group (fresh samples) and each of the oth er three storage protocols were calculated to further compare the influence of the storage protocol on the nonlinear mechanical behavior of skin ( Table 1 . The NRMSD between each storage protocol and the control group for cyclic loading stress - stretch curves and rupture loading stress - stretch curves. were calculated . Results showed that Protocol II had the lowest total NRMSD (7.11%), followed by Protocol III (13.75%) and Protocol I (16.97%). Since the influence of different storage conditions on the mechanical properties of skin was not significant, samples from Group I and Group II animals were put together to study the infl uence of location, orientation, and sex on skin mechanical properties. 36 Figure 15 . Influence of storage condition on initial slope, maximum slope, UTS , rupture strain, and toughness. Values shown in the bar graph are means and standard deviations. Values shown below each bar graph are p - values of one - way ANOVA comparing the mechanical parameter difference between four different storage conditions . Fig ure 16 ( a ) Averaged cyclic loading curves and ( b ) averaged rupture loading curves of each storage protocol. Horizontal bars on each curve represent standard errors . Table 1 . The NRMSD between each storage protocol and the control group for cyclic loading stress - stretch curves and rupture loading stress - stretch curves. C yclic loading R upture loading Total Protocol I 13.21% 3.77% 16.97% Protocol II 4.67% 2.44% 7.11% Protocol III 5.94% 7.81% 13.75% 37 2.3.2. L o cation, orientation and sex effect Mechanical parameters evaluated from rupture stress - stretch curves were grouped based on three different factors, i.e., location (upper vs. lower back), orientation (axial vs. transverse direction), and sex (male vs. female). Group means and standard devia tions based on different factors were shown in bar graphs in Figure 17 . Influence of location, orientation and sex on ( a ) initial slope , ( b ) maximum slope , ( c ) rupture strain , ( d ) UTS, and ( e ) toughness. M eans and standard deviations are shown in the bar graphs . P - values of the three - way ANOVA are shown below the bar graph for comparison within each factor, i.e., upper vs. lower, axial vs. transverse, male vs. female. Statistical significance ( <0.05) and significance ( <0.001) are single and double asterisks, respectively . P - values of the three - way ANOVA showing the effect of different factors on mechanical paramete rs were listed below the corresponding paired bar graphs. Factors that contributed statistical significance were marked using asterisks. No significant differences were found between sexes on mechanical parameters, except for sample thickness (male: , female: , ). However, significant differences were detected between locations, and orientations. The lower back had significantly higher values of initial slope, maximum slope, UTS, and toughness when compared to t he upper back. Orientation - wise, significant differences were detected between axial and transverse directions for the initial slope, maximum slope, rupture strain, and toughness. The axial direction had higher values of initial slope, rupture strain, and toughness. While transverse samples had higher values of maximum slope. 38 Figure 17 . Influence of location, orientation and sex on ( a ) initial slope, ( b ) maximum slope, ( c ) rupture strain, ( d ) UTS, and ( e ) toughness. M eans and standard deviations are shown in the bar graphs . P - values of the three - way ANOVA are shown below the bar graph for comparison within each factor, i.e., upper vs. lower, axial vs. transverse, male vs. female. Statistical significance ( and strong st atistical significance ( are marked by single and double asterisks, respectively . 2.3.3. Fiber s orientation, straightness and relative collagen content Orientation probability density of collagen fiber bundles for upper right back location were shown in Figure 18 a . Each bin represented a collagen family, where the width of the bin was the angle ran ge that a collagen family covered ( i.e., , and the height of a bin was the probability densit y of that collagen family. Probability density was calculated as , where was the number of fibers in the collagen family and was the total fiber numbers detected . As to the pattern of fiber orientation distribution, collagen fib er bundles aligned predominantly in two directions , 40° to 50° and - 40° to - 50° , with respect to the spine at all four locations . The collagen fiber bundles straightness for each location was shown in Figure 18 b . No significant difference was found between locations. All back locations had mean collagen fiber straightness equal to ~0.95. Relative collagen content was shown for upper and lower back in Figure 18 c . A Student 39 t - test showed that samples from the lower back have significantly higher relative collagen content when compared to samples from the upper back ( , with . Figure 18 . ( a ) upper left location on rat back . ( b ) Means and standard deviations of collagen fiber straightness at four back locations measured in PSR stained histology images . 50 fibers/5samples were measure d for each locations . UL upper left, UR upper right, LL lower left, and LR lower right. ( c ) Means and standard deviations of relative collagen content at upper and lower back ) is marked using asterisk . 40 2.3.4. Constitutive model parameters Representative experimental data and constitutive mode l with corresponding best - fit constitutive parameter s of back skin samples from 2 animals at upper and lower locations, and in axial and transverse directions were shown in Figure 19 . NRMSD ( for 64 samples) showed a good agreement between experimental data and constitutive model. A uniqueness study showed that the best - fit constitutive parameters of each sample were independent by the initial guess values. Best - fit constitutive parameters obtained from the parameter estimation process were then grouped based on three factors, i.e., location, orientation, and sex ( Table 2 ) . The influence of different factors were evaluated using a three - way ANOVA ( Figure 20 ). No significan t difference was found between sexes. Significant differences were detected between locations, as well as orientations. Upper back samples had significantly lower values of and , but significantly higher values of when compared to lower back samples. With respect to the orientation, samples tested in the axial direction had significantly higher values of and , and significantly lower values of when compared to the transverse direction. 41 Figure 19 . Rep resentative experimental data and constitutive model with corresponding best - fit constitutive parameter s at upper and lower rat back locations . Axial experimental data are from animal 5, while transverse experimental data are from animal 10. Best - fit constitutive parameters of the simulated stress - stretch curves shown in the figure are as follows. Upper back sample in axial direction (animal 5) , , , ; and transverse direction (animal 10) , , , . Lower back sample in axial direction (animal 5) , , , ; and transverse direction (animal 10) , , , . Table 2 . Means and standard deviations of best - fit constitutive parameter s , , , and , grouped based on lo cation, orientation, and sex. upper lower axial transverse male female 42 Figure 20 . Influence of location, orientation, and sex on best - fit constitutive parameters, ( a ) , ( b ) , ( c ) , ( d ) . M ean s and standard deviation s are shown in the bar graphs . P - values of a three - way ANOVA are shown below the bar graph for comparison within each factor , i.e., upper vs. lower, axial vs. transverse, male vs. female. Statistical sign ificance ( and strong statistical significance ( are marked by single and double asterisks, respectively . 2.4. Discussion and conclusions Differences in the mechanical properties of human skin associated with sex have been reported in literatur e [82], [92] . Similar studies on murine skin were rare [93] . Considering the wide application of rat skin as an alternative model for human skin [94] [97] , it is im portant to investigate the influence of sex on the mechanical properties of rat skin. Results of our study showed that sex ha s no significant effect on any of the evaluated parameters of rat back skin. However, male animals showed significantly higher skin thickness when compared to females , which confirmed what had been observed in humans previously [98] . When we investigated the effect of location on the mechanical parameters, w e found that the lower back had higher values 43 o f initial slope, maximum slope, UTS, and toughness when c ompared to the upper back . Tensile test data collected by Haut [99] showed that the lower back h ad higher UTS compared to the upper back o f 2 month - old rats, which agreed with our study. However, our results did not confirm the rupture strain difference between the upper and lower back that were reported in the same study [99] . When it comes to anisotropy, Karimi, et al. [76] reported that the transverse direction of mice back skin had higher maximum slo pe compared to the axial direction, which was supported by our results. Our results also confirm ed the higher failure strain of samples tested in the axial direction compared to the transverse direction reported by Haut on 2 month - old rats [99] . However, our findings did not confirm the UTS difference between the axial and transverse direction of 2 m onth - old rats reported by Haut. Although there is no clear consensus in the literature regarding the influence of location and orientation on the mechanical parameters of rat skin, which may be due to different experimental conditions (preconditioning, tes ting rate, animal age etc.), the authors are confident in stating that heterogeneity and anisotropy are crucial characteristics of skin mechanical properties and should be addressed when testing this tissue. One limitation of the tests conducted here was t hat we kept the samples hydrated throughout the test by spraying saline solution (i.e., HBSS) directly on them, at 5 minutes intervals. This may however alter the tissue mechanical behaviors compared to the in vivo condition. Moreover, differences between hydration conditions could be one of the factors that contribute d to the discrepancies reported in literature. Numerous studies have investigated the storage condition effect on skin tissue cell survival rates [100] [102] for the purpose of skin graft storage in clinical practice . Sample storage is also usually required for skin biomechanics study due t o the time interval between sample collecting and sample testing. However, very few studies have reported the storage condition effect on mechanical properties of skin. Marangoni et al. [103] recorded the influence of three storage 44 conditions , i.e. (1) 1.6 °C, (2) 37.8 °C, and (3) - 92.8 °C preceded by flash freezing, on mechanical properties of skin samples collected from guinea pigs during the time span of 6 hours. Their study concluded that - 92.8 °C storage condition following flash freezing maintained values of maximum stress, maximum strain, maximum stiffness and maximum work input of the skin samples, while 1.6 °C and 37.8 °C significantly altered the values of those evaluated parameters. Foutz et al. [104] reported that strength, ultimate strength . Although conclusions from previous studies are not consistent, it can be generalized that low temperature storage preceded by flash freezing may maintain mechanical properties of skin tissues. Our study compared the influence of three different storage protocols on skin mechanical parameters ( Figure 15 ), and we found no significant difference between stored samples and fresh samples. When comparing averaged stress - stretch curves of different storage conditions to those of the controls, we f ound that Protocol II (i.e. flash freezing, then stored in - 78.5 °C, and thawed at 4 °C for 6 hours before testing) led to the smallest NRMSD among the three protocols. The structure of the collagen network in the dermis is the link that connects the macr oscopic level mechanical response of skin to its microscopic structure, and the quantification of its characteristics remains a challenging task. Specifically, collagen fiber bundles orientation distribution, relative collagen content (either mass or volum e), and collagen fiber undulation, have been previously shown to play an important role in defining a successful structural constitutive modeling framework [43], [50], [105] . For example, Jor et al. [37], [49] quantified collagen orientation distribution in the cross - section of porcine abdomen skin, and N í Annaidh et al. [38] measured collagen orientation distribution in the parallel section of human back skin. In the current study, we developed an autom ated image analysis routine based on [38] . In order to improve fiber 45 identification, we used PSR hist ological stain ing instead of Van Gieson . One reason for this choice is that collagen fiber bundles can be more easily identified in PSR stain ed tissue when imaged under polarized light , because of the absence of other constituents in the image, which are i nstead still visible with Van Gieson stains. Another reason is that f ibers with different thickness es exhibit different colors in a PSR stain under polarized light, which effectively reduce s fiber cross - linking and increase s the efficiency of recognizing i ndividual collagen fiber bundles, by group ing fiber s through color segmentation. Also, N í Annaidh et al. [38] isolated collagen fiber bundles by method, which could misidentify some collagen pixels as background. In this study, collagen fiber bundles of different thicknesses were first isolated by their color ( Figure 14 f - h ), then those exact collagen pixels were turned into binary image s . This method decrease d the amount of collagen pixels los t during image processing. Our results also showed that in - plane collagen fiber bundles were distributed mainly in two directions, ±40° ~ 50° with respect to spine, although a certain percentage of fibers w ere distributed equally in every direction. A similar distribution of collagen fiber bundles has been reported previo usly in humans [38] and pigs [37] . On one hand, we found no difference in collagen orientation distri bution ( Figure 18 a ) or straightness ( Figure 18 b ) between the four locations. On the other hand, we observed a significantly higher relative collagen content in the lower back when compared to the upper back ( Figure 18 c ) , which could be a contributing factor to the difference in mechanical parameters observed here ( Figure 17 ). In the micro structural constitutive model presented, the collagen orientation distribution was described by a number of discrete collagen families, which led to a closed - form analytical solution for a collagen fiber bundle. Similar approach es in the literature include the models from Flynn et al. [105] , Karimi, et al [41] , and Groves, et al. [40] . The use of d iscrete collagen families avoids 46 the sometimes costly numerical integra tion compared to models that employ continuous fiber orientation distribution [38], [49] . The mass fraction not only described the distribution characteristics of collagen fa milies, but weighed the contribution of ground matrix and collagen. The location - specific measurement of relative collagen content , along with the fiber bundles orientation distributio n measurement, ensured that best - fit constitutive parameters could accou nt for the real contribution of each constituent in skin to its macroscopic level mechanical response. Therefore, constitutive parameters from different samples can be further compared to investigate the influence from external factors. The influence of lo cation and orientation detected on the mechanical parameters ( Figure 17 ) were also reflected in the constitutive parameters ( Figure 20 ), which cover ed both high - stress and lower - stress range s of skin mechanical properties, respectively. In our model, the parameter values of ( Figure 20 a ) were directly related to values of initial slope ( Figure 17 a ) since we hypothesize d that the ground ma trix dominated the response in lower stress range in our constitutive model. It should be noticed that the standard deviations of initial slope and were large. This could be due to experimental errors at the early phase of the testing, whose force r esponse were relatively small. The i nfluence of location and orientation detected in best - fit values differed from the manually measured results. Manually measured values ( Figure 18 b ) show ed no location difference, while best - fit values from the constitutive model ( Figure 20 d ) show ed significant difference between the upper and lower back, and also significant difference between axial and transverse directions. This discrepancy could come from the following reasons. (1) Limited number of manually measured fiber bundles: only 10 fibers per sample were manually measured, while the best - fit values reflected the averaged straightness of all fiber bundles within a sample. (2) Different reference configuration : manually measured collagen fiber bundles were found in deformation - fre e samples, while best - fit values for the 47 constitutive model were for samples that had been through preconditioning and cyclic - loading. (3) Inter - specimen variation . However, the averaged value of manually measured had good agreement with that obt ained from the parameter estimation process, being 0.95 and 0.96, respectively. 48 A SUBJECT SPECIFIC MODEL THAT ADDRESSES NONLINEAR MECHANICAL BEHAVIOR OF HUMAN THIGH INFORMED BY IN VIVO EXPERIMENTAL DATA 3.1. Overview The goal of this chapter is to develop a subject - specific FE model that addresses the nonlinear mechanical behavior of thigh soft tissue. A parameter estimation process that matched nonlinear force - deflection behavior of FE model to in vivo data was developed. Effects of modeling region dimension (3D vs semi - 3D) and constitutive model (neo - Hookean, Mooney - Rivlin, first Order Ogden, and Fung orthotropic) were studied using a subject - specific thigh FE model. After the best two constitutive model s were determined, a total of 20 thigh FE models were developed and corresponding best - fit constitutive parameter s were obtained using the two constitutive model s. Effects of sex and location on thigh soft tissue were also detected on those thigh mechanical parameters. 3.2. M ethods 3.2.1. Mechanical testi ng We recorded force - deflection data sets from the leg for 20 individuals : 10 males (average age: 20.6 years old , average weight: 76.2 kg, average height: 170.4 cm) and 10 females (average age: 20.9 years old , average weight: 66.4 kg, average height: 162.8 cm), following a previously published protocol [106] . Briefly, we first recorded standard anthropom etric measurements, such as height, weight, seated height, buttocks width , and leg dimensions for each participant. Each participant was seated on the ischial tuberosity so that the soft tissue of the thigh was unloaded and undeformed. Then, a custom made hand - held compression device was used to gather force - 49 deflection data on the proximal, middle, and distal locations on the posterior side of the thigh. Load was applied until a biological barrier was reached. 3.2.2. Geometry generation The participant specific g eometry was recorded for one representative participant (female) employing the following protocol. The participant was seated on the ischial tuberosity so that the thigh was not deformed and the knee was flexed at . A Sense TM V2 3D scanner was used to record the CAD model of the thigh ( Figure 21 b ). A femur geometry obtained from GrabCAD [107] was placed within the model of the thigh, and we use d anthropometric measurements of the participant and CT images of thigh cross - sections from literature [108] to ensure the correct size and the accurate position of the femur. The geometry of the remaining 19 participants was obtained by employing a scaling technique starting from t he geometry of the one representative participant. Specifically, the singular thigh CAD model developed was scaled to that of another participant by matching two key participant - specific measurements collected at the time of mechanical testing, namely the length from hip joint to knee joint and the circumference of the middle thigh. All participants weight and height fell into the US mid - size male and female groups previously published [109] . This enabled us to use a constant diameter for the femur within the same sex group, while the diameter of the femur of male participants was increased by 4.8% with respect to that of female participants [110] . 3.2.3. Finite element model Mesh generation . Each CAD model was cut along two transverse planes, one passing through the joint between buttock and thigh and ano ther passing through the joint between thigh and knee ( Figure 21 c ). Then, the model was meshed (HyperMesh ver.13.0) using 10 node s tetrahedral elements of size 10 mm. Th e FE models had an average of 32 , 062 elements and 47 , 518 nodes. A 50 of size 10 mm (i.e., 27,724 elements, 41,227 nodes), elements of size 7.5 mm (i.e., 54,568 elements, 79,830 nodes), and elements of size 5 mm (i.e., 141,083 elements, 203,742 nodes). We compared the normalized root - mean - sqaure deviations (NRMSDs , E q uation ( 17 ) ) between force - deflection curves of the FE model with 10 mm - sized elements and the FE model with 7.5 mm and 5 mm - sized elements, at each tested location. All NRMSDs are b e low 2.2%. This confirms that mesh convergence was achi eved using elements of size 10 mm. The soft tissue in the model was then divided into three locations, corresponding to the three testing sites. The same homogenous constitutive model but different constitutive parameter s described the lumped mechanical pr operties of soft tissues (i.e., skin, fat, and muscle) at each location, and the femur was modeled as a rigid body. 3D model . Three areas were identified at the proximal, middle, and distal locations on the posterior side of the thigh ( Figure 21 d , top row), each corresponding to a location tested during the in vivo experiment. Boundary conditions were applied to the femur (i.e., constrained against disp lacement and rotation in any direction), to the two cut planes at proximal and distal thigh (i.e., nodes were constrained against longitudinal displacements along the femur longitudinal direction), and to the three testing areas (i.e., nodes were restricte d to move only in an anterior - posterior direction, to prevent instability during deformation). Additionally, participant - specific compression forces at each location, as recorded in vivo , were applied to one testing area at a time. Compression forces were applied in 20 evenly spaced steps from zero to the maximum value recorded in the experiment, and the nodal displacements of the compressed area 51 at each step were recorded and averaged. Simulated force - deflection curves were then compared to the experimenta l data to estimate the best - fit constitutive parameter s . Figure 21 . ( a ) Schematic of the seating positioning for a participant during 3D scanning. ( b ) Thigh CAD model obtained from 3D scanning of the thigh of a representative participant. ( c ) Representative 3D FE model. ( d ) Dimensions and positioning of testing areas for the 3D FE model, on top, and the semi - 3D FE model, on the bottom. From left to right, proximal, middle, and distal thigh testing locations are shown in red. ( e ) Flow chart of the constitutive parameter optimization protocol . Semi - 3D model . Additionally, we developed three semi - 3D FE models of one representative participant. We cut the thigh CAD model at the proximal, middle, and distal tested areas respectively, using two transverse planes placed at from one other. These three long (in the femur longitudinal direction) CAD models were then meshed using the same element type and element size as the 3D FE model described above ( Figure 21 d , bottom row). The boundary conditions applied to the femur, cut planes, and tested areas, the load application procedure, and the constitutive parameter estimation were the same a s discussed above for the full 52 3D model. We then compared the best - fit constitutive parameter s of 3D and semi - 3D FE model s that described the same set of experimental data , as well as the simulated force - deflection data of 3D and semi - 3D FE model estimated employing the same set of constitutive parameter s. 3.2.4. Constitutive model s and parameter optimization Neo - Hookean model . The strain energy function for the neo - Hookean model [111] was defined as , ( 9 ) where was a constitutive parameter with the dimension of a stress, and were respectively, the first and third invariant of the right Cauchy - Green strain tensor, , where was the deformation gradient. In the principal reference system, one could write , where and represent ed the s tretches in the principal directions. The coefficient represented the bulk modulus - like penalty parameter, which was defined as , where represented [20] to ensure a nearly incompressible material description. Mooney - Rivlin model . The strain energy function for the Mooney - Rivlin model [112] was defined as , ( 10 ) where and were constitutive parameter s with the dimension of a stress, and were the first and third invariant of , as defined above, and was the 53 second invariant of the tensor . The bulk modulus - like penalty parameter in this case was defined as . First order Ogden model . The strain energy function for the first order Ogden model [113] was defined as , ( 11 ) where and were constitutive parameter s with the dimension of a stress and dimensionless, respectively, was the third invariant of and and represent ed the stretches in the principal directions, as defined above. The bulk modulus - like penalty parameter in this case was defined as . We defined also the parameter initial stiffness as , ( 12 ) which corresponded sticity theory, i.e., represented the slope of the stress - strain curve under uniaxial loading when the material underwent small deformation. Fung orthotropic model . With the assumption of isotropic symmetry, the strain energy function of the Fung orthotropic model [114] was defined as , ( 13 ) where was a constitutive parameter with the dimension of a stress, was the third invariant of , and was the bulk modulus - like penalty parameter, where represented the we defined the exponential parameter as 54 , ( 14 ) where w as the Green strain tensor, was the identity tensor, and the dimensionless constitutive parameter was defined as . ( 15 ) Constitutive parameter optimization . Constitutive parameter s at proximal, middle, and distal thigh locations were optimized employing the simplex search method implemented in MATLAB (i.e., function fminsearch) . The FE simulation in each optimization iteration was performed in FEBio. The flow chart of the optimizat ion process is shown in Figure 21 e . The goal of the optimization process was to minimize the difference between the FE simulated force - deflection curves and the experimen tal data at each location. The objective function to be minimized was , ( 16 ) where , and were normalized root - mean - square deviation (NRMSD) estimated at each location as ( 17 ) where was the number of loading steps (i.e., ), and were the deflection values for the location from the FE simulation and experimental measurement respectively, and represented the maximum experimental values of deflection for each location (i.e., the deflection measure for the maximum load). The initial guess co nstitutive parameter s for each constitutive model were obtained from literature, as shown in Table 3 , and all material parameter values were constraint to be positive. 55 Table 3 . Initial guess constitutive parameter s for each constitutive model . Constitutive model Constitutive parameter Initial guess parameter value Reference neo - Hookean [ kPa ] 8.50 Linder - Ganz et al. [19] Mooney - Rivlin [ kPa ] 1.65 Ver ver et al. [21] [ kPa ] 3.35 Ogden [ kPa ] 3 Oomens et al. [16] 30 Fung orthotropic [ kPa ] 10 Pilot investigation of current study 3.2.5. Statistical analysis After the best - fit constitutive parameter s of all 20 participants were obtained, statistical analysis was carried out to determine the influence of sex and location on constitutive parameter s. - test was utilized to determine the sex influence on constitutive parameter s at each thigh lo cation, while one - way repeated measures ANOVA, followed by Holm - Sidak post - hoc test, was utilized to determine the influence of location on constitutive parameter s within each sex group. 3.3. R esults 3.3.1. Effects of the choice of constitutive model and of modeling r egion dimension All results in this s ection have been evaluated considering the geometry and mechanical properties collected experimentally for one representative participant. Effect of the choice of constitutive model . Constitutive model effects were investigated using a 3D FE model. Figure 22 showed experimentally measured (symbols) and numerically calculated (lines) compression force - deflection curves f or each thigh location and for each constitutive model considered. All sections of Figure 22 a - d show ed the same three sets of experimental data collected in vivo . Each section of Figure 22 also showed the numerical results for a 3D FE model, each of these models had the same geometry, and the soft tissue was described employing one of t he four constitutive model s . The numerical results shown have been obtained employing the 56 optimized constitutive parameter s, and the results were presented for each testing location, namely proximal (solid), middle (dashed), and distal thigh (dotted). Values of best - fit constitutive parameter s for each constitutive model and each location were also reported in Table 4 . The NRMSD, Equation ( 17 ) , between the 3D FE simulation and the experimental data was calculated fo r each tested location and for each constitutive model adopted, in addi tion to the total NRMSD, Equation ( 16 ) . The neo - Hookean model had the largest total NRMSD (35%), followed by the Mooney - Rivlin model (26%), the Fung model (18%), and finally the first order Ogden model (8%). In the analysis of the datasets of the overall 20 participants we employed fun orthotropic and the Ogden material models. Figure 22 . Compression force - deflection curves from experimental measures (symbols) and model predictions (lines) for each location and each constitutive model considered. ( a - d ) The same experimental data are shown in every section, repre senting the mechanical behavior for representative participant F6 at the proximal (triangles), middle (circles), and distal (asterisks) locations. Also shown are the model predictions after parameter optimization for each location, namely proximal (solid l ine), middle (dashed line), and distal (dotted line). The models considered are ( a ) neo - Hookean, ( b ) Mooney - Rivlin, ( c ) Ogden, and ( d ) Fung orthotropic . Effect of the choice of modeling region dimension . First, we compared best - fit constitutive parameter s optimized for 3D and semi - 3D FE mode ls for the three thigh regions, employing the same experimental data set. Percentage differences between the semi - 3D and 3D models were 57 shown in Table 4 . The largest difference occurred at the middle thigh location for each constitutive model : 84% for neo - Hookean, 253.2% for Mooney - Rivlin, 385.7% for Ogden, and 99.3% for Fung. We then compared force - deflection mechanical behavior of the 3D and semi - 3D FE models employing the best - fit constitutive parameter s evaluated for the 3D FE model shown in Table 4 . Figure 23 shows deflection differences between the semi - 3D FE model and 3D FE model (i.e., ) for each constitutive model at each location ( Figure 23 a - c , respectively). The NRMSD between the semi - 3D FE model and the 3D FE model at each location for each constitutive model were also calculated. A total value of NRMSD is reported here, estimated as the sum of the values at each location. The Mooney - Rivlin constitutive model showed the largest total NRMSD (81%), followed by the neo - Hookean (64%), Ogden (38%), and Fung (29%). It is important to notice that the effe ct of the modeling region dimension choice (semi - 3D vs. 3D) affects the predicted behavior more dramatically than the choice of constitutive model , based on the value of total NRMSD. Table 4 . Best - fit constitutive parameter s at proximal, middle, and distal thigh location for each constitutive model (i.e., neo - Hookean, Money - Rivlin, Ogden, Fung orthotropic). Experimental data collected for representative subject F6 have been used to inform the 3D and semi - 3D FE model. The perc entage difference of best - fit constitutive parameter values indicate semi - 3D model have higher parameter values, and negative values indicate 3D mo del have higher parameter values). Constitutive model Constitutive parameter Proximal thigh Middle thigh Distal thigh 3D Semi - 3D % 3D Semi - 3D % 3D Semi - 3D % neo - Hookean [ kPa ] 8.50 11.90 40.0 3.12 5.74 84.0 3.02 3.61 19.5 Mooney - Rivlin [ kPa ] 2.18 2.22 1.8 1.49 1.74 16.8 1.73 1.84 6.4 [ kPa ] 6.00 9.69 61.5 1.11 3.92 253.2 0.62 1.84 196.8 Og den [ kPa ] 0.49 2.30 369.4 0.49 2.38 385.7 1.70 3.43 101.8 35.78 12.97 - 63.8 11.85 6.00 - 49.4 5.38 3.50 - 34.9 Fu ng orthotropic [ kPa ] 25.30 18.25 - 27.9 3.00 5.98 99.3 2.70 3.22 19.3 [ kPa ] 3.20 0.95 - 70.3 16.70 13.41 - 19.7 23.30 16.31 - 30.0 58 Figure 23 . Mechanical behavior comparison between 3D and semi - 3D FE thigh models of participant F6 at ( a ) proximal thigh, ( b ) middle thigh, and ( c ) distal thigh. Horizontal axis represents forces applied to the soft tissue of thigh at each loading step, while the vertical axis represents deflection difference between semi - 3D FE model and 3D FE model at each loading step, i.e. , . 3.3.2. Ef fect of sex and location Best - fit constitutive parameter s for Ogden and Fung models were obtained for all 20 participants, employing fully 3D geometries scaled to match the anatomical measures of each participant. Figure 24 showed box plots for the constitutive parameter s for both models, grouped by sex (i.e., male and female) and location (i.e., proximal, middle, and distal thigh). Specifically, Figure 24 showed the parameters , , and for the Ogden model on the left, where was calculated employing Equation ( 12 ) . The parameters , , and for the Fung orthotropic model are shown on the right, where was calculated employing Equation ( 15 ) . Effect of sex . The constitutive parameter s for each model were compared between males and females at each location. The results of the statistical analysis were presented in Table 5 . For the - test found that males had a significantly higher value of at the proximal thigh location ( p < 0.05), as well as a higher value of at the pr oximal and middle thigh locations ( p < 0.05 for both locations) when compared to females. For the Fung orthotropic model, males had a higher value of at the middle thigh location ( p < 0.05), and a higher value of at the proximal thigh location ( p < 0.05). 59 Effect of location . Since we detected multiple differences in constitutive parameter s between male and female participants, we analyzed the effect of location on the constitutive parameter s within the same sex. One - way repeated measures ANOVA f ound that location had a significant effect on (1) and for the Ogden model ( p < 0.001 for both parameters and for both male and female group) and (2) all constitutive parameter s , and for the Fung orthotropic model ( p < 0.001 for both male and female group for all parameters). To better understand which location comparisons contribute d to the differences in the constitutive parameter s, we carried out a Holm - Sidak post - hoc pairwise test for the following locations : proxima l vs. middle, proximal vs. distal, and middle vs. distal. The results of the statistical analysis were presented in Error! Reference source not found. . Table 5 . P - - test comparing sex difference on Ogden and Fung orthotropic constitutive parameter s at different thigh locations. The gray areas highlight the statistical significance, i.e., p < 0.05. Male vs Female Proximal thigh Middle thigh Distal thigh 0.038 0.429 0.865 0.682 0.235 0.547 0.023 0.021 0.846 0.170 0.502 0.223 0.035 0.252 0.105 0.101 0.004 0.297 Table 6 . P - values of Holm - Sidak post - hoc pairwise tests comparing model parameter difference between different thigh locations , in male and female subjects group, for both Ogden and Fung constitutive models, respectively. The gray areas highlight the statistical significance, a lighter gray corresponds to statistical significance ( p < 0.05), and a darker gray corresponds to strong statistical significance ( p < 0.001). M ale F emale P rx vs M id Prx vs Dis Mid vs Dis P rx vs M id Prx vs Dis Mid vs Dis <0.001 <0.001 <0.001 <0.001 <0.001 <0.001 <0.001 <0.001 <0.001 <0.001 <0.001 0.052 <0.001 <0.001 0.020 <0.001 0.001 0.050 <0.001 <0.001 0.186 <0.001 <0.001 0.017 <0.001 <0.001 <0.001 <0.001 <0.001 0.001 60 Figure 24 . Best - fit Ogden model and Fung orthotropic model parameters, i.e., ( a ) , ( b ) , ( c ) , ( d ) , ( e ) , and ( f ) for 20 tested subjects, at proximal, middle, and distal thigh locations, respectively. Parameter values are group ed based on location (marked at the top of each plot) and sex (marked at the bottom of each plot). The central mark within the box indicates the group median, and the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively. The whiskers extend to the most extreme data points not considered outliers, and the outliers are plotted individually using the '+' symbol. 3.4. D iscussion and conclusions Nonlinearity is a unique characteristic of soft biological tissues and failing to include t hat within models could significantly affect the accuracy of the predictions. However , to be able to capture non - linear behavior in large deformations, one needs to have access to experimental data that (1) record both forces and deformations simultaneousl y and (2) expand over an appropriate loading range. In other words, there is a need for data from a mechanical test that captures the in 61 vivo mechanical behavior of soft tissues under large deformations, in an appropriate full - body configuration (i.e., sea ted vs. laying down). Such data sets have not been available until recently [106] . Previously published data still made an effort to describe the nonlinearity of the thigh and buttocks soft tissues. For example, some studies obtained constitutive parameter s from literature and later validated the FE model by reporting the error between the FE simulation and the experimental data [19] [21] . Other studies performed an optimization to estimate the best - fit constitutive parameter s that minimized the error between the FE simulation and the experimental data in one specific configuration [27], [28] . The experimental data available for both approaches were , however , limited and described the deformations corresponding to one loading configuration (i.e., weig ht - bearing configuration in upright sit t ing position). Due to the non - linear behavior of soft tissues, we propose that to have access to data in only two configurations, namely an unloaded position and a loaded seated position, is not sufficient to describ e accurately the mechanical behavior of the tissues. This is a pressing issue, especially when determining the stress distributions within the tissue that cannot be directly experimentally validated. The optimization process in this study used the simplex search algorithm ( fminsearch function in MATLAB). A major reason of using this algorithm was that it was a derivative - free method. Algorithms like trust - region would require gradient information of the objective function. The objective function in our stu dy had no analytic expression and perform like a black box: inputs were the constitutive parameter s and outputs were the NRMSD s between the FE simulation and experimental data. Due to the lack of analytic expression s , the convexity of the objective functio n cannot be confirmed. Therefore the optimized constitutive parameter s could be solutions to local minima instead of global minima. However, the optimized constitutive parameter s were of the same magnitude as those in literature ( Table 3 ) and provid e d good match to experimental data. 62 They could be seen as a subset of reasonable solutions, if multiple solutions exist ed . It is important to notice that fminsearch does not provide constraints for the optimized parameters . Therefore, constraints were added in the optimization process to give artificial high error values for parameter searching attempts that fell out of constrained area ( i.e., constitutive parameter s being positive ) . The optimization procedure in this study compared deflection estimated by FE simulation to deflection from a unique experimental dataset collected in vivo for a value of loading that gradually increased from 0 to a maximum load [106] . Because of the continuous nature of the experimental dataset, this work had the capability of accurately describing the mechanical behavior of the thigh soft tissues for a large range of loading configurations . For this reason, we had the opportunity to compare the accuracy of four different materials models in describing the behavior of soft tissues. We show ed that 3D FE models employing the neo - Hookean and the Mooney - Rivlin material models (i.e., NRMSD = 35% and 26% respectively) describe d the thigh soft tissue mechanical behavior of the representative subject less accurate than Fung orthotropic and first orde r Ogden material models (i.e., NRMSD = 18% and 8% respectively). As to the form of discrepancy between FE simulation and in vivo data, we observe d that neo - Hookean and Mooney - Rivlin material model largely underestimate d deformations at the low force level and overestimate d deformations at the high force level ( Figure 22 a - b). In other words, neo - Hookean and Mooney - Rivlin material models were - linear orthotropic model overestimate d deformations at the low force level and underestimate d deformations at the high force level ( Figure 22 d). In other words, Fung orthotropic material model was - rder Ogden material model offered the lowest error ( Figure 22 c). The strain energy function formulation of each material model also matched the mechanical behavior characteristics we 63 observe d above. Being a polynomial function with exponents values ( ) adjustable, first order Ogden material model can describe non - linear mechanical behavior in a broad range. Being special cases of Ogden material model (i.e., neo - Hookean is a special case of first order Ogden material model, Mooney - Rivlin is a special case of second order Ogden material model), the exponents values of neo - Hookean and Mooney - Rivlin material models are fixed and therefore the non - linear mechanical behavior they can describe are largely limited. Similarly, the mechanical behaviors that Fung orthotropic material model can accurately describe are limited by the natural exponential function. With the results of first order Ogden model being more accurate than Fung orthotropic mat erial, it suggests that the force - deflection curves in this dataset show a more polynomial than natural exponential behavior, overall. After examining the force - deflection behavior characteristics of data from the remaining 19 subjects, we determine d the n eo - Hookean and Mooney - Rivlin material model were not suitable to describe the highly non - linear in vivo data in this study. Therefore, Fung and first order Ogden model were employed for the remaining 19 FE models. It is important to notice that neo - Hookean and Mooney - Rivlin material models are still commonly used in soft tissue modeling [115], [116] , our results however suggest that these material model s should be used with care, along with proper experimental data, parameter optimizat ion process, and detailed micro structural information [117] . In this study we also quantified the differences between 3D FE model and semi - 3D FE model in two ways: (1) the constitutive parameter differences of the two types of model when f itting the same set of experimental data and (2) the force - deflection mechanical behavior difference when employing the same set of constitutive parameter s. Regarding the first approach, as shown in Table 4 , semi - 3D models, when compared to 3D models, could overestimate certain mechanical parameters values (in the case of neo - Hookean, Mooney - Rivlin, the parameter for Ogden, 64 and the at middle and distal thigh location for Fung) and underestimate others ( for Ogden, and for all locations and at proximal thigh l ocation for Fung). This suggested that constitutive parameter s provided by semi - 3D model depict ed experimental data in a stiffer fashion compared to the 3D model for in vivo data in current study , in terms of initial stiffness (i.e. , , , , and given in Equation ( 15 ) , with at proximal thigh being an exception ). Regarding the second approach, as shown in Figure 23 , a positive deflection differe nce value indicated that semi - 3D models show ed a more compliant behavior when compared to the 3D FE model, while a negative value indicate d a less compliant behavior. When the neo - Hookean or the Mooney - Rivli n constitutive laws were employed, the deflection differences were positive for all values of load at all locations, suggesting that semi - 3D models behave d in a more compliant way when compared to the 3D FE model for in vivo data in current study . The same behavior was observed for the Fung model at the middle and distal locations. In this case, the use of constitutive parameter s estimated for a 3D model to perform a semi - 3D analysis would lead to underestimating the load needed to cause a certain displacement and ultimately to underestimating the stress distribution within the thigh. On the other hand, when the Ogden model was employed in all locations and the Fung orthotropic model was employed at the proximal location, the deflection difference have positive values for small loads and negative values for large loads. This suggest ed that, when using parameters that have been optimized for a 3D model to perform a semi - 3D FE analysis , one would overestimate the values of loads needed to apply a larg e deflection. This would likely lead to overestimating the stress state within the tissue, specifically for the seated condition. We propose that the effect due to a choice of modeling region dimension (i.e., semi - 3D vs. 3D) can largely affect mechanical prediction. The best - fit constitutive parameter s can be 65 significantly different between a semi - 3D model and a fully 3D model when optimized using the same experimental dataset (e.g. , at middle thigh location, best - fit values of semi - 3D model was 3 .9 times that of the 3D model). Further more , when using constitutive parameter s optimized for a 3D model to perform a semi - 3D analysis , the error in the prediction can be as high as 44.56%. While we do not think that the specific conclusions based on a sin gle subject study are generalizable, due to the complexity of subject - specific soft tissue geometries and boundary conditions, we want to highlight that significant differences between semi 3D model and full 3D model exist and need to be taken into account when comparing results from different studies. Furthermore, even within the single - subject study, the difference between semi - 3D and 3D model are inhomogeneous and unpredictable. For instance, the semi - 3D model can, with respect to the 3D model, both unde restimate and overestimated the values of , based on location as shown in Table 4 , and it can both overestimate and underestimate predicted displacement, based on material model as shown in Figure 23 . Figure 23 . Mechanical behavior comparison between 3D and semi - 3D FE thigh models of participant F6 at ( a ) proximal thigh, ( b ) middle thigh, and ( c ) distal thigh. Horizontal axis represents forces applied to the soft tissue of thigh at each loading step, while the vertical axis represents deflection difference between semi - 3D FE model and 3D FE model at each loading step, i.e. , . These results suggest ed that, when describing the same in vivo mechanical behavior, constitutive parameter s obtained from one modeling condition can be significantly different than parameters obtained from a different modeling condition ; so comparison betw een them might be misleading. These results also suggest ed that one needs to be careful when using constitutive parameter s that have been estimated in a specific modeling condition to perform an analysis in a different modeling condition. After showing tha t Ogden and Fung orthotropic were the most accurate models among the ones considered here for one 66 representative participant, the best - fit constitutive parameter s were evaluated for the remaining 19 participants. A statistical analysis was carried out to h ighlight potential differences in mechanical properties between sexes and between locations. Male participants were found to have consistently higher initial tissue stiffness than female participants, regardless of constitutive model , at the proximal thigh ( and ) and middle thigh ( ) location ( Figure 24 ). This might be caused by the different compositions of cross - sections between men and women, specificall y women were found to have a higher fat content when compared to men [118] . Also, we detected a high location - dependence of consti tutive parameter s across constitutive model s. Specifically, the proximal location showed to be significantly stiffer (higher values of and ) when compared to the middle and distal location for both men and women ( Error! Reference source not found. ). F urther more , the middle location proved to be significantly stiffer when compared to the distal location for both men and women ( Error! Reference source not found. ). This suggest ed that the t high tissue increased its compliance moving distally from the buttocks. This result confirmed what have been reported in [119] , where soft tissues in the thigh - buttock area were modeled using a linear elastic isotropic material, with four differe modulus values, respectively. The material parameters estimated in that study show ed a similar stiffening trend across thigh locations as what we have reported here. The sex and location differences on initial stiffness found in the current study employing Ogden and Fung orthotropic constitutive law agree d with a previous study from this same group that used a simplified uniaxial compression model employing Mooney - Rivlin constitutive law to describe the same ex perimental dataset [106] . The sex and location influences on the constitutive parameter s shown above indicate d the importance of individual and location specific in vivo measurement of biological soft tissue in FE modeling. Further more , the present results show ed that this modeling approach has 67 the capability of detecting mechanical differences between g roups, in a consistent way, independent of the non - linear model of choice (e.g., Fung or Ogden model). We propose that, in the future, this modeling approach, supported with experimental data, will have the capability of detecting changes in mechanical properties due to other conditions, such ageing or disease. The current study has some limitations. First, the thigh geometry generation lack ed subject - specific data. While we have collected the specific surface geometry for the representative subject, the geometry of the femur and the determination of femur relative position within the thigh were based on literature so urces. Furthermore, the FE models of the remaining 19 subjects also lack ed surface geometry details, however we ensure d that the femur size and its relative position were within anatomically reasonable ranges, and the hip - joint lengths and the middle thigh circumferences of the other 19 subjects were matched to subject - specific data. Second, the thigh soft tissues were modeled by the use of a bulk material, thus, the material parameters were evaluated based on different locations rather than soft tissue typ e s . Having access to a more precise anatomy of the soft tissues within the thigh, possibly through medical imaging, will give us the capability of estimating mechanical properties correlated to specific tissues rather than locations . It has been shown by p revious MRI study that the muscle strain is more sensitive to changes in seat support surface and load distribution than strain in subcutaneous fat [120] , which could be an important characteristic to validate an anatomically accurate FE model and further increase its predicting capability. This will be the focus of future investigatio ns. Third , the bulk properties of soft tissues in the present study were modeled using isotropic material models, although the mechanical properties of soft tissues (e.g., muscle and skin) have been known to be anisotropic. Due to the lack of multi - axial i n vivo experimental data and to the choice of modeling the soft tissues as a bulk material, however, using isotropic material models is a reasonable option [115], 68 [116] . Fourth , the optimization process in this study used the simplex search algorithm ( fminsearch function in MATLAB). The use of this algorithm was motivated by the fact that it was a derivative - free method. Algorithms like trust - region would require gradient information of the objective function, however the objective function in our study had no analytic expression and performed like a black box: inputs were the material parameters and outputs were the NRMSDs between the FE simulation and experimental data. Due to the lack of analytic expressions, the convexity of the objective function cannot be confirmed, which opened the possibility of the optimized material parameters to be loca l minima instead of global minima. In the next chapter , a FE model of the buttock - thigh region with accurate participant - specific overall geometry and subdermal anatomical structure will be developed. Constitutive parameter s of different soft tissue type, i.e. , muscle, fat, skin will be evaluated. 69 A NONLINEAR HUMAN THIGH MODEL WITH HIGH ANATOMICAL AND MECHANICAL FIDELITY 4.1. Overview The goal of this chapter is to develop a FE model of the thigh that is accurate in the description of both anatomical de tails and nonlinear mechanical behaviors . A 3D thigh FE model that contained accurate geometry of muscle, fat, and skin was developed based on the MRI images of a subject. Using the optimization process described in the previous c hapter , the best - fi t constitutive parameters for each tissue component were obtained to match the compressive force - deflection behavior observed in vivo . 4.2. Method s 4.2.1. Mechanical behavior of soft tissues in vivo Thigh i ndentation test . F orce - deflection data sets from the leg of o ne subject (male, 26 years old, 185 cm, 185 lb) . The right thigh under a quadruped and a prone posture , respectively. An indenter was manually a pplied to each of the tested area on the thigh using a constant load rate. Force and deflection measurements were taken via a six - axis load cell and a linear potentiometer, respectively, calibrated to report zero deflection and force when the indenter was flush with the undeformed thigh . Load was applied until a biological barrier was reached. Three locations on the posterior side of the right thigh were tested (i.e., proximal, middle, and distal location). Skin suction test . Skin mechanical data, collected in vivo , were obtained from a previous study [121] . Briefly , to measure the mechanical properties of the s kin on the anterior surface of the forearm , Diridollou and colleagues employed a small cylinder device that could create partial 70 vacuum caus ing deformation of the skin by suction ( Figure 25 a ). An ultrasound transducer measured the resulting vertical displacement of the skin . The end of the cylinder device was attached to skin using doubled sided adhesive tape and t he central aper ture, through which the skin was under suction pressure, measured 6 mm in diameter. In this study, we extracted 20 data points from the pressure - displacement curve presented in [121] ( Figure 25 b ) and we employed this mecha nical data to inform optimization process and obtain the constitutive parameters for skin . Figure 25 . ( a ) A cylinder device that applies suction pressure on skin and measures the displacement of skin by an ultrasound transducer. . ( b ) The ascendant part of pressure - displacement data of skin under suction pressure. The solid line represents experimental data, while the dotted line represents theoretical data . ( c ) FEM g eometry to reproduce the suction test. The thickness of the skin is 0.95 mm. The diameter of the outer circle is 30 mm, corresponding to the cylinder suction device shown in ( a ) ( ). The diameter of the inner circle is 6 mm, corresponding to the aperture size shown in ( a ) ( . ( d ) Mesh of the FEM of the skin. ( a ) and ( b ) are from [121] . 4.2.2. Geometry generation Thigh model . A total of 38 MRI images that ranged from the distal end of the buttock to the proximal end of the patella of the right thigh of the subject were obtained ( Figure 26 a ) . Image 71 processing for tissue isolation was done in Adobe Photoshop CS6 . Pixels of individual tissue types were recognized by their different shading values ( Figure 26 b ) . The boundary of each tissue types was then traced , smoothed, and exported ( Figure 26 c ). The cross - secti on boundary lines of all tissues were then stacked in the 3D space in Siemens NX 12.0 to guide the generation of 3D geometry of each tissue type , i.e., fat, muscle, femur ( Figure 26 d - f ) . Skin layer was added on top of the fat layer in the meshing process. For the most part, all muscle groups were lumped together into a single section, but eventually the muscle groups sepa rate so much that it was not reasonable to continue lumping them . Toward the distal end of the thigh , it was found that the muscle started splitting from the generally homogenous whole into five separate groups. So, up to this point, the muscle groups were lumped, and after this point the muscle was split into five sections. The split in the geometry occurred at 216 mm from the proximal end of the thigh, where the separation between muscles exceeded 30 pixels (13.45 mm) in diameter. Geometry for the fat and bone were single solids for the entire length of the thigh. Skin suction model . We developed a skin model that described the experimental setup in [121] . Specifically, to represent the skin area tested , we developed a model of a cylindrical area of the skin with a height of 0.95 mm , corresponding to skin thickness , and an out er diameter of 30 mm, corresponding to the diameter of the cylinder suction device. The diameter of the center circular area was 6 mm, corresponding to the aperture size ( Figure 25 a ) . 72 Figure 26 . ( a ) A cross - section MRI image of the right thigh of the subject. ( b ) Isolated pixels of the muscle tissues from the cross - section MRI image. Femur geometry will be subtracted from the muscle geometry when the 3D models of femur and muscle are formed. ( c ) Boundary of the muscle tissues recognized from the cross - section muscle tissues. ( d ) Boundary lines of the muscle tissues stacked in the 3D s pace. ( e ) 3D geometry of the muscle tissues generated from the cross - section boundary lines. ( f ) 3D geometry of the thigh with fat, muscle, and femur. ( g ) M esh of the thigh 3D geometry. 4.2.3. Finite element model Mesh generation and boundary conditions of the th igh model . Each component of the 3D thigh geometry mentioned above (i.e., fat, muscle, femur) was meshed (HyperMesh ver.13.0) using 10 nodes tetrahedral elements with an average size of 7 mm ( Figure 26 g ) . A skin layer mesh was generate d on the thigh outer surface using 6 nodes triangular shell elements with thickness of 1.64 mm [122] . The FE model had 63,921 elements and 85, 272 nodes. Following the procedure described in c hapter 3 , we identified t hree areas at the proximal, middle, and distal locations on the p osterior side of the thigh, corresponding to the tested regions of the in vivo indentation test ( Figure 27 a - c ) . Boundary conditions were applied to the femur (i.e., constrained against displacement and rotation in any direction), to the two cut planes at proximal and distal thigh (i.e., nodes were constrained against longitudinal displacements along the femur longitudinal direction), and to the three testing areas (i.e., nodes were restricted to move only in an anterior - posterior direction ). Additionally, subject - spec ific compression forces at each location, 73 as recorded in vivo , were applied to one testing area at a time. Compression forces were applied in 20 evenly spaced steps from zero to the maximum value recorded in the experiment, and the nodal displacements of t he compressed area at each step were recorded and averaged. Simulated force - deflection curves were then compared to the experimental data to estimate the best - fit constitutive parameter s . Mesh generation and boundary conditions of the skin suction model . T he skin model was meshed (Preview 1.19) using 8 nodes hexahedral elements ( Figure 25 d ) . There were 4 layers of elements in the thickness direction. The averaged element size for the center part of the skin under pressure was 0.35 mm. The size of the elements gradually grew when it approached the edge. The FE model has 1536 elements and 2005 nodes. We applied b oundary conditions to represent surface suction pressure to th e central circular area with diameter of 6 mm, we also constrained displacements and rotations in every direction for the elements in the outer surrounding area of the skin surface , to represent the tied contact with the loading cylinder . The suction press ure was applied in 20 evenly spaced steps from zero to the maximum value of 10 k Pa. Simulated press - displacement curves were then compared to the experimental data to estimate the best - fit constitutive parameters. Figure 27 . indentation test area s at ( a ) proximal thigh, ( b ) middle thigh, and ( c ) distal thigh. ( d ) The central part of skin that is under suction pressure . 74 4.2.4. Constitutive model s and parameter estimation Constitutive models of f at and muscle . The mechanical behavior of fat and muscle were described by the first order Ogden Model. The strain energy function for fat and muscle were as follows, respectively. , ( 18 ) , ( 19 ) where the meaning of each parameter is the same as explained in s ection 3.2.4 . Constitutive model of s kin . The mechanical behavior of skin was described by the Gasser - Ogden - Holzapfel (GOH) model [50] , which was a solid mixture model that applied to incompressible solids with two preferre d directions aligned along the unit vector and . The strain energy function of skin was written as, , ( 20 ) where was a constitutive parameter with the dimension of a stress, , and were as described in Equation ( 9 ) , was a constitutive parameter with the dimension of a stress, and was a dimensionless constitutive parameter. was the structure tensor described as, ( 21 ) where was the dispersion facto rs from the normalized - periodic von Mises distribution , was the preferred fiber direction of collagen family ( Figure 28 a ) , and was the identity tensor. We consider ed two collagen fiber families in this constitutive model. The dispersion of each collagen family follow ed the normalized - periodic von Mises distribution [38] . The two collagen f amilies were assumed to share the same constitutive parameters, i.e., , , 75 . The bulk modulus - like penalty parameter was defined as to enforce i n compressibility . Figure 28 . ( a ) Structure of two crossing collagen fiber families, where and represent the preferred collagen fiber directions of collagen family 1 and 2, respectively. The bisector of the direction and direction is defined as the preferred fiber orientation of the skin, and is the angle between the direction and the preferred fiber orientation of skin. ( b ) The wrinkle line of the human leg that indicates the mean fiber direction. At the anterior side of t he thigh, wrinkle lines run ~37° with respect to the circumferential direction. At the posterior side of the thigh, wrinkle lines run in transverse direction. Around the knee area, wrinkle lines run in transverse direction. ( b ) is from [57] . Preferred fiber orientation of skin . The wrinkle lines [57] ( Figure 28 b ) were used to determine the preferred fiber orientation vector for each skin element. As described by Kraissl in his original paper, the wrinkle line pattern was transverse direction was measured as (ImageJ). For posterior aspect of the thigh and knee area, . In the FE model, the preferred fiber orientation vectors were calculated for each skin element individually. Together with the angle , the preferred direction of each collagen family could be determined for each skin element ( Figure 29 ). In the skin suction model, however, the preferred fiber orientation vectors of all elements were set in the same direction within the skin 76 surface plane since the represented skin area in the model is small (i.e., 30 mm in diameter). We assume d that the preferred fiber orientation of skin do es not vary in such a small area. Figure 29 . Preferred fiber orientation vector of each skin shell element at ( a ) anterior side of the thigh, ( b ) posterior side of the thigh, ( c ) knee area. 4.2.5. P ara meter estimation Best - fit constitutive parameters for fat and muscle . We estimated the best - fit constitutive parameters for fat and muscle using the thigh FE model without the skin layer. The FE model depend ed on four constitutive parameters (i.e., , , , and ). They were determined in the optimization process, employing the simplex search method implemented in MATLAB (i.e., function fminsea rch) . The objective function to be minimized was described by Equation ( 16 ) , ( 17 ) . The optimization iterations would stop if one of the two following stopping criteria was met: (1) the solver attempted to take a step that was smaller than , which was set to be , and the change in the val ue of the objective function during a step was smaller than , which was set to be ; (2) the objective function value was smaller than the threshold, which was set to be . The initial guess values for the parameters to start the parameter estimation process were set to be , 77 , , [28] . The best - fit constitutive parameters for fat and muscle were determined using the in vivo force - deflection data collected in quadruped and prone pos ture, respectively. Skin constitutive parameters . The skin constitutive parameters in vivo were evaluated using the skin suction FE model and the in vivo data shown in Figure 25 b . The skin suction model had five constitutive parameters in total. Two of them (i.e., , ) were structural parameters informed by literature [38] . Three of them (i.e., , , ) were determined in the optimization process, employing the simplex search method mentioned above. The FE simulation , in each optimization iterat ion , was performed in FEBio. The objective function to be minimized was described by ( 22 ) where was the number of loading steps (i.e., ), and were the displacement values of the central location of skin under pressure from the FE simulation and experimental measurement , respectively, and represented the maximum experimental values of displacement (i.e., the displacement measure ment for the maximum pressure ). The optimization iterations would stop if one of the two following stopping criteria was met: (1) the solver attempted to take a step that was smaller than , which was set to be , and the change in the value of the objective function during a step was smaller than , which was set to be ; (2) the objective function value was smaller than the threshold, which was set to be . The initial guess values to start the par ameter estimation process were set to be , , . 78 Effect of in vitro and in vivo skin constitutive parameters on thigh mechanical behavior . The effect of in vitro vs. in vivo skin constitutive parameters on the overall thigh FE model mechanical behavior during the indentation test was investigated using the thigh model with the skin layer. The skin structure parameters (i.e., , , ) were informed by literature [38], [57] , and the paramete rs for fat and muscles ( , , , and ) were t he best - fit constitutive parameters obtained from the parameter estimation process matching the in vivo indentation test data as mentioned in the beginning of this section. The n the mechanical behavior of the thigh model with in vitro skin constitutive par ameters (i.e., , , ) , which were extracted from [38] , was compared to that of the thigh model with in vivo skin constitutive parameters, which was obtained from the skin suction model parameter estimation process informed by in vivo data mentioned above. The difference was evaluated using , , ( 23 ) where , and were normalized root - mean - square deviation (NRMSD) estimated at each location as ( 24 ) where was the number of lo ading steps (i.e., ), and were the deflection values for the location from the thigh model with in vivo skin constitutive parameters and the thigh model with in vitro skin constitutive parameters, respectiv ely, and represented the maximum experimental values of deflection for each location (i.e., the deflection measure ment for the maximum load). 79 4.3. Results 4.3.1. Best - fit constitutive parameters for fat and muscle The best - fit constitutive parameters for fat and muscle tissues matching the in vivo indentation test data in quadruped and prone postures were obtained using the thigh FE model without the skin layer. The best - fit constitutive parameters were listed in Table 7 . The in vivo indentation test data were represented by red markers in Figure 30 a - b (i.e., dotted markers for proximal thigh, circular markers for middle thigh, and cross markers for distal thigh). The force - deflection behavior of FE simulation with best - fit constitutive parameters were represented by lines (i.e., solid line for prox imal thigh, dashed line for middle thigh, and dahs - dotted line for distal thigh). The normalized RMSD between FE simulation and in vivo indentation test data are right below 15% for both quadruped posture and prone posture, as shown in Table 7 . Figure 30 . ( a ) Force - deflection behavior of thigh soft tissue from in vivo indentation data in quadruped posture and the FE simulation. ( b ) Force - deflection behavior of thigh soft tissue from in vivo indentation data in prone posture and the FE simulation. 80 Table 7 . Best - fit constitutive parameters for fat and muscle tissues that describe the in vivo indentation test data collected in quadruped and prone postures. Testing posture of in vivo data Total NRMSD Quadruped 1.80 k Pa 5.75 0.36 k Pa 8.43 15.0% Prone 0.26 k Pa 9.91 3.13 k Pa 53.97 14.9% 4.3.2. Skin constitutive parameters in vivo The best - fit constitutive parameters for the in vivo skin suction test, obtained through the parameter estimation process , were , , . Figure 31 a showed the displacement contour map of the tested skin area under the maximum suction pressure of . Figure 31 b show ed the pressure - displacement curve of the central point of tested area during the suction pressure load, for both in vivo data (represented by circular markers) and FE simulation data (represented by solid line). The normalized RMSD be tween in vivo data and FE simulation was 1.15%. Figure 31 . ( a ) Displacement contour map of the skin suction model with best - fit constitutive parameters with maximum pressure load. The displacement unit in the figure is millimeter. ( b ) Pressure - displacement curves of experimental data collected in vivo and FE simulation data with best - fit constitutive parameters. 4.3.3. Effect of in vi tr o and in vivo skin constitutive parameters on thigh m echanical behavior The force - deflection behavior of thigh FE model with in vitro and in vivo skin constitutive parameters were compared to that of the thigh in vivo indentation test data. The best - fit constitutive 81 parameters for fat and muscle obtained for quadruped posture in s ection 4.3.1 were applied to the thigh models. In Figure 32 a - b , t he in vivo indentation test data were represented by red markers, specifically dotted markers for proximal thigh, circular markers for middle thigh, and cross markers for distal thigh . In addition, t he force - deflection behavior of FE simulation were represented by lines, specifically solid line for proximal thigh, dashed line for middle thigh, and dahs - dotted line for distal thigh . Figure 31 a showed model results for the FE model of the t high employing in vitro skin parameters , and Figure 31 b show ed model results for the FE model of the thigh employing in vivo skin parameters ( T able 8 ). The normalized RMSD between FE simulation and in vivo indentation test data collected in the quadruped posture are shown in T able 8 . Figure 32 . Force - deflection behavior of thigh soft tissue collected from in vivo indentation test and FE simulation. The in vivo force - deflection curve in ( a ) and ( b ) are the same data from the quadruped posture . The FE simulated force - deflection curves in ( a ) are from thigh FE model with in vitro skin constitutive parameters. The FE simulated force - deflection curves in ( b ) are from thigh FE model with in vivo skin constitutive para meters. 82 T able 8 . The constitutive parameters of the thigh FE model with in vitro skin constitutive model and the thigh FE model with in vivo skin constitutive model. The NRMSD values show the normalized root - mean - square deviations between the FE simulated force - deflection curve and the in vivo indentation data collected in quadruped posture. Thigh FE model with in vitro skin constitutive model Thigh FE model with in vivo skin constitutive model NRMSD 157.75% 25.69 % 100 . 7 k Pa 7 .6 k Pa 24 53 0 k Pa 135 . 5 k Pa 0.1327 18.4 41° 0.15 1.80 k Pa 5.57 0.36 k Pa 8.43 4.4. Discussion and conclusions Previous studies that model ed thigh - buttock area soft tissue mechanical behaviors either model ed the bulk properties of all soft tissues ( i.e., skin, fat, muscle) together [21], [27] , or model ed skin and fat as one component [19], [28] . The only few studies that modeled the skin separately used simple phenomenological constitutive models [16], [20] . In this study, we constructed a thigh model with an added skin layer and included the collagen fiber structure in the skin layer to account for the nonlinear and anisotropic behavior of skin . The use of a location - specific continuous fiber distribution to model the skin layer, resulted in a significantly high computational cost. The FE simulation of each set of constitutive parameters took more than 3 hours, and usually hundreds of sets of constitutive parameters need ed to be evaluated to obtain the best - fit results. To keep the computational time to a reasonable length was important in the framework of this study. I n fact , an excessively high computational cost would prevent this model to be used in a clinical setting. With this goal in mind, first we performed the parameter estimation 83 for fat and muscle constitutive parameters only, employing the subject specific FE model with no skin layer, and then we quantified the effect of introducing a skin layer with material parameters estimated from in vivo and in vitro data collected separatel y. Although our FE model was intended to simulate the thigh - buttock area while seated, indentation tests could be difficult while in the seated posture, especially for people with disabilities. Thus, a quadruped and a prone posture were adopted for indentation tests. The quadruped posture had the same hip and knee flexion as the seated posture , while the prone pos ture closely resembles the supine position in which MRI images were taken. The best - fit constitutive parameters obtained from different in vivo posture differ ed significantly (e.g., and obtained from prone posture data were more than 8 an d 6 times of those obtained from quadruped posture data, respectively). We propose that the different reference configurations of soft tissues in the FE model and in vivo indentation tests could be a contributing factor. The thigh FE model was developed ba sed on MRI images taken in supine posture, while the indentation test was taken in quadruped and prone posture. The hip and knee flexion would change the reference configuration of soft tissues for different postures, which was not represented correspondin gly in the FE model. This could potentially be solved by building the FE model using MRI images taken at the same posture in which in vivo data were collected. However, the lack of data of a fully relaxed configuration for all soft tissues could still lead to discrepancy in best - fit constitutive parameters that match in vivo data collected in different postures. Our best - fit constitutive parameters of fat and muscle for in vivo data of different postures strongly suggest ed the importance of the proper refer ence configuration for FE model development and in vivo data collection. It is common for simulation studies to extract constitutive parameters from literature to apply to the FE model directly [20], [21] . However, soft tissues tested in different environment 84 (i.e., in vivo vs. in vitro) exhibit very d ifferent mechanical properties [38], [121] . The GOH model adopted in this c hapter for skin mechanical properties was also used in [38] , which obtained the best - fit GOH model parameters for human back skin samples tested in vitro . However, when these constitutive parameters were applied to the skin layer of the thigh m odel presented here , along with the best - fit in vivo constitutive parameters for fat and muscle, the FE model output differ ed from in vivo indentation test drastically ( Figure 32 a ), with normalized RMSD being 157.75%. For this reason , we developed a skin FE model ( Figure 27 d ) to search for constitutive parameters of skin that match in vivo experimental data. Pressure suction, as an important means for in vivo skin mechanical properties testing [121], [123], [124] , was chosen as the source of in vivo data. The best - fit constitutive parameters for skin estimated from in vivo data, differ ed from those estimated from in vitro data of several orders of magnitude, as shown in T able 8 . When the constitutive parameters for skin estimated from the skin suction model were used to the thigh FE model, along with same in vivo constitutive parameters for fat and muscle used previously, the difference between FE model output and in vivo data from indentation test dropped to a relative low level ( Figure 32 b ), with normalized RMSD being 25.69%. Our results here show that (1) the skin layer has significant contribution to the thigh mechanical properties , and (2) i n vitro skin constitutive parameters might not be suitable to model soft tissue mechanical behavior in vivo . As to the limitations of the current study, the interface interaction between different soft tissues were not counted in current FE model . The deep fascia that separated muscles into differ ent fascia compartments and allow ed relative movement between muscle - bone interface [125] were not described in current model. Another known relative movement that could happen at the soft tissue interface is between skin layer and the underlying fascia, which h as been shown by ultrasound transducer during skin suction test [121] . However, in vivo data regarding the interface 85 movement (i.e., friction coefficients, relative movement limits) are still rare in literature. The expensive computational cost of friction contact in the FE model could also make solutions fo r parameter estimation process hard to achieve. In additio n, it should be noted that best - fit constitutive parameters obtained for both the skin suction model and the thigh indentation model were not the sole solutions to describe the in vivo experimental data. The first factor affecting the solution values was the stopping criteria of the para meter estimation process (i.e., , , and in s ection 4.2.5 ). We chose the stopping criteria in consideration of both accuracy (i.e., difference between simulation output and in vivo data) and computational cost (i.e., calculating time). Different stopping criteria c ould lead to different sets of constitutive parameters. Even the stopping criteria was the same, different initial guesses or fitting approach might lead to different best - fit solutions. For instance, in a different parameter estimation approach for the sk in suction model, the value was fixed at 0.1327 [38] , with initial guess , and the same stopping criteria as mentioned in s ection 4.2.5 . The resulting best - fit values were , . This meant the following two sets of parameters: , , and , , are deemed to represent the same pressure - displacement behavior in the current skin suction model setup. More data that describe ground matrix properties (i.e., ) and collagen properties (i.e., , ) separately could help identify the unique solution. With the current available in vivo data for skin suction model and thigh indentation model, no alternative solutions, if found, can be rejected. 86 CONC LUSIONS AND FUTURE WORK This research was aimed at proposing a finite element modeling approach that can accurately describe the nonlinear mechanical behavior of soft tissues in thigh - buttock area, so that the stress/strain distribution within soft tissue can be described precisely for a be tter understanding of pressure ulcers formation mechanism. First, we conducted a comprehensive investigation of skin biomechanics using rat skin samples . D ifferences in mechanical properties between the upper and lower back , and the axial and transverse direction were detected. The collagen network structure was quantified through automated image analysis routines and manual measurements from histology images . Results show ed that: (1) Collagen fiber bundles on the back o f rats distribute d about two predominate oblique - crossed directions: to spine. (2) Collagen content of skin at the lower back was higher than the upper back , with vs. . A micro structurally based c onstitutive model that descr ibed the nonlinear, anisotropic mechanical behavior of skin was proposed and showed good agreement with experimental data. This study contribute d valuable primary data to the skin mechanics research community and contribu ted to improved understanding of va riation and similarities in skin mechanical properties across different species. Second, we looked into finite element modeling parameter estimation regarding non - linear mechanical behavior of human thigh soft tissue s . We proposed an optimization process that addressed the non - linear characteristic of in vivo force - deflection data, which has been neglected in the literature. For the first time, we compared the ability of four widely used constitutive model s in describing soft tissue non - linear behavior. Re sults show ed that first order Ogden and Fung orthotropic model were more suitable to describe the nonlinear mechanical behavior of soft tissues, compared to neo - Hookean and Mooney - Rivlin model. Also, for the first time in literature, we investigated the ef fect of the choice of modeling region dimension (3D vs. semi - 3D). Results 87 show ed that the mechanical properties between 3D and semi - 3D model differ ed significantly . Therefore, semi - 3D FE models were not recommend for soft tissues modeling. We also reported constitutive parameter range for Ogden and Fung orthotropic models based on in vivo force - deflection data for 20 p articipants. This study provided deep insights and reliable data for researchers interested in non - linear mechanical behavio r of human soft tissues. In the last c hapter, we focus ed on the in vivo mechanical properties of different soft tissue groups. We d eveloped a human thigh model with detailed geometric details using MRI images. A state - of - the - art FE model of skin layer was created and the best - fit constitutive parameters were obtained using in vivo data in literature . Fat and muscle mechanical behavior were described using first order Ogden model and b est - fit constitutive parameters were obtained using subject - specific in vi vo experimental data collected in different postures (i.e. quadruped and prone posture). Results show ed that testing postures in which in vivo data were collected could affect the constitutive parameters of soft tissues. By comparing the force - deflection b ehavior of the thigh model with in vitro and in vivo skin constitutive parameters, we conclude d that (1) the skin layer had significant contribution to the thigh mechanical properties, and (2) in vitro skin constitutive parameters might not be suitable to model soft tissue mechanical behavior in vivo . Future work will focus on (1) modification/simplification of skin FE model to achieve shorter solvin g time, (2) stress/strain field within soft tissues introduced by sitting action, with consideration of the sitting interface condition (i.e., friction, sitting surface contours), and (3) the interaction/relative movement between different soft tissue groups. 88 BIBLIOGRAPHY 89 BIBLIOGRAPHY [1] Ostomy Wound Manage , vol. 54, no. 10, pp. 26 8, 30 5, Oct. 2008. [2] M. Reddy, S. S. Gill, and P. A. Ro JAMA , vol. 296, no. 8, pp. 974 984, Aug. 2006. [3] age - sex specific all - cause and cause - specific mortality for 240 causes of death, 1990 - 2013: Lancet , vo l. 385, no. 9963, pp. 117 171, Jan. 2015. [4] Medical News Today , Dec - 2017. [Online]. Available: https://www.medicalnewstoday.com/articles/173972.php. [Accessed: 10 - Apr - 2018]. [5] T. J. Ry Pressure Sores - Clinical Practice and Scientific Approach , Palgrave, London, 1990, pp. 141 152. [6] decubitus ulce Journal of Biomechanics , vol. 14, no. 12, pp. 879 881, Jan. 1981. [7] surface pressure - Journal of Rehabilitation Research and Development; Washington , vol. 36, no. 2, pp. 109 20, Apr. 1999. [8] Arch Phys Med Rehabil , vol. 62, no. 10 , pp. 492 498, Oct. 1981. [9] Archives of Physical Medicine and Rehabilitation , vol. 84, no. 4, pp. 616 619, Apr. 2003. [10] A. A. M transcutaneous oxygen level characterizations in human skin with changes in normal and shear loads Clinical Biomechanics , vol. 25, no. 8, pp. 823 828, Oct. 2010. [11] Clinical Biomechanics , vol. 28, no. 5, pp. 574 578, Jun. 2013. 90 [12] Annals of Internal Medicine , vol. 94, no. 5, pp. 661 666, 1981. [13] D. M. Brienza, P. E. Karg, M. Jo Geyer, S. Kelsey, and E. Tr pressure ulcer incidence and buttock - seat cushion interface pressure in at - risk elderly Archives of Physical Medicine and Rehabilitation , vol. 82, no. 4, pp. 529 533, Apr. 2001. [14] e - Br J Nurs , vol. 6, no. 6, pp. 323 328, Mar. 1997. [15] E. Linder - - induced pressure sores in rat Jour nal of Applied Physiology , vol. 96, no. 6, pp. 2034 2049, Jun. 2004. [16] C. W. J. Oomens, O. F. J. T. Bressers, E. M. H. Bosboom, C. V. C. Bouten, and D. L. Bader, Computer Methods in Biomechanics and Biomedical Engineering , vol. 6, no. 3, pp. 171 180, Jun. 2003. [17] muscle tissue with reference to pressure ul [18] - Element Biomechanical Model for Evaluating Buttock Bioengineering Research of Chronic Wounds , pp. 181 205, 2009. [19] E. Linder - Ganz conditions in sub - dermal tissues during sitting: A combined experimental - MRI and finite Journal of Biomechanics , vol. 40, no. 7, pp. 1443 1454, 2007. [20] M. Makhsous, D. IEEE Transactions on Neural Systems and Rehabilitation Engineering , vol. 15, no. 4, pp. 517 525, D ec. 2007. [21] M. M. Verver, J. van Hoof, C. W. J. Oomens, J. S. H. M. Wismans, and F. P. T. Baaijens, Computer Methods in Biomechanics and Biomedical Engineering , vol. 7, no. 4, pp. 193 203, Aug. 2004. [22] D. Brienza et al. Journal of the American Geriatrics Society , vol. 58, no. 12, pp. 2308 2314, Dec. 2010. [23] G. N. He time and cardio - metabolic biomarkers in US adults: NHANES 2003 Eur Heart J , vol. 32, no. 5, pp. 590 597, Mar. 2011. 91 [24] B. G. Raijmakers, M. G. Nieuwenhuizen, H. Beckerman, American Journal of Physical Medicine & Rehabilitation , vol. 94, no. 2, pp. 101 113, Feb. 2015. [25] A. V. Rowlands et al. entary behavior with the GENEActiv: introducing the [26] BMJ , vol. 332, no. 7539, pp. 472 475, Feb. 2006. [27] R. Ragan, T. W. Kernozek, M. Bidar, and J. W. Matheson - interface pressures on Archives of Physical Medicine and Rehabilitation , vol. 83, no. 6, pp. 872 875, Jun. 2002. [28] R. M. A. Al - Dirini, M. P. Reed, J. Hu, and D. Thewli Ann Biomed Eng , vol. 44, no. 9, pp. 2805 2816, Sep. 2016. [29] The handbook of material s behavior models , vol. 3, pp. 1049 1063, 2001. [30] bed - The Journal of Pathology and Bacteriology , vol. 66, no. 2, pp. 347 358, Oct. 1953. [31] I. British Journal of Dermatology , vol. 89, no. 4, pp. 383 393, Oct. 1973. [32] - dimensional elastic properties of human Medical Engineering & Physics , vol. 17, no. 4, pp. 304 313, Jun. 1995. [33] F. Xu and T. Lu, Introduction to Skin Biothermomechanics and Thermal Pain . Berlin, Heidelberg: Springer Berlin Heidelberg, 2011. [34] - Biomedical engineering , vol. 4, no. 7, pp. 322 327, 1 969. [35] P. P. van Zuijlen et al. superior to multi - J. Pathol. , vol. 198, no. 3, pp. 284 291, Nov. 2002. [36] H. J. C. de Vries, D. N. H. Enomoto, J. van Marle, P. P . M. van Zuijlen, J. R. Mekkes, and Laboratory Investigation , vol. 80, no. 8, pp. 1281 1289 , Aug. 2000. 92 [37] Skin Research and Technology , vol. 17, no. 2, pp. 149 159, May 2011. [38] A. Ní Annaidh et al. Ann Biomed Eng , vol. 40, no. 8, pp. 1666 1678, Mar. 2012. [39] atural Tension in the Mechanical Testing J Investig Dermatol , vol. 69, no. 3, pp. 310 314, Sep. 1977. [40] model for skin: Experimental measurements , finite element modelling and identification of Journal of the Mechanical Behavior of Biomedical Materials , vol. 18, pp. 167 180, Feb. 2013. [41] tion of the rat Bioengineered , vol. 6, no. 3, pp. 153 160, May 2015. [42] - strain relationships in flat J ournal of Biomechanics , vol. 12, no. 6, pp. 423 436, Jan. 1979. [43] Journal of biomechanics , vol. 16, no. 1, pp. 1 12, 1983. [44] P. Tong and Y. - - strain relationship Journal of Biomechanics , vol. 9, no. 10, pp. 649 657, 1976. [45] Finite Journal of Biomechanics , vol. 3, no. 1, pp. 111 124, Jan. 1970. [46] M. Comninou a - strain nonlinearity of connective Journal of Biomechanics , vol. 9, no. 7, pp. 427 433, Jan. 1976. [47] J. W. Y. Jor, M. D. Parker, A. J. Taberner, M. P. Nash, and P. M. F. and experimental characterization of skin mechanics: identifying current challenges and WIREs Syst Biol Med , vol. 5, no. 5, pp. 539 556, Sep. 2013. [48] implementation of Computer Methods in Applied Mechanics and Engineering , vol. 135, no. 1 2, pp. 107 128, Aug. 1996. [49] erial parameters Biomech Model Mechanobiol , vol. 10, no. 5, pp. 767 778, Oct. 2011. 93 [50] with dist Journal of The Royal Society Interface , vol. 3, no. 6, pp. 15 35, Feb. 2006. [51] Plastic and Reconstruct ive Surgery , vol. 27, no. 6, p. 597, Jun. 1961. [52] Plastic and Reconstructive Surgery , vol. 104, no. 1, p. 208, Jul. 1999. [53] Clinical Anatomy , vol. 27, no. 2, pp. 162 168, 2014. [54] G. Dupuytren, Traité théorique et pratique des blessures par armes de guerre . H. Dumont, 1835. [55] British Journal of P lastic Surgery , vol. 31, no. 1, pp. 3 8, Jan. 1978. [56] Operative Techniques in General Surgery , vol. 4, no. 3, pp. 199 206, Sep. 2002. [57] LINES FOR ELECTIVE Plastic and Reconstructive Surgery , vol. 8, no. 1, p. 1, Jul. 1951. [58] Plast Reconstr Surg , vol. 73, no. 1, pp. 144 150, Jan. 1984. [59] A - plasties on scars, and British Journal of Plastic Surgery , vol. 15, pp. 242 254, Jan. 1962. [60] and Stress Analysis of Proceedings of the IMechE , vol. 208, no. 1, pp. 9 17, Mar. 1994. [61] an integrative experimental Clinical Biomechanics , vol. 15, no. 3, pp. 217 219, 2000. [62] E. Linder - Ganz, N. Shabshin, Y. Itzchak, Z. Yizhar, I. Siev - stresses in sub - dermal tissues of the buttocks are greater in paraplegics than in healthy Journal of Biomec hanics , vol. 41, no. 3, pp. 567 580, 2008. [63] new sitting concept designed for prevention of pressure ulcer on the buttock using finite Med Bio Eng Comput , vol. 45, no. 11, pp. 1079 1084, Nov. 2007. 94 [64] J Biomech Eng , vol. 100, no. 2, pp. 79 87, May 1978. [65] E. M. H. Bosboom, M. K. C. Hesselink, C. W. J. Oomens, C. V. C. Bouten, M. R. Drost, Journal of Biomechanics , vol. 34, no. 10, pp. 1365 1368, Oct. 2001. [66] te elements model of the human body: geometry and non - Proceedings of the TMCE , 2002, vol. 2002. [67] K. K. Ceelen et al. With MR Tagging: A Contribution to Pre J Biomech Eng , vol. 130, no. 6, pp. 061015 - 061015 8, Oct. 2008. [68] E. Linder - Ann Biomed Eng , vol. 35, no. 12, pp . 2095 2107, Dec. 2007. [69] R. M. A. Al - Clinical Biomechanics , vol. 30, no. 7, pp. 662 668, Aug. 2015. [70] S. E. Sonenblum, S. H. Sprigle, J. M. Cathcart, and R Journal of Tissue Viability , vol. 24, no. 2, pp. 51 61, May 2015. [71] - dimensional computer model of the human buttocks, Journal of Rehabilitat ion Research and Development; Washington , vol. 31, no. 2, pp. 111 9, 1994. [72] - elastic Biomech Model Mechanobiol , vol. 4, no. 2 3, pp . 190 199, Nov. 2005. [73] - Layer Structural Model for Arterial Walls With a Fung - J Biomech Eng , vol. 126, no. 2, pp. 264 275, May 2004. [7 4] - dimensional mechanical properties of rabbit skin II. Journal of biomechanics , vol. 7, no. 2, pp. 171IN9175 174182, 1974. [75] Journal of Biomedical & Pharmaceutical Engineering , vol. 2, no. 1, pp. 22 28, 2008. [76] axial and circumferential mechanical properties of the skin tissue using exp erimental testing Computer Methods in Biomechanics and Biomedical Engineering , vol. 18, no. 16, pp. 1768 1774, Dec. 2015. 95 [77] IEEE Engineering in Medicine and Biology Magazine , vol. 10, no. 2, pp. 55 57, Jun. 1991. [78] E. Balli et al. Ecotoxicology and environmental safety , vol. 72, no. 3, pp. 889 894, 2009. [79] Journal of the Mechanical Behavior of Biomedical Materia ls , vol. 41, pp. 241 250, Jan. 2015. [80] Arch Dermatol Res , vol. 277, no. 6, pp. 484 488, 1985. [81] A. Delalleau, G. Josse, J. - M. Lagarde, H. Zahoua ni, and J. - Skin research and Technology , vol. 14, no. 2, pp. 152 164, 2008. [82] human skin in Skin Res Technol , vol. 20, no. 2, pp. 127 135, May 2014. [83] of the anisotropic mechanical prop Journal of the Mechanical Behavior of Biomedical Materials , vol. 5, no. 1, pp. 139 148, Jan. 2012. [84] Microstructurally Motivated Model of Ar terial Wall Mechanics with Mechanobiological Ann Biomed Eng , vol. 42, no. 3, pp. 488 502, Mar. 2014. [85] Arteries in Health and Disease: Advantages of and M Ann Biomed Eng , vol. 41, no. 7, pp. 1311 1330, Jul. 2013. [86] C irculation , vol. 71, no. 4, pp. 740 744, Apr. 1985. [87] Braz J Morphol Sci , vol. 22, no. 2, pp. 97 104, 2005. [88] P . Arun Gopinathan, G. Kokila, M. Jyothi, C. Ananjan, L. Pradeep, and S. Humaira Nazir, Scientifica , 2015. [Online]. Available: https://www.hindawi.com/journals/scientifica/2015/802980/abs/. [Accessed: 27 - Sep - 2017]. 96 [89] R. Rezakhaniha et al. the arterial adventitia using confocal laser scanning microsc Biomech Model Mechanobiol , vol. 11, no. 3 4, pp. 461 473, Mar. 2012. [90] S. Zeinali - Davarani, Y. Wang, M. - Collagen Fiber Undulation to Regional Biomechanical Properties Along Porcine Thoracic J Biomech Eng , vol. 137, no. 5, p. 051001, May 2015. [91] S. Baek, R. L. Gle Potential utility in computations of fluid Computer Methods in Applied Mechanics and Engineering , vol. 196, no. 31, pp. 3070 3078, Jun. 2007. [92] F. Aur - time Acta Derm Venereol , vol. 73, no. 5, pp. 344 347, Oct. 1993. [93] M. J. Muñoz et al. amage and Journal of Biomechanics , vol. 41, no. 1, pp. 93 99, 2008. [94] DRM , vol. 208, no. 2, pp. 112 119, 2004. [95] H. Takeuchi et al. Experimental Animals , vol. 60, no. 4, pp. 373 384, 2011. [96] J. C. J. Wei, G. A. Edwards, D. J. Martin, H. Huang, M. L. Crichton, and M. A. F. Kendall, - medical Scientific Reports , vol. 7, no. 1, p. 15885, Nov. 2017. [97] Biomed Eng Online , vol. 12, p. 125, Dec. 2013. [98] C. Escoffier, J. de Rigal, A. Rochefort, R. Vasselet, J. - L. Lévê - Journal of Investigative Dermatology , vol. 93, no. 3, 1989. [99] High and J Biomech Eng , vol. 111, no. 2, pp. 136 140, May 1989. [100] Journal of Experimental Biology , vol. 29, no. 3, pp. 454 468, Se p. 1952. [101] Ann Surg , vol. 120, no. 4, pp. 431 448, Oct. 1944. 97 [102] Anat. Rec. , vol. 89, no. 1, pp. 75 85, May 1944. [103] R. D. Ma rangoni et al. Annals of the New York Academy of Sciences , vol. 136, no. 16, pp. 441 453, Nov. 1966. [104] anical properties Am J Vet Res , vol. 53, no. 5, pp. 788 792, May 1992. [105] International Journal for Nu merical Methods in Biomedical Engineering , vol. 27, no. 11, pp. 1793 1811, Nov. 2011. [106] the in vivo material properties of the seated human buttocks and thig International Journal of Non - Linear Mechanics , Oct. 2018. [107] GRABCAD , 22 - Aug - 2016. [Online]. Available: https://grabcad.com/library/james - human - body - for - scale - 1. [Accessed: 10 - Apr - 2018]. [108] S. Strandberg, M. - tomography measurements in assessment of thigh muscle cross - sectional area and BMC Medical Imaging , vol. 10, p. 18, Aug. 2010. [109] C. C. Gordon et al. 1989. [110] Gender J Bone Miner Res , vol. 16, no. 7, pp. 1291 1299, Jul. 2001. [111] Phil. Trans. R. Soc. Lond. A , vol. 240, n o. 822, pp. 459 490, Jan. 1948. [112] Phil. Trans. R. Soc. Lond. A , vol. 242, no. 845, pp. 173 195, Dec. 1949. [113] R. W. Og on the correlation of theory and Proc. R. Soc. Lond. A , vol. 326, no. 1567, pp. 565 584, Feb. 1972. [114] - invariant fo Journal of Biomechanics , vol. 42, no. 6, pp. 781 785, Apr. 2009. 98 [115] B. Manafi - Khanian, L. Arendt - Nielsen, and T. Graven - - based leg model used to simulate biomechanical phenomena during cuff algometry: a fin Med Biol Eng Comput , vol. 54, no. 2, pp. 315 324, Mar. 2016. [116] Biomechanical Properties of the Pubovisceralis Muscle Using a Genetic Algorithm and the J Biomech Eng , vol. 141, no. 1, pp. 011009 - 011009 11, Oct. 2018. [117] K. M. Myers et al. Journal of Biomechanics , vol. 48, no. 9, pp. 1533 1540, Jun. 2015. [118] - sectional area and Europ. J. Appl. Physiol. , vol. 68, no. 2, pp. 148 154, Feb. 1994. [119] C. Mergl, T. Anton, R. Madrid - SAE Technical Paper 2004 - 01 2132, Jun. 2004. [120] R. M. A. Al - Dirini, J. Nisyrios, M. P. Reed, and D. Th - static response to loading of sub - dermal tissues in the human buttock using magnetic Clinical Biomechanics , vol. 50, pp. 70 77, Dec. 2017. [121] S. Diridollou et al. Skin Research and Technology , vol. 6, no. 4, pp. 214 221, 2000. [122] A. Laurent et al. ultrasound to as Vaccine , vol. 25, no. 34, pp. 6423 6430, Aug. 2007. [123] Bed Sore Biomech anics , R. M. Kenedi and J. M. Cowden, Eds. Macmillan Education UK, 1976, pp. 109 117. [124] F. M. Hendriks, D. Brokken, J. T. W. M. V. Eemeren, C. W. J. Oomens, F. P. T. Baaijens, - experimental method to characterize th e non - linear Skin Research and Technology , vol. 9, no. 3, pp. 274 283, 2003. [125] J Biomech Eng , v ol. 132, no. 8, pp. 084502 - 084502 4, Jun. 2010.