35 Aponeuroses are stiff sheath -like components of the muscle -tendon unit that play a vital role in force 36 transmission and thus locomotion . There is clear importance of the aponeurosis in musculoskeletal 37 function, but there have been relatively few studies of aponeurosis material properties to dat e. The goals 38 of this work were to: 1) perform tensile stress -relaxation tests, 2) p erform planar biaxial tests, 3) e mploy 39 computational modeling to the data from 1 and 2, and 4) perform scanning electron microscopy to 40 determine collagen fibril organization for aponeurosis tissue. Viscoelastic modeling and statistical 41 analysis of stress -relaxation data showed that while relaxation rate differed statistically between strain 42 levels (p=0.044 ), functionally the relaxation behavior was nearly the same. Biaxial testing and associated 43 modeling highlighted the nonlinear (toe region of ~2 -3% strain) and anisotropic (longitudinal direction 44 linear modulus ~50 MPa, transverse ~2.5 MPa) tensile mechani cal behavior of aponeurosis tissue. 45 Comparisons of various constitutive formulations showed that a transversely isotropic Ogden approach 46 balanced strong fitting (goodness of fit 0.984) with a limited number of parameters (five), while damage 47 modeling param eters were also provided. Scanning electron microscopy showed a composite structure of 48 highly aligned, partially wavy collagen fibrils with more random collagen cables for aponeurosis 49 microstructure. Future work to expand microstructural analysis and use t hese data to inform 50 computational modeling would benefit this work and the field. 51 52 1. Introduction 53 Aponeuroses are stiff sheath -like components of the muscle -tendon unit that play a vital role in force 54 transmission and thus locomotion (Eng and Roberts, 2018; Bojsen -Møller and Magnusson, 2019) . While 55 aponeurosis has traditionally been viewed as an extension of the tendon that mechanically acts in series 56 with muscle, this has been shown to not be the case, and has led to improper assumptions regarding the 57 function of aponeurosis (Herzog, 2019) . Various experimental and computational studies have shown a 58 complex relationship between aponeurosis deformation and muscle force , both active and passive 59 (Arellano et al., 2016) . Specifically, incisions of the aponeurosis in lateral gastrocnemius muscles of wild 60 turkeys altered the gearing and force relationship in the muscle (Eng and Roberts, 2018) . Varying the 61 width, length and thickness of the aponeurosis in a 3D muscle model has also shown to have a significant 62 impact on peak stretch location and magnitude of muscle (Rehorn and Blemk er, 2010) . However, the 63 specific contributions that each fibrous tissue has on muscle -tendon unit mechanics is still not fully 64 understood (Herzog, 2017) . 65 66 There is clear importance of the aponeurosis in musculoskeletal function , but there have been relatively 67 few studies of aponeurosis material properties to date (Azizi et al., 2009; Bojsen -Møller and Magnusson, 68 2019; Shan et al., 2019) . Determining the material properties of the aponeurosis will provide evidence of 69 its structure and function to identify how it contributes to the movement of the musculoskeletal system. 70 Azizi et al measured uniaxial tensile mate rial propertie in turkey gastrocnemius aponeurosis , exhibiting 71 the highly anisotropic nature of the tissue due to the aligned collagen structure in the longitudinal or 72 along -fiber direction (Azizi et al., 2009) . While some studies have explored structural mechanical 73 properties such as stiffness (Kitaoka et al., 1994; Scott and Loeb, 1995) , this leaves a distinct lack of data 74 on the materi al properties of aponeurosis such as modulus , especially in comparison to other orthop aedic 75 tissues. While there is some debate regarding the similarities and differences of the aponeurosis and 76 tendon, one of the main physiological differences between the tendon and aponeurosis is that the tendon 77 is free from muscle fascicle s with the exception of the myotendinous junction , while the aponeuros is 78 spreads over a large portion of the muscle belly. It has been shown that during muscle force production , 79 transverse stretch in the aponeurosis exceeds stretch in the longitudinal direction by roughly a factor of 80 four (Azizi and Roberts, 2009) . Therefore, it can be inferred that the aponeurosis is subject to biaxial load 81 during muscle contraction, unlike tendon which is primarily subject to uniaxial load. However, studies of 82 aponeurosis material properties have been completed under uniaxial stre tch only (Azizi et al., 2009; Shan 83 et al., 2019) , despite a clear need to investigate biaxial aponeurosis stretch. 84 85 Additionally, m any soft tissues exhibit time dependent (viscoelastic ) behavior such as stress relaxation 86 (Johnson et al., 1994; Provenzano et al., 2001) , and while aponeu rosis is also viscoelastic (Pavan et al., 87 2011) , literature data are limited . Aponeurosis stores and returns energy during locomotion (Wager and 88 Challis, 2016) , which is affected by tissue viscoelasticity . To the authors ’ knowledge there has been one 89 study on stress relaxation of plantar aponeurosis (Pavan et al., 2011) , which has a distinctly different role 90 and deformation state than typical aponeurosis that acts as a sheath over muscle . While skeletal muscle 91 has shown some degree of nonlinear viscoelasticity – a dependence of stress relaxation behavior on the 92 level of tensil e strain (Wheatley et al., 2016a, 2016c) , tendon is generally observed to be quasi -linear 93 viscoelastic – where stress relaxation is independent of strain level (Pioletti and Rakotomanana, 2000) . It 94 remains unknown if sheath aponeurosis exhibits nonlinear or quasi -linear viscoelastic behavior. Thus, 95 there exists a need to better understand the time dependent material properties of aponeurosis tissue. 96 97 While determining material properties of aponeurosis can provide quantitative insight into muscle -tendon 98 unit function, they are most useful in developing and employing computational models of the 99 musculoskeletal system. Approaches to model muscle -tendon unit s range from classical Hill -type models 100 (Ackland et al., 2012) , more geometrically complex finite element simulat ions (Tsui et al., 2004; 101 Wheatley et al., 2018) , and larger scale simulations of whole human movement (Rajagopal et al., 2016) . 102 These simulations can aid in understanding how a healthy musculoskeletal system functions (Rajagopal et 103 al., 2016) and how disease and impairment can compromise this function (Ong et al., 2019) . To properly 104 employ these modeling techniques, however, we must have a comprehensive understanding of how 105 aponeur osis behaves as a mechanical component . Thus, there exists a strong need to further explore the 106 tissue -level material properties of aponeurosis, specifically in response to biaxial stretch and under stress - 107 relaxation conditions. 108 109 In addition to characterizing and modeling material properties of orthopaedic tissues, microstructural 110 analysis and imaging provides vital insight into tissue morphology . Scanning electron microscopy is a 111 viable method to investigating microstructure in or thopaedic tissues such as tendon (Provenzano and 112 Vanderby, 2006; Herod et al., 2016) and muscle (Mohammadkhah et al., 2017; Sleboda et al., 2020) . 113 However, there is a distinct lack of microstructural analysis through use of scanning electron microscopy 114 for aponeurosis tissue, particularly in relation to s tructure -function relationships. While aponeurosis 115 macrostructure shows highly aligned collagen bundles (Azizi et al., 2009) , we hypothesi ze this is also 116 observed at the microstructural level of collagen fibrils. 117 118 Based on the above, we have identified four aims of this research, with the overarching goal of 119 developing a better understanding of aponeurosis mechanical properties: 120 1) Measure the stress -relaxation behavior of aponeurosis tissue, and determine if it is dependent on 121 strain level (nonlinear viscoelastic or quasi -linear viscoelastic) 122 2) Perform planar biaxial tensile testing to evaluate aponeurosis tissue stress -stretch respon se to 123 multi -axial load 124 3) Employ computational modeling based on the above to develop the appropriate material 125 properties of aponeurosis 126 4) Perform scanning electron microscopy to determine collagen fibril organization and alignment in 127 aponeurosis 128 129 We hypothesize that this approach of materials testing, imaging, and modeling will further elucidate the 130 mechanics of aponeurosis tissue. Future use of these data in computational models of muscle -tendon units 131 can aid in understanding locomotion and neuromusc ular impairments. 132 133 2. Methods 134 2.1 Sample Preparation 135 For all testing, p orcine shoulder samples were obtained from a local abattoir and the triceps brachii 136 aponeurosis (Figure 1A) was dissected from the muscle using standard dissection tools . Long and steady 137 cuts using sharp dissection blades were used to ensure minimal damage to the aponeurosis. As much 138 muscle tissue as possible was removed without tearing or damaging the aponeuros is (Azizi et al., 2009; 139 Shan et al., 2019) . All samples were kept moist by applying PBS throughout the dissection process and 140 before testing (Azizi and Roberts, 2009; Azizi et al., 2009) . Each sample was speckled with graphite 141 powder to enhance digital image correlation strain tracking during testing (Lynch et al., 2003; Luyckx et 142 al., 2014) . Graphite washout was prevented by gently dabbing the sample surface with dry gauze before 143 applying the speckle. Postmortem effects were assumed negligible due to similar assumptions made for 144 tendon (Evans and Barbenel, 1975) . Each sample’s thickness was measured using a dissection microscope 145 by taking images across the entire length of the aponeurosis and averaging all of the measurements 146 (Figure 1B). We differentiated the aponeurosis tissue from the muscle tissue under the microsco pe by its 147 distinctive silver coloration compared to muscles red coloration. To acquire thickness measurements 148 across the entire length of the sample, nine images were obtained, and four measurements were taken off 149 each image with image acquisition software (ThorCam, Thorlabs, Inc.). 150 151 152 Figure 1 . Overview of experimental workflow. A) Triceps brachii muscle with attached aponeurosis 153 (sliver sheath) . B) Light microscope image for thickness measurements (white arrows). C) Dog bone 154 sample , uniaxial material testing instrument, and example digital image correlation image with 155 highlighted region of interest (red rectangle) . D) Cruciform sample , planar biaxial material testing 156 instrument, and example digital image correlation image with with highlighted region of interest (red 157 square) . 158 159 2.1 Stress Relaxation Testing 160 Four porcine shoulders were obtained from a local abattoir to extract the triceps brachii muscle for 161 testing. The aponeurosis samples were cut into dog bone samples approximately 2 2 mm in height by 5 162 mm in width (Figure 1 C) (Azizi et al., 2009) . All samples were cut parallel to the muscle fiber direction 163 (n=20). Tensile stress relaxation was performed on each sample with a 5965 Instron uniaxial testing 164 system with a 5 kN load cell (Figure 1C). Two dimensional i mages were acquired approximately every 165 0.25 seconds throughout each test to determine the precise strain of each sample using a digital image and 166 correlation software (Azizi et al., 2009; Wheatley et al., 2016b) . The employed test procedure was as 167 follows: 20 cycles of preconditioning (between 0.5 N and 1 N), then two consecutive stress relaxation 168 steps of approximately 1 -2% sample engineering strain at a rate of 1 -2%/sec and 300 seconds of 169 relaxation (Komatsu et al., 2007) . Engineering stress (Equation 1) was calculated from force 𝐹 and initial 170 specimen area 𝐴0. Engineering s train for each stress relaxation ramp was determined with digital image 171 correlation over a rectangular region of interest comprising the width of the sample and approximately 5 172 mm in height above and below the centerline (Figure 1C). Stress data were t hen fit to a power law 173 (Equation 2) , which characterized the rate of relaxation with a 𝑏 parameter. 174 175 σ=𝐹 𝐴0 (1) 176 σ=𝑎𝑡𝑏+𝐶 (2) 177 178 2.2 Biax ial Tensile Testing 179 Similar to fascia samples (Eng et al., 2014) , aponeurosis samples were cut into cruciform shapes 180 approximately 30 mm in arm length and 20 mm in arm width by a custom die (Figure 1 D). All samples 181 (n=9) were punched such that the aponeurosis fibers were best aligned with one orientation of the 182 crucifo rm, which was possible due to the highly aligned nature of the tissue and size of the porcine 183 muscle. Samples were speckled with graphite powder similar to stress -relaxation testing (Figure 1D) . 184 Biaxial tensile testing was completed in a planar biaxial mat erial testing system (ADMET, Inc.) with 22 0 185 N load cells (Figure 1 D). An initial gauge length of approximately 30 mm for each orientation was 186 selected, and t he longitudinal direction of the tissue was preloaded to 0.5 N while the transverse direction 187 of the tissue was preloaded to 0.15 N to remove slack . Equibiaxial constant rate tensile testing was then 188 performed at a grip -to-grip strain rate of 0.05%/sec until failure. During testing, images were taken 189 approximately every second and digital image correlation was performed on a 10 mm x 10 mm region of 190 interest centered on the initial image of each sample (Figure 1D). For two -dimensional digital image 191 correlation, longitudinal and transverse stretch were determined for constitutive modeling (Szczesny et 192 al., 2012) (Correlated Solutions, Inc.). Thickness measurements were measured post-hoc from tested 193 samples due to the inability to obtain accurate thickness measurements of the aponeurosis in the region of 194 interest without damaging the tissue . The cruciform samples were cut at 45-degree angles between each 195 arm making an “X” in the middle of the sample. Five measurements were taken from each cut and all of 196 the measurements were averaged together to obtain the overall thickness in the region of interest . Regions 197 with clear damage were avoided for thickness measurements. This method was chosen due to the failure 198 mechanism of aponeurosis tissue (rupture versus plastic deformat ion) and the similarities in measured 199 thickness between stress -relaxation and biaxial sample s. 200 201 Nominal stress in both the longitudinal and transverse direction was determined for each sample using the 202 force values from MTESTQuattro and the thickness and a rm length measurements . The nominal stress 203 value in the arm of each sample was multiplied by a cofactor to determine the stress in the center of the 204 cruciform sample similar to (Labus and Puttlitz, 2016) . The cofactor was determined via a finite element 205 modeling protocol similar to previous modeling approaches of biaxially stretched, soft, fibrous tissues 206 (Jacobs et al., 2013) . A quarter -symmetric cruciform finite eleme nt model with a 25 mm2 square region of 207 interest was developed in FEBio (Figure 2A) (Maas et al., 2012) . The average Cauchy stress data in the x 208 or y direction from the region of interest was compared to the nominal stress in the corresponding 209 cruciform arm. By employing a transversely isotropic neo -Hookean material model, a graph of correction 210 factor v ersus stretch was created to determine the trend between the correction factor values and the 211 various strain levels in both the x and y direction (Figure 2B) (Jacobs et al., 2013) . To determine the 212 appropriate correction factor functions for the longitudinal and transverse directions, piecewise functions 213 were fit to model correction factor data points (Figure 2 B). 214 215 216 217 Figure 2. A) Quarter symmetric finite element model for correction factor modeling. B) Correction factor 218 model outputs (circles) and approximated correction factor piecewise functions (solid curves). 219 220 2.3 Modeling 221 2.3.1 Viscoelastic Modeling 222 To characterize relaxation behavior of the aponeurosis, a three term P rony series quasi -linear viscoelastic 223 model (Equations 3-4) was fit to both a normalized hold phase stress from both the stretches and an 224 average hold phase stress from both the stretches (Chen, 2000) . 225 226 𝜎(𝜀,𝑡)=∫𝐸(𝑡−𝜉) 𝑑𝜀(𝜉) 𝑑𝜉𝑑𝜉𝑡 0 (3) 227 𝐸(𝑡)=𝐸0(1−∑𝐸𝑖[1−exp(−𝑡 𝑖)])3 𝑖=1 (4) 228 229 Here stress 𝜎 is calculated from the convolution integral (Equation 3), which includes the Prony series 230 reduced relaxation function 𝐸(𝑡), nominal strain 𝜀, and an integration parameter 𝜉. This includes three 231 relaxation coefficients 𝐸𝑖, three time constants 𝑖, and the instantaneous modulus 𝐸0. As the purpose of 232 this analysis was to compare relaxation behavior only, 𝐸0=1 was fixed and all data were normalized. 233 Parameter determination was performed in two steps: a Monte Carlo simulation followed by a nonlinear 234 least-squares deterministic optimization ( lsqnonlin in MATLAB) (Vaidya and Wheatley, 2019) . In the 235 Monte Carlo simulation, the six parameters ( 𝐸1−3 and 𝜏1−3) were randomly varied for 100,000 236 simulations, ensuring 0<𝐸1+𝐸2+𝐸3<1. The set of parameter s minimizing percent error between the 237 normalized model and experimental hold stress was used as initial guesses for the deterministic 238 optimization, which optimized percent difference between normalized model and experimental hold 239 stresses. 240 241 2.3.1 Hyperelastic Modeling 242 Adjusted biaxial stress -stretch data were fit to a range of hyperelastic constitutive models through inverse 243 finite element analysis in FEBio (Maas et al., 2012; Wheatley et al., 2015) . A finite element 244 representation of the 10 mm x 10 mm region of interest was d eveloped using quarter symmetry. This 245 simple geometric model was subject to longitudinal and transverse stretch boundary conditions from the 246 mean digital image correlation data. Longitudinal and transverse model stress were output by averaging 247 the Cauchy s tress for each element in each orientation. The model stress -stretch outputs were then fit to 248 experimental mean stress -stretch data by varying model parameters (specific parameters outlined below) 249 within lsqnonlin in MATLAB similar to the viscoelastic modeling approach. A total of five constitutive 250 models were identified and implemented that ranged in nonlinearity, fiber definition, and material 251 damage. For hyperelastic modeling , the constitutive relation is define d by the strain energy density 252 function Ψ (Equation 5, where 𝝈 is the Cauchy stress, 𝐽 is the volume ratio of the solid, 𝑭 is the 253 deformation gradient, and 𝑪=𝑭𝑇𝑭 is the right Cauchy -Green deformation tensor). The strain energy 254 density function for anisotropic materials such as aponeurosis can be decomposed into isotropic and 255 anisotropic components (Equation 6). The isotropic component is typically referred to as an isotropic 256 ground matrix while the anisotropic component generally represents collage n fibers in soft tissues. 257 258 𝝈=2𝐽−1𝑭𝜕Ψ 𝜕𝑪𝑭𝑇 (5) 259 Ψtot(𝑪)=Ψiso(𝑪)+Ψaniso(𝑪) (6) 260 261 For this work, we have identified two isotropic strain energy density terms commonly used for biological 262 and engineered materials subject to large strains, the Mooney -Rivlin formulation (ΨMR) and a first -order 263 Ogden formulation ( ΨOG). The Mooney -Rivlin (MR) model is a function of the first and second 264 invariants 𝐼1 and 𝐼2 of 𝑪 and includes two shear modulus -like parameters 𝑐1 and 𝑐2 (Equation 7). The 265 Ogden model is a summation function of the principal stretches 𝜆𝑖 of the deformation and requires 266 specification of moduli -like parameters 𝑐𝑖 and nonlinear parameters 𝑚𝑖 (Equation 8). For simplicity and 267 reduction of parameters, in this work we have chosen a first -order Ogden model ( 𝑁=1). 268 269 ΨMR(𝐼1,𝐼2)=𝑐1(𝐼1−3)+𝑐2(𝐼2−3) (7) 270 ΨOgden(𝜆1,𝜆2,𝜆3)=∑𝑐𝑖 𝑚𝑖(𝜆1𝑚𝑖 +𝜆2𝑚𝑖+𝜆3𝑚𝑖−3)𝑁 𝑖=1 (8) 271 272 The specific formulations presented in Equations 7-8 are generic forms of these models, and for finite 273 element implementation they must be adapted to include the material response to volumetric 274 deformations. These complete, isotropic coupled equations are provided for the Mooney -Rivlin model 275 (ΨisoMR Equation 9) and Ogden model ( ΨisoOgden Equation 10) and each include a bulk -like modulus 276 parameter 𝑘. This parameter dictates the response to volumetric deformations, as characterized by the 277 volume ratio 𝐽. 278 279 ΨisoMR(𝐼1,𝐼2,𝐽)=𝑐1(𝐼1−3)+𝑐2(𝐼2−3)−2(𝑐1+2𝑐2)ln𝐽+𝑘 2(ln𝐽)2 (9) 280 ΨisoOgden(𝜆1,𝜆2,𝜆3,𝐽)=∑𝑐𝑖 𝑚𝑖(𝜆1𝑚𝑖 +𝜆2𝑚𝑖+𝜆3𝑚𝑖−3−𝑚𝑖ln𝐽)𝑁 𝑖=1 +𝑘 2(ln𝐽)2 (10) 281 282 For material anisotropy, nonlinear tension -only fibers were embedded into the isotropic ground matrix. 283 The constitutive relation for these fibers (Ψanisof, Equation 11) is a piecewise function that enacts this 284 tension -only response and is a function of the fiber stretch 𝜆𝑓. The material parameters for this relation 285 include a transition stretch 𝜆𝑓∗ where the fiber transitions from a nonlinear po wer law response to a linear 286 response, a power law parameter 𝛽 and a fiber modulus 𝐸. It should be noted that parameters 𝐴, 𝐵, and 𝐶 287 are calculated to enforce continuity of strain energy density. This fiber material model can be included in 288 the overall material response in FEBio using the solid mixture material (Maas et al., 2012) . 289 290 Ψanisof(𝜆𝑓2)= { 0 ,𝜆𝑓≤1 𝐴 2𝛽(𝜆𝑓2−1)𝛽,1<𝜆𝑓<𝜆𝑓∗ 𝐸(𝜆𝑓∗−𝜆𝑓)+𝐵(𝜆𝑓2−𝜆𝑓∗2)+𝐶,𝜆𝑓∗≥𝜆𝑓 (11) 291 292 The final constitutive formulation that was explored in this work was reactive damage modeling (Safa et 293 al., 2018) . In continuum damage modeling applications, the strain energy density function is reduced 294 linearly with accumulated material damage 𝐷 (Equation 12). For this work, m aterial damage was 295 specified with respect to two conditions – the initiation of damage and the accumulation of damage. Here 296 initiation was chosen as the normal Lagrange strain for the reinforcing fibers and assuming that damage 297 can only accumulate (this work does not consid er the effects of healing and regeneration), we defined the 298 damage accumulation with a log -normal cumulative distribution function (Equation 13, where Ξ is the 299 damage criterion, in this case normal Lagrange strain, μ dictates the value of Ξ where 50% damag e has 300 occurred, and 𝜎 dictates the rate of damage accumulation ). We also assumed that the maximum amount of 301 damage was 100%, or 𝐷𝑚𝑎𝑥 =1. 302 303 Ψr(𝑪,𝐷)=(1−𝐷)Ψ (12) 304 𝐷(Ξ)=1 2erfc[−ln(Ξ μ) 𝜎 2] (13) 305 306 A total of five different constitutive models were fit to experimental data. These included 1) transversely 307 isotropic Mooney -Rivlin, 2) transversely isotropic Ogden, 3) Mooney -Rivlin with two fiber families, 4) 308 transversely isotropic Ogden with longitudinal fiber damage, and 5) Moo ney-Rivlin two fiber families 309 with longitudinal damage. For all models the bulk -like modulus was fixed at 𝑘=100 MPa. For Mooney - 310 Rivlin models, a 𝑐1=𝑐2 constraint was applied to ensure unique solutions for optimization fitting. 311 312 2.4 Statistical Analysis 313 The experimental stress relaxation data were analyzed by comparing the power law coefficient 𝑏 that 314 resulted from the first and second relaxations through a paired two sample t -test for means ( significance 315 set at p<0.05) to determine whether the relaxation after the first and second strains of each sample were 316 statistically significant . 317 318 The goodness of fit (GoF) for all fits were evaluated with the goodnessOfFit function in MATLAB 319 (Equation 15). Here 𝑃𝑖𝑚𝑜𝑑 and 𝑃𝑖𝑒𝑥𝑝 are the model and experimental stress values, respectively, at 320 the ith data point and N is the total number of data po ints. Fits range from -∞ (worst) to 1 (perfect). The 321 overall percent error, peak stress percent error, and normalized root mean square error (NRMSE) were 322 also determined (Equation 16). 323 324 𝐺𝑜𝐹=1−∑[𝑃𝑖𝑚𝑜𝑑−𝑃𝑖𝑒𝑥𝑝 𝑃𝑖𝑚𝑜𝑑−𝑚𝑒𝑎𝑛(𝑃𝑖𝑒𝑥𝑝)]2 𝑁 𝑖=1 (15) 325 𝑅𝑀𝑆𝐸=√∑(𝑃𝑖𝑚𝑜𝑑−𝑃𝑖𝑒𝑥𝑝)2𝑁 𝑖=1 𝑁 (16) 326 327 2.5 S canning Electron Microscopy 328 Two porcine samples were obtained from a local abattoir to extract the triceps brachii aponeurosis for 329 scanning electron microscopy ( SEM ). The triceps brachii aponeurosis was dissected by the same 330 procedure as the stress relaxation and biaxial tests, and s ubject to a fixing a preparation protocol for 331 imaging (Carson and Cappella no, 1997) . The samples were cut into approximately 4 mm x 8 mm 332 rectangle s then fixed in 10% formaldehyde for 24 hours at 4 degrees Celsius with a 1:20 volume to 333 fixative ratio. Samples were then subject to a washing alcohol series that included 15 minu tes in 70% 334 ethanol, 15 minutes in 90% ethanol, and then a series of 15 , 30, then 45 minutes in 100% ethanol. To 335 prepare for dehydration, samples were transferred to a 1:2 solution of HMDS to ethanol for 20 minutes 336 followed by 20 minutes in a 2:1 HMDS to ethanol solution, and finally two cycles of 20 minute s in 100% 337 HMDS solution. Appr oximately half of the HMDS in each vial was then discarded and the vials were 338 placed without caps in a vacuum for desiccation . Samples that were dehydrated could then be properly 339 sputter coated with gold approximately 48 to 72 hours late r. Scanning electro n microscopy was conducted 340 on a Hitachi SU5000 Field Emission Gun microscope. Relevant images were analyzed with the 341 Directionality plugin found in Fiji (Schindelin et al., 2012) to determine collagen fibril alignment. The 342 number of bins was set to nineteen and the histogram range was set to 0 -180°. 343 344 3. Results 345 3.1 Stress Relaxation 346 The average thickness of the 20 aponeurosis samples were 0.506 mm ± 0.034 (standard error) . The 347 average power law constant of the 20 samples during the first relaxation phase was -0.1013 ± 0.0 057 348 (standard error ) and during the second relaxation phase the average power law constant was -0.0893 ± 349 0.0018 (standard error ). The paired t -test compari ng the first and second relaxation phases were 350 statistically significant (p=0.044). However, further evaluation of 𝑏 parameter versus strain percentage 351 data shows minimal variation between the first relaxation (Ramp 1) and the second relaxation (Ramp 2) 352 (Figure 3 ), as an R2 value of 0.0291 with a linear slope of 0.0026 was fit to these data . Thus, any strain 353 dependent relaxation may be statistically significant, but minimal. 354 355 Figure 3. Strain v ersus power law coefficient for each relaxation phase with a linear regression fit. 356 357 Comparison between the experimental viscoelastic behavior data and the viscoelastic model shows 358 similar comparable behaviors (Figure 4). The global stochastic Monte Carlo simulation in conjunction 359 with the deterministic optimization yielded excellent fits between the three -term linear Prony series 360 viscoelastic model and normalized experimental stress data (average percent error = 0.29 ± 0.0 1%, 361 average RMSE = 0.40 ± 0.02, average GoF = 1.00) (Figure 4) . For the viscoelastic model, the parameters 362 show a similar amount of stress relaxation behavior (𝐺) and similar rate of relaxation (𝜏) between the 363 first, second, and average relaxation s (Table 1) , suggesting that the use of a linear or quasi -linear 364 viscoelasticity is appropriate for aponeurosis tissue . Additionally, aponeurosis exhibited greater relaxation 365 during initial relaxation as ev ident by the larger 𝐺1 parameters in comparison to 𝐺2 or 𝐺3. 366 367 368 Figure 4 . A) Relaxation behavior of average experimental data and model data for (Ramp1, Ramp2, and 369 ModelAvg) over the entire 300 second relaxation periods . B) Highlighted r elaxation behavi or over the 370 first three seconds of relaxation . 371 372 Table 1 . Optimized viscoelastic parameters for the hold phase of stress relaxation after first, second, and 373 average of first and second stretches. 374 Parameter Ramp One Ramp Two Average 𝑮𝟏 [MPa] 0.420 0.426 0.423 𝑮𝟐 [MPa] 0.135 0.127 0.130 𝑮𝟑 [MPa] 0.101 0.094 6 0.0978 𝝉𝟏 [s] 0.542 0.541 0.551 𝝉𝟐 [s] 8.59 8.83 8.86 𝝉𝟑 [s] 98.0 101 100 375 3.2 Biaxial Stretch 376 Tensile planar equibiaxial stress -stretch data shows consistent anisotropy of aponeurosis tissue (Figure 5). 377 The average thickness of the nine aponeurosis samples were 0.626 mm ± 0.028 ( standard error ). 378 Qualitatively, the aponeurosis is shown to be stiffer in the longitudinal direction compared to the 379 transverse direction and s train more in the transverse direction compared to longitudinal (Figure 5A) . 380 Quantitatively, these results are supported by the finding that over the course of the experiment the stress 381 ratio leveled out at approximately ten and the strain ratio leveled out at approximately 0.5 between 382 longitudinal and transverse orientations (Figure 5B) . 383 384 Figure 5. A) Planar equibiaxial tensile stress -stretch data. Raw data are provided in light background 385 curves for longitudinal (black) and transverse (blue) direction . Mean and standard error are shown with 386 open circles and error bars. B) Mean data ratios for stress (blue) and strain (o range) between the 387 longitudinal and transverse directions as a function of data points. Data presented are 388 longitudinal:transverse ratios. 389 390 Finite element model fits to experimental data were generally strong as observed both visually (Figure 6) 391 and statistically ( Table 2). The model with the worst fit was the transversely isotropic Mooney -Rivlin 392 formulation due to the linear stress -stretch response in the transverse direction (Figure 6B). The best fit 393 was the two fiber families model with longitudinal fiber damage, which produced a root mean square 394 error of 0.00655 MPa . However, this model was statistically nearly identical to the transversely isotropic 395 Ogden damage model, and the non -damage versions of these m odels had strong fit statistics as well 396 (above a 0.95 goodness of fit value (Table 2)). 397 398 399 Figure 6. Comparison between biaxial experimental data and model in A) longitudinal direction and B) 400 transverse direction . Note the significant model overlap. 401 402 Table 2. Statistical evaluations of model fits to experimental data. 403 Model Percent Error Goodness of Fit RMSE [MPa] TI Mooney -Rivlin 87.9 0.945 0.0233 TI Ogden 33.7 0.952 0.020 4 Two Fiber Families 25.2 0.953 0.0198 TI Ogden Damage 24.3 0.984 0.00665 Two Fiber Damage 19.0 0.984 0.00655 404 Model parameters for all five formulations (Table 3) show a longitudinal fiber modulus value of 405 approximately 50 MPa and a transverse modulus of approximately 2.4 MPa. The 𝜆∗ values also suggest 406 that the longitudinal fibers are fully recruited at stretch values between 1.02 and 1.03, while the transverse 407 response exhibited a larger toe -region with full fiber recruitment occurring at approximately 1.06 -1.065 408 (Table 3). The transversely isotropic Ogden model exhibited both strong fitting statistics and a low 409 number of fitting parameters (five), while the two fiber families with damage formulation exhibited the 410 best fitting statistics and enables the use of damage modeling, but requires nearly twice the parameters 411 (nine). 412 413 Table 3. Biaxial finite element m odel parameter values. 414 Model 𝒄𝟏 [MPa] 𝒎𝟏 [-] 𝑬𝒇 [MPa] 𝜷 [-] 𝝀∗ [-] 𝝈 [-] 𝝁 [-] TI Mooney -Rivlin 0.0493 - 43.3 5.17 1.0233 - - TI Ogden 0.185 41.0 44.0 4.82 1.023 4 - - Two Fiber Families 0.0172 - 44.3 4.76 1.0237 - - 2.35 5.21 1.0619 - - TI Ogden Damage 0.157 44.2 69.8 4.02 1.0327 0.509 0.118 Two Fiber Damage 0.0135 - 55.1 4.06 1.0296 0.316 0.122 2.37 4.34 1.0643 - - 415 3.3 Scanning Electron Microscopy 416 Scanning electron microscopy (SEM) of aponeurosis at different magnifications shows a composite 417 hierarchical structure with highly aligned collagenous sheets and collagen cables (Figure 7). Superficial 418 collagen cables (Figure 7A) appear to have some degree of random organization, while a collagen fibril 419 sheet (Figure 7B) is comprised of highly aligned collagen fibrils. In addition to variation in structure, 420 collagen fibril sheets appear to have a variation in waviness (Figure 7C). Characterization of fiber 421 dispersion shows that the directional alignment in aponeurosis can very from highly aligned and straight 422 fibers (dispersion standard deviation of ~7 degrees, Table 4) to a more dispersed case (~30 degrees, Table 423 4). 424 425 Table 4. Gaussian distribution fit dispersion (standard deviation of Gaussian fit) and go odness of fit for 426 analyzed scanning electron microscopy images. 427 Image Dispersion [deg] Goodness of Fit [-] Straight Image 6.94 0.97 Wavy Image 27.4 0.95 428 429 Figure 7. Scanning electron microscopy of aponeurosis tissue and analysis of fiber alignment. A) Image 430 showing aligned collagen sheet under collagen cables. B) Higher magnification of aligned fiber sheet with 431 straight fibrils (“Straight”) taken from a region in A . C) Higher magnification of aligned fiber sheet with 432 wavy fibrils (“Wavy”) taken from a region not shown in A . D) Comparison of collagen fibril alignment 433 image analysis data and normal distribution fits for the images in B and C. Note that data are plotted as a 434 function of angle from the transverse direction, or that a value of 90° indicates the longitudinal direction. 435 436 4. Discussion 437 While the broader goal of the study was to elucidate the mechanic al properties of aponeurosis tissue, the 438 four specific goals were: 1) measur e the stress -relaxation behavior of aponeurosis tissue and determ ine if 439 there was strain level dependence (nonlinear viscoelastic or quasi -linear viscoelastic) , 2) perfor m planar 440 biaxial tensile testing to evaluate aponeurosis tissue response to multi -axial load, 3) employ 441 computational modeling based on the above to develop the appropriate material properties of aponeurosis , 442 and 4 ) perform scanning electron microscopy to visually determine collagen organization in aponeurosis. 443 444 4.1 Stress Relaxation 445 When co nsidering both statistical analysis (Figure 3) and viscoelastic fitting (Table 1 and Figure 4), we 446 conclude that from approximately 1 -4% tensile strain , aponeurosis relaxation behavior is comparable . 447 While a statistically significant result between power law parameters for the two relaxation steps was 448 found (p=0.044), this represents a case where the data may be statistically different but not functionally 449 different. These results suppor t the claim tha t the aponeurosis exhibits quasi -linear viscoelastic behavior 450 (Pavan et al., 2014) , similar to other fibrous connective tissues like tendon (Taylor et al., 1990; Johnson et 451 al., 1994; Provenzano et al., 2001) . Qualitative examination of digital image correlation strain plots 452 suggests that t he aponeurosis does not strain in a homogenous manner due to the non -uniform recruitment 453 of collagen fibers and fascicles during uniaxial tensile testing. These findings are similar to a previous 454 study that showed that the aponeurosis -tendon unit strained nonuniformly when submaximal voluntary 455 contractions occurred in vivo (Finni et al., 2003) . While characteriz ation of this inhomogeneous strain was 456 outside of the scope of this work, future efforts to measure regional aponeurosis strain would be a benefit 457 to the field. One limitation of the material tests performed in this work is the lack of tissue immersion in a 458 buffering solution. This was chosen for comparative purposes to other studies of aponeurosis material 459 properties that tested in non -immersed conditions (Azizi et al., 2009; Shan et al., 2019) , to improve digital 460 image correlation analysis, and to introduce as few variables into the testing protocol as various buffer 461 solutions can affect measured mechanical properties (Safa et al., 2017) . However, future work to eva luate 462 the effect of various buffer solutions, graphite powder speckling and dehydration, and of immersed versus 463 non-immersed conditions could further improve the translation of in vitro data to in vivo models and 464 considerations. 465 466 Viscoelastic fitting of the three -term linear Prony series model to the normalized experimental data show 467 an excellent ability of this approach to fit stress -relaxation data . As highlighted above, the similar 468 viscoelastic model parameters optimized for bot h stress -relaxation ramps and the accuracy of a fit to 469 mean data suggest the parameters provided in Table 1 could be reasonably used for future modeling 470 efforts. One l imitation of this work is that only normalized stress relaxation data was used in the 471 viscoelastic model , which was completed for simplicity . Additionally, these viscoelastic parameters were 472 not validated against independent data, and as such they should only be used for simulations of similar or 473 slower strain -rate conditions. Future studies could incorporate and couple the stress -relaxation and stress - 474 stretch data presented in this study to other viscoelastic models and perform validation experiments . 475 476 4.2 Biaxial Stretch 477 To the best of the authors’ knowledge, there exists only two previous studies that have directly measure d 478 stress -stretch material properties of aponeurosis (Azizi et al., 2009; Shan et al., 2019) . Azizi et al ., 2009 479 measured the anisotropic line ar modulus and failure properties of turkey gastrocnemius aponeurosis and 480 found that the longitudinal direction had a higher tensile modulus (~700 MPa) value compared to the 481 transverse direction (~100 MPa) . Shan et al. studied the morphological and mechanical characteristics of 482 the human triceps surae aponeuroses in various regions of the gastrocnemius and soleus , determining that 483 thickness varied depending on location and a similarly higher modulus in the longitudinal ( 64-145 MPa) 484 versus transverse (0.3-1.5 MPa ) directions (Shan et al., 2019) . Our linear region modulus results 485 presented here were generated by constitutive model fitting via inverse finite element analysis. The 486 findings of a longitudinal modulus of ~50 MPa and transverse modulus of ~2 MPa agree well with Shan 487 et al. In comparing these values to tensile moduli for tendon and muscle, we see that aponeurosis tissue 488 compares similarly to tendon and is approximately 1 -2 orders of magnitude stiffer than muscle (Table 5). 489 This suggests that aponeurosis functions from a materials perspective as an extension of tendon and not as 490 a transition from tendon to muscle. However, this requires further extensive material testing of 491 aponeurosis and tendon from the same muscle -tendon unit to confirm. Despite the range of reported 492 moduli values for aponeurosis in literature (Table 5), our findings corroborate the observation that the 493 aponeurosis is stiffer in the longitudinal direction compared to the transverse direction and that the 494 transverse direction experienced more stretch during the experiment (Figure 5) (Azizi and Roberts, 2009; 495 Azizi et al., 2009) . 496 497 498 Table 5. Comparison between anisotropic linear -region tensile elastic modulus values reported in 499 literature and those presented in this study for tendon, aponeurosis, and muscle tissue. Note that * denotes 500 some approximate values for skeletal muscle stress -stretch curve tangents, as muscle exhibits a great 501 degree of nonlinearity. 502 Tissue Longitudinal Modulus [MPa] Transverse Modulus [MPa]