19 Computational modeling, such as finite element analysis, is employed in a range of biomechanics 20 specialties, including impact biomechanics and surgical planning. These models rely on accur ate material 21 properties for skeletal muscle , which comprises roughly 40% of the human body . Due to surrounding 22 tissues, compressed skeletal muscle in vivo likely experiences a semi -confined state . Nearly all previous 23 studies investigating passively compres sed muscle at the tissue level have focused on muscl e in 24 unconfined compression . The goals of this study were to (1) examine the stiffness and time -dependent 25 material properties of skeletal muscle subjected to both confi ned and unconfined compression (2) d evelop 26 a model that captures passive muscle mechanics under both conditions and (3) determine the extent to 27 which different assumptions of volumetric behavior affect model results. Muscle in confined compression 28 exhibite d stiffer behavior, agreeing with pr evious assumptions of near -incompressibility . Stress relaxation 29 was found to be faster under unconfined compression , suggesting there may be different mechanisms that 30 support load these two conditions . Finite element cal ibration was achieved through nonlin ear optimization 31 (normalized root mean square error <6%) and model validation was strong (normalized root mean square 32 error <17%) . Comparisons to commonly employed assumptions of bulk behavior showed that a simple 33 one pa rameter approach does not accurately simulate confined compression . We thus recommend the use 34 of a properly calibrated, nonlinear bulk constitutive model for modeling of skeletal muscle in vivo . Future 35 work to determine mechanisms of passive muscle stiffne ss would enhance the efforts present ed here. 36 1 Introduction 37 Skeletal muscle comprises approximately 40% of the mass of the human body [1]. Computational modeling 38 of passive skeletal muscle is thus essential to simulations of impact biomechanics [2]–[8], rehabilitation 39 engineering [9], [10] , surgical planning [11], [12] , and bed sore development [9], [1 3]. These models rely 40 on accurate material properties for skeletal muscle , which have been shown to be anisotropic [14], [15] , 41 time dependent [3]–[5], [16], [17] , non -linear [3], [17] , and asymmetric in regards to tension and 42 compression [18], [19] . However, the compressive behavior of skeletal muscle is not fully understood, 43 particularly regarding the differences in muscle response to in vivo loading conditions [3]–[5], [20] . 44 It is likely that in vivo muscle experiences a variation between confined and unconfined volumetric 45 boundary conditions [3]–[5], [21] , where semi -confinement is created by tissues surrounding the muscle. 46 Nearly all previous studies investigating passively compressed muscle at the tissue level have focused on 47 muscle in unconfined compression (UC) [3]–[5], [22] –[24], where the sample is free to expand laterally 48 when loaded . One group has investigated muscle under anisotrop ic semi -confined compression, and the 49 specific confinement was found to affect both muscle structural deformation and mechanical response [14], 50 [15]. Traditionally, muscle has been model led as a n early incompressible hyperelastic material [14], [24] – 51 [26]. However, to the best of the authors’ knowledge, there have been no investigations of skeletal muscle 52 in fully confined compression (CC), where volumetric strain is applied and the assumption of near 53 incompress ibility ca n be directly tested. This gap in understanding the effects of volumetric boundary 54 conditions (UC and CC) on the compressive properties of skeletal muscle affect the subsequent models 55 derived to predict skeletal muscle and whole body behavior. 56 The time d ependent nature of muscle can be observed in significant stress relaxation following compressive 57 deformation [3]–[5], [16], [17] . Stress relaxation tests have been thus used extensively to characterize the 58 stress -strain and stress -time behavior of the tissue, and are typically accompanied by various viscoelastic 59 modelling approaches [3], [4], [17], [27] –[29]. Inverse finite elem ent methods are effective in determining 60 material properties through parameter optimization to experimental data [27], [30] –[32]. In previous 61 continuum mechanics based modeling of skeletal muscle, the assumption of near incompressibility leads to 62 a decoupling of the volumetric (volume changing) and isochoric (shape changing) responses of the 63 hyperelastic model [23], [33] . We developed a non -linear hyper -viscoelastic f inite element model to 64 simulate both UC and CC testing conditions concurrently that was calibrated using inverse finite element 65 analysis through a nonlinear optimization protocol . This model was then used to investigate various 66 assumptions about muscle com pressibility, and what the most appropriate modelling approaches may be 67 for passivel y compressed skeletal muscle. 68 The goals of this study were to (1) examine the stiffness and time -dependent material properties of skeletal 69 muscle subjected to two boundary conditions (UC and CC) (2) develop a computational model that captures 70 the behavior of muscle subject to these dif ferent volumetric boundary conditions and (3) determine the 71 extent to which different assumptions of volumetric behavior affect model results. We hypothesize that the 72 material properties of skeletal muscle differs in confined versus unconfined compression. Since in vivo 73 muscle behavior is likely to be semi -confined , we also believe that by considering both conditions we 74 improve the accuracy of f inite element models of skeletal muscle. 75 76 2 Methods 77 2.1 Sample Preparation 78 Whole tibialis anterior (TA) muscles were isolated from seven female porcine hind limbs following 79 sacrifice [14], [15], [23], [24], [26], [34] . Partici pants did not handle live animals as all tissue was acquired 80 from a local abattoir. Connective tissues and fat were removed from TA using standard dissection 81 instruments. Two sample geometries were identified for experimentation in this study: transverse o riented 82 cuboids for UC and transverse cylindrical plugs for CC. The cuboids (height = 7.2 ± 0.9 mm, cross -sectional 83 area=149.6±25.7 mm2, n = 15 for fast compression and n = 14 for slow compression) [14] were acquired 84 using scalpels and a custom tissue slicer . Sample height was measured via micrometer and area was 85 calculated through image analysis prior to testing . The cylindrical plugs (height = 7.1 ± 0.6 mm, n = 16 for 86 fast compressio n and n = 15 for slow compression ) were obtained using metal hole punch (ø = 10mm) and 87 scalpels. Since passive skeletal muscle is an incredibly soft tissue, using a larger metal punch ensured that 88 the sample tightly fit into the compression well. Four to five cuboids and cylindrical plugs were excised 89 from each muscle and the samples were not paired between unconfined and confined compression. All 90 samples were taken from muscle midbelly and were kept hydrated by phosphate buffered saline throughout 91 testing [24], [26] , [34] . To limit effects of rigor mortis, all testing was completed within eight hours of 92 sacrifice [3]–[5], [17], [29], [35] . Tissue damage was controlled at two points in experimental protocol. 93 Firstly, duri ng dissection, a custom slicer and high -profile histology blades or surgical scalpels were used 94 to cut the sample a s few times as possible. Secondly, after each UC and CC tests, the sample was visually 95 checked for damage, and any damaged sample was discard ed. 96 97 2.2 Experimentation 98 Custom instrumentation was fabricated to perform UC and CC using a uniaxial tabletop Inst ron 3366 tensile 99 testing system. A lightweight delrin top and a fixed stainless steel bottom platen were used for UC testing 100 (Figure 1A). An Al2O3 porous plunger (diameter = 6.4 mm, length = 25.5 mm) was used along with an 101 impermeable steel well (diameter = 6.9 mm, depth = 8 mm) for CC (Figure 1B). 102 103 Figure 1: Schematics showing experimental set up for ( A) unconfined compression and (B) confined 104 compression. Associated finite element geometry and boundary conditions for (C) unconfined compression 105 and (D) co nfined compression. Unconfined finite element geometry used quarter symmetry and confined 106 geometry used axisymmetric. 107 Two stress relaxation testing conditions were employed for both UC and CC conditions: fast and slow 108 compression stress relaxation [5], [29], [36] . All tests were c ompleted under transverse compression to 109 simulate the most common uniaxial physiological loading orientation [15], [22], [23], [34], [37] . For UC, 110 all samples were strained to 40% compressive no minal strain at either 40% s-1 (fast compression, n=15) or 111 5% s-1 (slow compression, n=14) [15], [38], [39] . For CC, all samples were strained to 1 5% compressive 112 nominal strain at either 15% s-1 (fast compression, n=16) or 1.5% s-1 (slow compression, n=15). All samples 113 were subject to a 400 seconds stress -relaxation hold [4]. Data was acquired by either a 10N or 100N Instron 114 load cell (2350 series) at 100 Hz. Time (seconds), extension (mm), and load (Newtons) were recorded. All 115 model calibrat ion (determination of model parameters) was completed with fast -compression data, while 116 slow -compression data were only used for model validation. 117 118 2.3 Data Analysis and Viscoelastic Modelling 119 First Piola -Kirchhoff (PK) stress 𝑃, nominal strain 𝜀, and peak stress were determined through original 120 specimen dimensions (Equation 1, where 𝐹 is the measured load and 𝐴𝑜 is the original specimen cross - 121 sectional area ) [17], [28] . Peak modulus ( 𝑃𝑝𝑒𝑎𝑘 𝜀𝑝𝑒𝑎𝑘⁄ ) was also calculated. Three relaxation ratios (1 -5s, 122 6-105s, 106 -400s) were determined to evaluate the amount of relaxation associated with various relaxation 123 times (Equation 2, where 𝑅𝑅 is the relaxation ratio and 𝑃𝑖 and 𝑃𝑗 are the first PK stress at time points 𝑖 and 124 𝑗) [29]. These three time periods were chosen as they generally characterized “short”, “medium”, and “long” 125 term relaxation for the samples tested in this study, and are similar to time periods previously used for 126 skeletal muscle in tension [29]. 127 𝑃=𝐹 𝐴0 (1) 128 𝑅𝑅=𝑃𝑖−𝑃𝑗 𝑃𝑖 (2) 129 To more finely characterize relaxation behavior, a three term Prony series quasi -linear viscoelastic model 130 (Equations 3 -4) was fit to normalized hold phase stress from all testing grou ps [17], [28], [39] –[41]. 131 𝑃(𝜀,𝑡)=∫𝐸(𝑡−𝜏) 𝑑𝜀(𝜉) 𝑑𝜉𝑑𝜉𝑡 0 (3) 132 𝐸(𝑡)=𝐸0(1−∑ 𝐸𝑖[1−exp (−𝑡 𝑖)])3 𝑖=1 (4) 133 Here stress 𝑃 is calculated from the convolution integral (Equation 3), which includes the Prony series 134 reduced relaxatio n function 𝐸(𝑡), nominal strain 𝜀, and an integration parameter 𝜉. This includes three 135 relaxation coefficients 𝐸𝑖, three time constants 𝑖, and the instantaneous modulus 𝐸0. As the purpose of this 136 analysis was to compare relaxation behavio r only, 𝐸0=1 was fixed and all data were normalized. The 137 model also accounted for the experimental overshoot in strain applied by the Instron during fast 138 compr ession. Parameter determination was performed in two steps: a Monte Carlo simulation followed by 139 a nonlinear least -squares deterministic optimization ( lsqnonlin in MATLAB) [27], [28], [30], [42] –[45]. In 140 the Monte Carlo simulation, the six parameters ( 𝐸1−3 and 𝜏1−3) were randomly varied for 100,000 141 simulations, ensuring 0<𝐸1+𝐸2+𝐸3<1 [30]. The set of parameters mi nimizing percent error between 142 normalized model and experimental hold stress was used as initial guesses for the deterministic 143 optimization, which optimized percent difference between normalized model and experimental hold 144 stresses . This approach used the global stochastic Monte Carlo method in conjunction with the precision of 145 a local deterministic approach. All modelling was performed in MATLAB (The Mathworks, Inc.). 146 147 2.4 Finite Element Modelling 148 Two finite element models of UC a nd CC geometries were dev eloped and simulated with an implicit finite 149 element approach in Abaqus/Standard (Dassault Syst èmes). (Figure 1C-D). The UC geometry was reduced 150 to a quarter of the sample by symmetry, with 175 first -order 8 -node hexahedral elements with a cubic 151 volume of 1 mm3 (type C3D8RH). This model was compressed by coupling the top surface to a reference 152 node and displacing the node to follow expe rimental dis placement. The reaction force on the reference 153 point was divided by initial area to acquire first Piola -Kirchh off stress. The CC geometry was reduced to a 154 two-dimensional axisymmetric model of cylinder with 96 first -order 4 -node quadrilateral e lements with a 155 rectangular area of 3.75 mm2 (type CAX4RH). A convergence study was performed by doubling and 156 halving the m esh densit ies, and the model outputs were virtually identical to the outputs of the original 157 model s. Displacement was prescribed for t op surface and as no lateral expansion occurred, first Piola - 158 Kirchhoff stress was determined directly from model axial str ess. Displacements in both fast compression 159 models simulated the slight experimental ov ershoot applied by the Instron. Four boundary c onditions were 160 applied to the UC model to ensure quarter symmetry while leaving exposed faces to expand due to the 161 Poisson effect (Figure 1C). Since very little sample sliding was noticed during pilot testing or 162 experimentation, sliding was not controlled and was not included in any of the models. An axisymmetric 163 rectang ular model was developed to simulate the cylindrical CC testing geometry with restricted exterior 164 faces to simulate the impe rmeable steel well (Figure 1D). Any initial lateral pressure exerted by the walls 165 of the CC well during insertion and prior to loading were small compared to the loading experienced by the 166 sample during tests. Thus, it was deemed appropriate to model the CC well as a fixed boundary with no 167 pre-stress. 168 A quasi -linear hyper -viscoelastic material formulation was chosen to model the behavior of skeletal muscle 169 subject to both CC and UC [3], [23], [26], [38] . The model utilized a Yeoh form [41] of a polynomial 170 hyperelastic strain energy density function Ψ(𝐂) (Equation 5). The initial shear modulus 𝐺0 and bulk 171 modulus 𝐾0 are given according to the 𝑁=1 hyperelastic material parameters (Equation 6) [38], [4 1]. 172 Ψ(𝐂)=∑ 𝐶𝑖03 𝑖=1 (𝐼̅1−3)𝑖+∑1 𝐷𝑖2 𝑖=1(𝐽−1)2𝑖 (5) 173 𝐺0=2C10,𝐾0=2 𝐷1 (6) 174 Here 𝐶𝑖0 and 𝐷𝑖 are material parameters that characterize the isochoric and volumetric re sponses, 175 respectively. 𝐼̅1 is defined as 𝐼̅1=𝜆̅1+𝜆̅2+𝜆̅3 where 𝜆̅𝑖=𝐽−1 3 𝜆𝑖 (𝜆𝑖 are the principle stretches), and 𝐽 is 176 the volume ratio. Due to the nonlinearity of the stress -strain curves for UC and CC data, three 𝐶𝑖0 terms 177 and two 𝐷𝑖 terms were used. A Prony series viscoelastic model (Equation 7) was applied to the decoupled 178 responses in Equation 4. Here 𝐾(𝜏) is the time dependent bulk modulus and 𝐺(𝜏) is the time depe ndent 179 shear modulus. 𝐾∞ and 𝐺∞ model long -term bulk and shear moduli, respectively. 𝜏𝑖𝐾 and 𝜏𝑖𝐺 are time 180 constants ( 𝜏1𝐺=𝜏1𝐾=0.05𝑠,𝜏2𝐺=𝜏2𝐾=1𝑠,𝜏3𝐺=𝜏3𝐾=20𝑠,𝜏4𝐺=𝜏4𝐾=400 𝑠) [17], [28], [30], [46] – 181 [48] 182 𝐾(𝜏)=𝐾∞+∑ 𝐾𝑖𝑒−𝜏 𝜏𝑖𝐾4 𝑖=1 𝐺(𝜏)=𝐺∞+∑ 𝐺𝑖𝑒−𝜏 𝜏𝑖𝐺4 𝑖=1 (7) 183 The fin ite element model consisted of thirteen parameters (five hyperelastic, eight viscoelastic). Model 184 calibration (determination of parameters) was again achieved in two steps: a Monte Carlo simul ation 185 followed by a nonlinear least -squares deterministic optimi zation ( lsqnonlin in MATLAB) [27], [ 28], [30], 186 [42]–[45]. For computational efficiency, the Monte Carlo simulation made 5000 random guesses for the 187 eight viscoelastic parameters (all hyperelastic parameters set to 1) and the error function (Equation 8) was 188 used to calculate the weighted d ifference between normalized model and experimental fast -compression 189 stress relaxation data for both UC and CC models simultaneously . This funct ion assigns greater weight 190 around the peak region where there are fewer time points, thus improving the fit thro ughout the optimization 191 procedure. The set of bulk -viscoelastic parameters minimizing the CC error and set of shear -viscoelastic 192 parameters mini mizing UC error were used as initial guesses for the deterministic optimization. Initial 193 guesses for all hyperel astic parameters were set to a value of one. Following calibration, UC and CC slow - 194 compression data were predicted by this optimized model as a means of validation. Due to the simple and 195 symmetric model geometries in this work, mesh convergence analysis sh owed virtually no difference in 196 model behavior as a function of element size. 197 𝑒𝑟𝑟𝑜𝑟 =∑ 𝑡∗(𝑃𝑚𝑜𝑑𝑒𝑙 −𝑃𝑒𝑥𝑝)𝑡𝑝𝑒𝑎𝑘 𝑡=0+∑(𝑃𝑚𝑜𝑑𝑒𝑙 −𝑃𝑒𝑥𝑝) 𝑡401 𝑡=𝑡𝑝𝑒𝑎𝑘 (8) 198 A second set of UC and CC finite element models wi th the same constitutive formulation were calibra ted 199 using only the UC fast compression data. This is to reflect the approach of assuming near -incompressibility 200 with a single -parameter bulk hyperelastic term, as is most common in finite element models of s keletal 201 muscle [23], [25], [33], [49], [50] . The volumetric parameters ( 𝐷𝑖) were assumed to be three, four, and five 202 orders of magnitude larger than the isochoric parameters ( 𝐶𝑖0) to reflect a range of assumpti ons. The time 203 dependent bulk ( 𝐾𝑖) and shear ( 𝐺𝑖) moduli were assumed to the same. This model repres ents the typical 204 approach for finite element modeling of skeletal muscle and was later used to predict the CC fast 205 compression data. 206 Finally, a semi -confined compression model (SC) was developed by surrounding the quarter -brick UC 207 model with a generic line ar elastic material to create a quarter disk (Figure 2 A). The outer boundar y of this 208 disk was restricted laterally , thus while the whole structure was su bject to confined compression, 209 modulating the Young’s modulus of this outer material enabled a semi -confined state for the muscle 210 geometry. This modulation then enabled the simulation of a transition from unconfined (Figure to 2B) to 211 confined (Figure 2C) compression during a 40% compressive strain ramp to mimic the UC model strain . 212 As the Young’s Modulus of the disk was varied, the S C model stress and volume ratio from the peak were 213 recorded. This parametric analysis was performed for the optimized parameters and parameters acquired 214 from the three models calibrated using UC fast compression data. 215 216 Figure 2: Finite element geometry an d boundary conditions of semiconfined compression . The muscle 217 component used the same unconfined quarter brick employing symmetry. The surrounding material 218 completes a quarter -disk of the same height ( 7 mm) and a radius of 20 mm. 219 To investigate the variabi lity of the experimental data collected, the finite element model was fit to the 220 upper and lower bounds of the standard deviation curves. In short, the full model optimized parameters 221 were scaled up or down to match experimental data in four additional cas es: both UC and CC plus one 222 standard deviation (UC+/CC+), both UC and CC minus one standard deviation (UC -/CC-), and the two 223 remaining cases with one plus standard deviation and one minus standard deviation( UC+/CC - and UC - 224 /CC+). The resulting parameters p rovide insight into how the variability of the data presented here affect 225 the observed incompressibility of passive skeletal muscle. Specifically, the initial Poisson’s ratio 𝜈 was 226 determined for each case based on the initial shear and bulk moduli (Equat ion 9) [51]. 227 𝜈=3𝐾0 𝐺0−2 6𝐾0 𝐺0+2 (9) 228 2.5 Statistical Analysis 229 All statistical comparisons between groups were performed using two -sample two -tailed t-tests, with 230 significance set to p<0.05. The goodness of fit (GoF) for all fits were evalu ated with 231 the goodnessOfFit function in MATLAB (Equation 10) [30]. Here 𝑃𝑖𝑚𝑜𝑑 and 𝑃𝑖𝑒𝑥𝑝 are the model and 232 experimental stress values, respectively, at the ith data point and N is the total number of data points. Fits 233 range from -∞ (worst) to 1 (perfect). The overall percent error, peak stress percent error, an d normalized 234 root mean square error (NRMSE) were also determined (Equation 1 1) [30]. 235 𝐺𝑜𝐹 =1−∑ [𝑃𝑖𝑚𝑜𝑑−𝑃𝑖𝑒𝑥𝑝 𝑃𝑖𝑚𝑜𝑑−𝑚𝑒𝑎𝑛 (𝑃𝑖𝑒𝑥𝑝)]2 𝑁 𝑖=1 (10) 236 𝑁𝑅𝑀𝑆𝐸 =√∑ (𝑃𝑖𝑚𝑜𝑑−𝑃𝑖𝑒𝑥𝑝)2 𝑁 𝑖=1 𝑚𝑒𝑎𝑛 (𝑃𝑒𝑥𝑝) (11) 237 238 3 Results 239 Despite lower strain levels, muscle in CC exhibited stiffer behavior than muscle in UC in both fast and slow 240 compression. Muscle in CC showed ~1200% higher peak modulus (p-value<0.001) (Figure s 3-4) (mean 241 UC fast peak modulus = 0.0524 ± 0.0340 MPa, mean CC fa st peak modulus = 1.856 ± 0.908 MPa). Muscle 242 in CC slow compression exhibited a ~860% higher peak stress than muscle in UC slow compression (p - 243 value<0.001) (Fig ures 3-4) (mean UC slow peak modulus = 0.041 ± 0.020 MPa, mean CC slow peak 244 modulus = 1.058 ± 0. 623 MPa). 245 246 Figure 3: Average experimental stress relaxation curves with standard deviation for (A) unconfined fast 247 (solid red) and slow (solid blue) compression, an d (B) confined fast (dashed red) and slow (dashed blue) 248 compression . 249 250 Figure 4: Peak m oduli of muscle samples for the four testing co nditions and t -test p -values. Unconfined 251 compression is represented by circles and confined compression by triangles. 252 253 Muscle in CC and UC fast relaxation showed different relaxation behavior (Figure 5 A). The thre e relaxation 254 ratios (RR1 for 1 -5s, RR2 for 6 -105s, and RR3 for 106 -400s) showed that differences in CC and UC fast 255 compression time dependence was more apparent at the early stage of relaxation. RR1 for UC 256 (0.648±0.043) was approximately twice RR1 for CC ( 0.344±0.191) (p<0.001). RR2 for UC (0.519±0.076) 257 was larger than RR2 for CC (0.403±0.180) (p=0.029), while RR3 was not different (p=0.52 1) (Figure 5B). 258 259 Figure 5: (A) Average normalized experimental stress relaxation curves for the hold phase plotted on 260 logarithmic scales for unconfined (solid red) and confined (dashed red) fast compress ion. (B) Relaxation 261 ratios RR1-3 for unconfined and confined fast compression data and t -test p -values. Unconfined 262 compression is represented by circles and confined compre ssion by triangles. 263 The global stochastic Monte Carlo simulation in conjunction with the deterministic optimization yielded 264 excellent fi ts between the three -term linear Prony series viscoelastic model and normalized experimental 265 stress data (average percen t error = 1.06 ± 0.13%, average NRMSE = 0.016 ± 0.002, average GoF = 0.99 8 266 ± 0.0002). Comparisons of fast compression viscoelastic param eters showed 𝐸1 for UC (0.902±0.036) was 267 larger than that for CC (0.604±0.192) (p<0.001), 𝐸2 for UC (0.043±0.015 ) was smaller than that for CC 268 (0.103±0.038) (p<0.001), and 𝐸3 for UC (0.033±0.019) was smaller than that for CC (0.141±0.078) 269 (p<0.001) (Fig ure 6A). Additionally, 𝜏1 for UC (0.132±0.029 s) was smaller than that for CC ( 0.331±0.190 270 s) (p<0.001) , 𝜏2 for UC (7.612 ±1.667 s) was smaller than that for CC (12.751±0.427 s) (p<0.001), 271 and 𝜏3 for UC (96.824±20.134 s) was smaller than that for CC (18 3.419±90.805 s) (p<0.001) (Figure 6B). 272 273 Figure 6: (A) Relaxation parameters 𝐸1−3 for unconfined and confined fast compression. (B) Time 274 constants 𝜏1−3 for unconfined and confined fast compression. Unconfined compression is represented by 275 circles and confined compression by triangles. 276 277 The global stochastic Monte Carlo in conjunction with the deter ministic optimization again yielded strong 278 concurrent fitting between the finite element model and experimental UC and CC fast data (Table 1) (Figure 279 7A-B). The initial shear and bulk modulus were calculated using the optimized hyperelastic parameters in 280 Table 2 as 0.0445 kPa and 18.89 kPa, respectively [52]. The model also exhibited very strong predictions 281 for non -linear ramp and relaxation for the slow compression data for both UC and CC (Table 1) (Figure 282 7C-D). 283 284 Figure 7: Average fast experimental data with standard deviation, finite element model calibrations, and 285 fits to the standard deviation curves for (A) unconfined compression (experiment in solid red, model in 286 solid black) and (B) confined compression (experiment in dashed red, model in solid black). Average 287 slow experimental data with standard deviation and finite elemen t model predictions for (C) unconfined 288 compression (experiment in solid blue, model in solid black) and (D) confined compression (experiment 289 in dashed blue, model in solid black). 290 Table 1: Overall percent error, peak error, normalized root mean square erro r (NRMSE ), and goodness of 291 fit (GoF) values for finite element model s calibrated to fast unconfined and confined compression data 292 concurrently and validated against slow unconfined and confined compression data concurrently . 293 Error Type Model Type UC CC Percent er ror Calibration 3.6% 5.9% Validation 12.1% 14.1% Peak error Calibration 0.8% -0.4% Validation 11.2% 35.4% NRMSE Calibration 4.6% 5.6% Validation 14.3% 16.8% GoF Calibration 0.99 0.96 Validation 0.88 0.64 294 Table 2: Hyperelastic and viscoelastic parameters of the finite element model calibrated using unconfined 295 and confined fast compression data concurrently. 296 Parameter Type Parameter Symbol Parameter Value Hyperelastic (MPa) 𝐶10, 𝐶20, 𝐶30 2.23e -05, 1.28e -04, 2.52e -05 Hyperelastic (MPa-1) 𝐷1, 𝐷2 105.9, 0.839 Shear Coefficients (-) 𝐺1, 𝐺2, 𝐺3, 𝐺4 0.741, 0.086, 0.093, 0.061 Bulk Coefficients (-) 𝐾1, 𝐾2, 𝐾3, 𝐾4 0.563, 0.150, 0. 108, 0.147 297 The finite element model with the 𝐷𝑖 parameters derived from the 𝐶𝑖0 and the bulk and shear coefficients 298 equal gave strong fits to the UC data (Table 3). However , these models provided poor prediction s for the 299 CC fast data , specifica lly missing the non -linear ramp, stiffness, and the relax ation behavior (Figure 8A- 300 B). The peak stress was seen to increase rapidly with decreasing volume ratio for the semiconfined ( SC) 301 model derived from the full optimized parameters (Figure 8C ). However , the peak stresses for the three, 302 four, and five orders of magnitude parameters increased more linearly and were thus comparatively too stiff 303 at low strains or too soft at high strains (Figure 8C-D). 304 Table 3: Overall percent error, peak error, normalized root mean square error (NRMSE) , and goodness of 305 fit (GoF) calibration values for finite element models fit with fast unconfined compression data only. 306 Error Type 3 orders of magnitude 4 orders of magnitude 5 orders of magnitude Percent error 5.0% 3.9% 5.1% Peak error -9.8% -0.2% 2.7% NRMSE 5.8% 4.3% 5.3% GoF 0.98 0.99 0.98 307 308 Figure 8 : Average confined compression fast experimental data with standard deviation (dashed red) and 309 finite element model predictions (solid gray) for (A) the three, four, and f ive orders of magnitude 310 predictions, and (B) adjusted y -axis to show only experiment and three orders prediction. (C ) Peak stres s 311 versus volume ratio for semiconfined compression models calibrated using unconfined data only (gray 312 curves ) and semiconfined c ompression model calibrated using both unconfined and confined compression 313 data (black squares) . (D) The same model results presented in ( C) shown at lower strain and stress values. 314 The four sets of bulk and shear moduli calculated from the four new fits ( Figure 7 A and B) show 315 the spread in the experimentally found moduli values have only a small effect o n initial Poisson’s 316 ratio (Table 4). Even with the variability observed and investigated in this study, muscle exhibited 317 nearly -incompressible behavior. 318 Table 4. Normalized root mean square error (NRMSE) , initial bulk and shear moduli, bulk -to-shear 319 modulu s ratio, and initial Poisson’s Ratio for the additional four fits to the mean experimental data plus or 320 minus standard deviation. 321 Model Type UC NRMSE CC NRMSE Bulk Modulus (kPa) Shear Modulus (kPa) Ratio (Bulk : Shear) Poisson’s Ratio Mean 4.6% 5.6% 18.89 0.045 423.45 0.499 UC+/CC+ 5.4% 6.6% 26.70 0.075 355.64 0.499 UC+/CC - 6.1% 6.2% 10.56 0.014 749.53 0.499 UC-/CC+ 8.9% 6.6% 10.56 0.078 135.63 0.496 UC-/CC- 8.8% 6.3% 26.72 0.014 1911.86 0.500 322 4 Discussion 323 This study aimed to characterize and compare the stress relaxation behavior of skeletal muscle subject to 324 two volumetric boundary conditions (unconfined compression or UC and confined compression or CC) by 325 comparing stress relaxation data and computational models. We chose to employ both UC and CC 326 conditions as in vivo muscle is likely to experience semi -confined compression (SC) where the muscle is 327 partially restricted by surrounding hard and soft tissues. While previous viscoelastic approaches 328 successfull y modelled stress relaxation of muscle in UC alone [3], [4] , this study simultaneously fit stress 329 relaxation data of muscle in CC and UC, thus providing greater accuracy for the vo lumetric behavior of 330 skeletal muscle. 331 4.1 Experimental Findings 332 Van Loocke et al. found that muscle exhibited a Cauchy stress of ~4 kPa in UC in cross -fiber direction at 333 strain of 30% applied at quasi -static rate of 0.05% s-1 [3]. This stress value is comparable to the First Piola - 334 Kirchhoff stress values from UC following relaxation (<5 kPa). Palevski et al. measured the short term and 335 long-term shea r modulus of porcine gluteus muscle in vitr o by rapid indentation tests [53]. They found the 336 shear modulus to be ~700±300 Pa. The hyperelastic finite element model developed in this study yielded a 337 shear modulus of 445 Pa, which agrees well with these previously publishe d data. A later Van Loocke et 338 al. study used a non -linear viscoelastic model and found the shear modulus to be 523 Pa, which is also in 339 close agreement with the shear modulus found in t his study [4]. 340 The linear viscoelastic model shows that there is approximately 97% relaxatio n of muscle instantaneous 341 modulus associated relaxation (Figure 6) for UC. This finding agrees with the find ings of Van Loocke et 342 al., whose viscoelastic model showed that muscle experiences about 80% relaxation in first 100 seconds of 343 the hold ph ase of st ress relaxation tests. Wheatley et al. also found that muscle exhibits up to 99% relaxation 344 in UC at long relaxation times, and also found muscle exhibits initial relaxation of ~75% , which is 345 consistent with what our viscoelastic model shows (Figure 6) [28]. On the other hand, muscle in CC shows 346 nearly half of initial relaxation that muscle in UC exhibited (Figure 5). This is supported by statistical 347 compari sons of viscoelastic paramet ers (Figure 6). In all, the time dependent behavior of skeletal muscle 348 has been shown to depend on the loading condition. 349 Muscle has been known to be ~70 -80% incompressible fluid [54] and some studies suggest this fluid plays 350 a significant role in the mechanical properties of the tissue [55]–[58]. In UC, when muscle is compressed, 351 fluid is free to redistribute within the tissue while maintaining a nearly constant volume and exude from the 352 sample from the pores on the lateral sides. Thus, it is unlikely to directly bear a significant portion of t he 353 load, leaving solid muscle constituents such as the extracellular matrix and myofibrils to perhaps resist 354 compression directly. When load is applied to muscle in CC, the fluid cannot escape laterally, and must 355 flow through the porous indenter. Fluid cou ld then be retained in interstitial and intracellular space, thus 356 pressurizing and supporting a greater load than in UC. Since fluid is incompressible, this effect can lead to 357 drastic differences in observed stiffness between UC and CC such as those observ ed here . This fluid 358 pressurization may also contribute the differences in relaxation behavior between UC and CC, as muscle 359 has been shown to have a non -negligible permeability [57]. However, this hypothesis remains untested and 360 future work should be completed to directly investigate fluid pressurization’ s effect on viscoelastic behav ior 361 of muscle in different loading conditions. 362 Two sample geometries were identified for experimentation in this study: transverse oriented cuboids for 363 UC and transverse cylindrical plugs for CC. To minimize tissue damage, hig h profile histology blades and 364 surgical scalpels were used for dissection. These geometries, in accordance with previous literature [59], 365 [60] were used due to the constraints from the experimental apparatus used (Figure 1 A and B). Specifically, 366 the differences in tissue stiffness and required boundary conditions based on testing condition made using 367 a single sample geometry between tests unreasonable. While previous work has shown that sample size can 368 affect the observed compressive modulus of pass ive skeletal muscle [61] sample size and dimension are 369 not likely to explain the major differences in compressive stiffness observed in this study. Additionally, the 370 size of muscle fibers (~50 -100 µm) relative to the sample si ze used here (multiple mm) suggest that the 371 specimens used in this study are representative of bulk muscle tissue. 372 While the stress -time data and standard deviation presented here shows that not all samples exhibited 373 identical passive material properties, this is not uncommon for biological soft tissues and i n particular 374 skeletal muscle [3], [4], [62] . These differences are often explained by natural variability of structure and 375 content of constituents such as collagen in the extracellular mat rix, fiber/fascicle size and organization, and 376 fluid c ontent from animal to animal and muscle to muscle. Despite the fact that clear and consistent 377 structure -function mechanisms in passive skeletal muscle are not fully understood, both the extracellular 378 matrix and muscle fibers are involved in passive load tr ansmission in skeletal muscle [63]–[65]. Other 379 sources of variability may be tissue hydration, although all samples were stored soa ked in phosphate 380 buffered saline prior to testing to limit this effect. 381 One limitation of this work is that muscle samples in UC and CC were compressed to different strain levels. 382 The strain levels of 4 0% in UC and 15% in CC ensured that enough load was ap plied to each sample without 383 damaging the samples. These strain levels were determined through extensive pilot stress -relaxation testing 384 to investigate strain level and tissue damage. Due to the highly soft nature of passive skeletal muscle under 385 unconfine d compression and the relatively stiff response in confined compression, these two strain levels 386 gave more comparable data than similar strain levels would . As muscle has been shown be nonlinearly 387 visco elastic in unconfined compression , the differences in relaxation behavior may vary somewhat with 388 strain. However, the major differences in tissue stiffness between testing conditions suggest different 389 mechanisms that support load under UC and CC . Additiona lly, the effectiveness of bulk and shear Prony 390 series viscoelastic terms employed in this study further support the notion of different mechanisms driving 391 the stiffness and time dependent responses in UC and CC. Future work to test muscle in UC and CC at 392 different strain levels would further clarify how relax ation behavior depends on strain level under these 393 conditions. 394 4.2 Model Findings 395 Two modeling approaches were used in this study: an analytical linear Prony series viscoelastic model and 396 an uncoupled Y eoh/Prony hyper -viscoelastic model. Both of these appr oaches had a similar two -step 397 optimization method, but each model served a different purpose. The relatively simple viscoelastic 398 analytical model generated sets of parameters for each individual sample that could be used to compare 399 relaxation behavior of m uscle between unconfined and confined compression. In comparison, the hyper - 400 viscoelastic finite element model was developed to concurrently characterize the behavior of both testing 401 conditions. Together , these two modeling approaches enabled statistical co mparison of testing conditions 402 as well as comprehensive characterization of tissue stress relaxation behavior. The finite element model 403 was calibrated using the fast compression data and used to predict the slow compression data as the fast 404 data encompasse s a more comprehensive time dependent data set; this is also a common practice in 405 viscoelastic modelling [4], [28] . Moreover, the current study aimed to characterize the larger differences 406 between the testing conditions and not the effect of more specific factors like strain rate, which could be 407 investigated in future studies. 408 Blemker et al. used a decoupled strain energy formulation to model the biceps branchii in which the 409 volumetric or bulk parameter is assumed to be five orders of magnitude larger than the isochoric or shear 410 parameters [33]. Simi larly, Calvo et al. and Grasa et al. take the only volumetric parameter to be between 411 two and three orders of magnitude la rger than the isochoric parameters [25], [66] . In this s tudy we collect 412 volumetric compression data (CC) and span the assumptio ns made by Blemker et al. and Calvo et al to 413 predict the CC data. The pre dictions are very poor (Figure 8), showing that this simple assumption is not 414 appropriate for representing the n onlinear stress -strain behavior unde r highly confined compression. We 415 present a finite element model in which the volumetric and isochoric responses are concurrently optimized. 416 Our models still agree that muscle is what would general ly be considered to be nearly -incompressible 417 (initial shear modulus of muscle is ~3 orders of magnitude smaller than initial bulk modulus), but that the 418 volumetric response is nonlinear as shown by the strong predictions of a two -term volumetric function 419 (Figure 7). We thus reco mmend using a higher -order volumetric term to better characterize compressed 420 muscle. 421 We believe that in vivo muscle experiences loading that is most similar to semi -confined compression. This 422 loading could vary between conditions that approach unconfined o r confined compression among different 423 muscles in the body, or even in different regions on the same muscle. This is supported by previous 424 magnetic resonance imaging of passively stretched human tibialis anterio r that observed volumetric strains 425 as high as 20% in one region and nearly 0% in another [67]. The three models representing literature 426 methods and the model that we concurrently optimized using UC and CC data provide very different results 427 for semi -confinement, particularly regarding nonlinea rity (Figure 8B). This may be a concern as passive 428 stiffness nonlinearity could act as a mechanism to prevent damage of bone and other tissues during high 429 impact loads in vivo . Additionally, the use of a single -term volumetric formulation is likely to be e ither too 430 stiff at low volumetric strains or too soft at high volumetric strains. 431 The study presents a finite element model that can concurrently characterize the unconfined and confined 432 compression conditions. The modeling and optimization approach employed here fit thirteen hyper - 433 viscoelastic param eters. This number of parameters were used because (1) the ramp phase for both 434 unconfined and confined compression are highly non -linear (2) muscle specimens relaxed for a four 435 hundred seconds following compresison at a relatively fast rate (up to 40%/sec) , thus enacting a wide range 436 of time dependence. Future work could reduce the number of optimized or varied parameters by locking 437 Prony terms or coupling terms together similar to previous work [30]. Alternatively, the model provides 438 strong predictions and future efforts could combine direction -dependent, contractile, and tensile mechanics 439 as well to create a more comprehensive model. It has previously been shown that fitt ing models to average 440 experimental data yields different parameters than fitting model individually to tests and then averaging the 441 parameters [68]. Since goal of the study was to characterize the bro ader differences between the testing 442 conditions, the authors found it appropriate to only fit the models to averaged data and acquire one set of 443 parameters. 444 One limitation of the finite element model provided in this st udy is that it cannot capture the com prehensive 445 properties exhibited by passive skeletal muscle such as tension -compression asymmetry, anisotropy, and 446 contractile properties. It is generally assumed , however, that skeletal muscle is primarily compressed in the 447 transverse direction in vivo, thus this is the direction of importance when considering in vivo muscle 448 deformation . A comprehensive model is likely to have quite a large number of parameters, thus increasing 449 the model complexity. However, previous studies [30] have explored how parameter coupling can reduce 450 the number of parameters in a m odel while still utilizing nonlinear optimization and statistical interpretati on 451 of model results. In the future, this approach could be employed to implement a model that not only 452 characterizes differences between muscle in unconfined and confined compres sion but also has other 453 established properties of skeletal muscle. We also cho se not to use a biphasic or poroelastic constitutive 454 approach in this work because v iscoelasticity what is commonly used to characterize skeletal muscle stress 455 relaxation and is a computationally efficient and stable approa ch to modeling time dependence. F luid 456 pressurization , however, requires the solution of an additional condition (either pressure equilibrium or 457 conservation of mass of the fluid) and is more unstable at high str ain rates . Future studies exploring the 458 mechanisms involved in tension and com pression and could employ mechanistic anisotropic models that 459 include components such as tension only fibers and saturating fluid. This type of a model would be greatly 460 beneficia l to the field of passive muscle mechanics. 461 462 5 Conclusion 463 In all, the study found that muscle in CC exhibits stiffer compressive behavior despite lower strains and 464 muscle in UC exhibits greater and faster relaxation. This study also showed that concurrent ly fitting 465 isochoric and volumetric hyper -viscoelastic parameters with these data improves model predictions, and is 466 recommended for cases where semi -confinement is likely. Future work to better understand the mechanisms 467 of force transmission in compressed skeletal muscle would greatly benefit the field. 468 6 Acknowledgements and Conflicts of Interest 469 The authors would like to thank the Program for Undergraduate Research (PUR) at Bucknell University 470 (Lewisburg, PA) for funding . Authors have no conflicts of in terest to report. 471 7