Ligaments and tendons undergo nonuniform deformation during movement. While deformations can be imaged, it remains challenging to use such information to infer regional tissue loading . Shear wave tensiometry is a promising noninvasive technique to gauge ax ial stress and is premised on a tensioned beam model. However, it is unknown whether tensiometry can predict regional stress in a nonuniformly loaded structure . The objectives of this study were to (1) determine the relationship between regional shear wave speed and regional axial stress, and ( 2) determine the sensitivity of regional axial stress and regional shear wave speed measurements to nonuniform load distribution and fiber alignment in the presence of variable tissue geometry and constitutive properties . We created a represent ative set of 12,000 dynamic finite element model s of a fibrous soft tissue with probabilistic variations in fiber alignment, stiffness, and aspect ratio. In each model, we applied a randomly selected , nonuniform load di stribution, and then excited a shear wave and tracked its Blank et al. , 2023 2 regional propagation. We found that regional shear wave speed was an excellent predictor of the regional axial stress (RMSE=0.55 MPa ) and that the nature of the wave speed -stress relationship was co nsistent with a tensioned beam model (R2 = 0.99). Variations in tissue geometry, nonuniform loading, fiber material properties, and fiber alignment did not substantially alter the wave speed -stress relationship, particularly at higher loads . Thus, these finding suggests that shear wave tensiometry could provide for a quantitative estimate of regional tissue stress in ligaments and tendons. 1. INTRODUCTION Ligaments and tendons undergo complex deformation during movement . For example, joint motion can induce nonuniform ligament stretch due to rotation of the ligament ’s attachments on adjacent bone segments (Robinson et al. 2004; Willinger et al. 2020) . Nonuniform deformation is also observed in tendinous tissue which can undergo heterogenous muscle loading (Franz et al. 2015) . The capacity of ligaments and tendons to change shape and effectively transmit loads is critical to normal biomechanical function. Interfibrillar sliding is a structural factor that may further facilitate nonuniform deformation (Szczesny and Elliott 2014 ), and t he degree of fiber alignment can affect the relative amounts of sliding and shear load transmission that occurs. For example, ligament fibers are relatively less aligned than in tendon (Amis 1998; Wymenga et al. 2006) , which is believed to provide for complex load transmission about a joint . Assessing nonuniform loading remains a challenge in biomechanical assessments of soft tissue mechanics. Currently, digital image correlation is the gold standard technique for gauging strain distributions (Mallett and Arruda 2017) . In combination with m aterial -informed finite element models, one can predict stress distributions in tendon and ligament (Readioff et al. 2020) . Such measurements clearly demonstrate the complex mechanics of the tissue in response to an anatomical load (J C Gardiner and Weiss 2003) . However, it remains challenging to verify these model -based stress estimates due to the p aucity of techniques for measuring regional stress . Blank et al. , 2023 3 Shear wave speeds are widely used to study musculoskeletal tissue mechanics (J. Blank et al . 2022) . A tensioned beam models suggests that t he speed of the propagating shear wave is proportional to the square root of axial tissue stress (Martin et al. 2018) . This general relationship is shown to hold for ligamen t (J. L. Blank, Thelen, and Roth 2020a) and tendon (Martin et al. 2019) under pure axial loading . However, it is uncertain whether the wave speed -stress relationship holds in a tissue undergoing nonun iform loading. We previously used a probabilistic, dynamic finite element model to investigate the effects of f iber alignment and material properties on shear wave speeds in ligaments and tendons (J. L. Blank et al. 2021) . This study extends that work by considering nonuniform deformations that can arise in a ligament due to joint rotation. Accordingly, t he objectives of this study were to, (1) determine the relationship between regional shear wave speed and regional axial stress, and ( 2) determine the sensitivity of regional axial stress and regional shear wave speed measurements to nonuniform load distribution and fiber alignment in the presence of variable tissue geometry and constitutive properties. 2. METHODS 2.1. Finite Element Model Overview We implemented a previously validated finite element model capable of simulating dynamic shear wave propagation in a fibrous soft tissue during nonuniform loading (J. L. Blank et al. 2021) (Fig. 1). The tissue wa s modeled as a 60 mm -long ellipse with a predefined width and thickness . The tissue was composed of a mesh with 3,600 hexahedral elements , which was prescribed under guidance of a prior mesh convergence study (J. L. Blank et al. 2021) . The proximal end of the tissue was attached to a rigid body capable of being loaded with unequal tension s on each side to create a nonuniform load distribution (Fig. 1a) . Blank et al. , 2023 4 Figure 1 : (a) We modeled the tissue as transversely isotropic with distributed fibers in the YZ plane. We first nonuniformly loaded the tissue in the static simulation using two point loads on either si de of the proximal rigid body. We then excited a shear wave at the bottom of the tissue and let it propagate in a dynamic simulation. (b) We monitored the transverse tissue motion at each point across the tissue width and length and measured a shear wave speed in each region (12 shear wave speeds total per simulation). 2.2. Constitutive Model We modeled the tissue as an incompressible, transversel y isotropic hyperelastic material with distributed fibers that represents a material with longitudinal fibers embedded in a matrix using the following combination of strain energy density functions (Eq. 1): Ψ=Ψ𝑚𝑎𝑡𝑟𝑖𝑥 +Ψ𝑓𝑖𝑏𝑒𝑟𝑠 (1) Here, Ψ𝑚𝑎𝑡𝑟𝑖𝑥 represents the contribution of the ground matrix , which governs the transverse shear modulus, and Ψ𝑓𝑖𝑏𝑒𝑟𝑠 represents the contribution of embedded fibers. The ground matrix was isotropic according to a Mooney -Rivlin constitutive model (𝐶1=𝐶2= 0.205 MPa) . The ground matrix, embedded fibers, and tissue’s volumetric response to loading were characterized using the following uncoupled strain energy density function, Ψ (Eq. 2) : Blank et al. , 2023 5 Ψ=𝐹1(𝐼̃1,𝐼̃2)+∫ 𝑃(𝜃) 𝐹2(𝜆̃[𝜃]) 𝑑𝜃𝜃𝑝+𝜋 2 𝜃𝑝−𝜋 2+𝐾 2(𝑙𝑛𝐽)2 (2) 𝐼̃1 and 𝐼̃2 represent the first and second invariants of the right Cauchy -Green deformation tensor, 𝜆̃ indicates the deviatoric part of the stretch ratio along the fiber direction, 𝐾 represents the bulk modulus, and finally 𝐽 represents the volume c hange of the deformation . The parameter 𝜃 is indicative of changes in fiber alignment according to the 2D unimodal distribution function, P( 𝜃) (Eq. 3) : 𝑃(𝜃)=1 𝜋𝐼0𝑘𝑓𝑒𝑘𝑓cos (2(𝜃−𝜃𝑝)) (3) The principal fiber orientation, 𝜃𝑝, represents the nominal fiber angle and was chosen to align with the longitudinal axis of the tissue (defined from positive z -axis in YZ -plane, Figure 1). 𝐼0 represents the modified Bessel function of the first kind, and the fiber concentration factor, 𝑘𝑓, determines the degree of fiber alignment within the tissue. Fiber concentration values from 0 to 100 were chosen to represent highly unaligned (i.e., isotropic) and highly aligned (i.e., transversely isotropic) fib ers, respectively (Stender et al. 2018) . 𝐹2(𝜆̃) represents the contribution of the collagen fibers to the strain energy density function (Eq. 4) : 𝜆̃𝜕𝐹2 𝜕𝜆̃=0 𝜆̃≤1 𝜆̃𝜕𝐹2 𝜕𝜆̃=𝐶3𝑒𝐶4(𝜆̃−1)−1 1<𝜆̃≤𝜆∗ 𝜆̃𝜕𝐹2 𝜕𝜆̃=𝐶5𝜆̃+𝐶6 𝜆̃>𝜆∗ (4) Here, 𝜆̃ represents the tissue stretch during loading, and 𝜆∗ is the stretch at which the fibers engage. Prior to this stretch level, the fibers are either slack or uncrimpin g, where the term C3 scales the exponential stress and the term C4 controls the strain -dependent rate at which the fibers uncrimp. The Blank et al. , 2023 6 term C5 corresponds to the elastic modulus of straightened fibers and C6 is determined using the requirement that the str ess-strain curve be continuous at 𝜆∗. For a transversely isotropic material, the stretch of the collagen fibers is the same as 𝜆̃, which is a function of deformation in the z -direction. For tissues with unaligned fibers, this response is scaled using the 2D unimodal distribution function . 2.3. Probabilistic Model Parameters We generated 1 2,000 models with variable tissue geometry (§2.3.1.), nonuniform loading (§2.3.2.), and fiber stiffness and alignment (§2.3.3.) as a representative cohort for variations in musculoskeletal soft tissues in human patients (Table 1). Table 1: Ranges of parameters varied in our cohort of models with probabilistic variations in geometr y, boundary conditions , hyperelastic material properties , and microstructure . Probabilistic M odel Category Parameter description Parameter Value Tissue geometry Tissue thickness w [mm] 3-6 Tissue width t [mm] 6-30 Nonuniform load Nonuniform load distribution ΔF [N] 0-100 Axial load F [N] 0-600 Fiber stiffness Toe region stress term C3 [MPa] 0.1 – 1.5 Toe region uncrimping term C4 30 – 70 Linear region modulus C5 [MPa] 300 - 900 Fiber alignment Fiber alignment factor kf 0-100 2.3.1. Tissue Geometry Parameters We evaluated the sensitivity of load distribution and regional shear wave speeds to tissue geometry by varying the tissue’s thickness from 3 to 6 mm and the widt h from 6 to 30 mm according to a random uniform distribution. The tissue was held to a constant width and thickness along the entire length. 2.3.2. Nonuniform Load Distribution Blank et al. , 2023 7 We generated nonuniform loads to evaluate the capability of regional shear wave speed to track regional axial stress during different loading conditions . We varied the axial load on the either side of the model from 0 to 300 N according to a random uniform distribution . We ensured that the difference in load between the two sides of the tissue did not exceed 100 N. 2.3.3. Fiber Alignment and Fiber Stiffness We systematically varied fiber alignment and stiffness to assess their effect s on regional loading and the resultant regional shear wave propagation. We varied fiber alignment by varying the value of k f from 0 to 100 according to a positive, one -sided norm al distribution (mean = 0, standard deviation = 25). Our one -sided normal distribution was selected to reflect the non -linear decreasing effect of fiber alignment on mechanics as k f is increased (J. L. Blank et al. 2021) . We also varied fiber stiffness using the constitutive model constants C3 (toe-region exponential stress coefficient), C4 (rate of fiber uncrimping), and C5 (linear region modulus). Fiber material properties were varied according to a random uniform distribution (Table 1) . 2.4. Simulation of Shear Wave Propagation under Nonuniform Loading We simulated shear wave speeds using a previously validated finite ele ment (FE) modeling framework (J. L. Blank, Thelen, and Roth 2020b; J. L. Blank et al. 2021) . All finite element simulations were performed in FEBio version 3.0 (Maas et al. 2012) . The bottom surface of the structure was held fixed. In the initial, static portion of the simulation, we loaded the rigid grip with unequal tension s on both sides by applying a predefined force to one node on each end of the grip (Fig. 1). During the dynamic portion of the simulation, we excited a shear wave in the model by applying a transverse displacement ramp of 10 microns at 2 kHz to nodes throug h the model’s cross -section (Fig. 1) . The excitation displacement and timing were chosen to match the excitation signal of current voice coil - Blank et al. , 2023 8 actuated tensiometer designs. Following the excitation, the shear wave was allowed to propagate along the structur e for 2.5 ms. We executed all our simulations on a high -throughput computing grid . We measured regional shear wave speeds according to prior finite element modeling studies (J. L. Blank et al. 2021) . Briefly, we monitored transverse tissue motion at locatio ns both along the length and width of the tissue (Fig. 1). The shear wave travel time was determined using a normalized cross - correlation of the transverse displacement profiles two consecutive measurement points (1 cm distance) along the structure’s long -axis. C osine interpolation was used to assess sub -sample alignment of the displacement profiles (Cespedes et al. 1995) . Transverse wave profiles having a correlation strength of less than 0. 75 were excluded due to signal decorrelation . Shear wave speed was then computed as the distance between measurement nodes divided by the wave travel time . We measured the corresponding regional axial stress as the average axial Cauchy stress at elements betw een the two measurement points. Shear wave speeds less than 15 m/s or greater than 160 m/s were considered out of range of our tensioned beam model and were excluded. We retained over 96.5% of the 12,000 simulation s for our analyses. 2.5. Statistical Analy sis For our first objective, we assessed the relationship between the regional shear wave speed and regional axial stress using a linear regression according to the following analytical relationship (Martin et al. 2018) (Eq. 5) : 𝜎=𝜌𝑒𝑓𝑓𝑐2−𝑘′𝜇 (5) Here, 𝜎 represents the axial tendon stress (in Pa), 𝜌𝑒𝑓𝑓 the effective density of the propagation mediu m (in kg/m3), 𝑐 the shear wave speed (in m/s), 𝑘′ a shear correction factor (unitless) , and 𝜇 the tangential shear modulus of the tissue (in Pa). We compared shear wave speeds to the tensioned beam Blank et al. , 2023 9 model evaluated using the material properties set for the FE model ( 𝜌𝑒𝑓𝑓 = 1500 , 𝜇 = 0.82 MPa , 𝑘′= 0.89 ), and also used the linear relationship of the tensioned beam model to determine the axial stress given shear wave speed and linear regressi on coefficients. For all linear regressions , the dependent variable was the regional axial stress, σ, and the independent variable was the squared regional shear wave speed, 𝑐2. Errors between the predicted and actual regional stress es are denoted as 𝜀 (in M Pa). To perform our sensitivity analyses stated in our second objective , we first performed a multiple linear regression with one -way interactions between the one of the response variable s (i.e., regional stress distributi on or stress prediction error) and two of the independent variable s (i.e., w, t, ΔF, F, C3, C4, C5, or k f from Table 1 ) (Eq. 6) : 𝑦=𝛼+ 𝛽1𝑥1+𝛽2𝑥2+𝛽3𝑥1𝑥2+𝜀 (6) Here, 𝑦 represents the one of the response variable s (i.e., either regional stress distribution or stress prediction error), 𝑥1 and 𝑥1 represent the two of the independent variables ( i.e., w, t, ΔF, F, C3, C4, C5, or k f ), and 𝛽1−3 represent the regression coefficients that we used to define sensitivi ty for each parameter . We also performed a linear regression between each independent variable (i.e., w, t, ΔF, F, C3, C4, C5, or k f ) and response variable (i.e., regional stress distribution or stress prediction error). Interactions that accounted for more variance in the response variable (i.e., had higher R2) were considered important. 3. RESULTS 3.1. Regional Shear Wave Speeds Predict Regional Axial Stress We observed regional stress distributions that were dependent on the tissue’s fiber alignment. More specifically, less well -aligned fibers corresponded to a more evenly distributed regional stress distribution, and tissues with highly aligned fibers had in creased differences in axial stress magnitudes Blank et al. , 2023 10 across the tissue width (Fig. 2a -b). Regional shear wave propagation speed scaled with increases in regional axial stress (Fig. 2c). Figure 2: (a) Images representing fiber alignment factors (k f) ranging from 0 (random alignment) to 500 (highly aligned) . (b) Axial stress distributions over the width of the structure when subject a 100 -N load differential from left to right and a total axial load of 300 N. (c) Shear wave propagation in response to a transverse ex citation of the bottom of the nonuniformly loaded structure. Wave speed across the width varied substantially as fibers became more aligned. Regional shear wave speed scaled with regional stress. We found that regional shear wave speed squared was an exce llent predictor of the regional axial stress (RMSE = 0.5 5 MPa, R2 = 0.9 9), and that regional shear wave speeds measured in the model agreed well with th ose predicted using the tensioned beam model with applied model parameters ( peff = 1500 kg/m3, μ = 0.82 MPa) (Fig. 3). Correspondingly, regional wave speed was not nearly as good at predicting average tendon stress (σaverage = F/A ), with prediction errors that double d when using regional shear wave speeds (RMSE = 1.11 MPa, R2 = 0.9 3) (Fig. 4 b). Generally higher errors in regional stress predictions were observed for Blank et al. , 2023 11 high nonuniform load distributions in low stress regions. High nonuniform load distributions resulted in higher average stress prediction errors independent of load level. Figure 3: Regional shear wave speed was a better predictor of the regional axial stress (a) (RMSE = 0.5 5 MPa, R2 = 0.9 9) than the average axial stress (b) (RMSE = 1.1 1 MPa, R2 = 0.9 3) in tissues undergoing increasing nonuniform load (as indicated by the colormap ). 3.2. Fiber Alignment Facilitate s Axial Stress Transmission Across the Tissue Our simulations revealed that nonuniform load distribution (ΔF) was an important factor for determining the regional stress distribution (R2 = 0.74) (Fig. 4). Among these factors, we determined that the interaction between nonuniform load distribution and fiber alignment accounted for 85% of the variability in regional stres s under nonuniform loading (Fig. 5) . We also determined that interactions between thickness and nonuniform load distribution accounted for 8 2% of the variability in regional stress (Fig. 3S ). We did not find that other parameters ( i.e., w, F, C 3, and C4) influenced Blank et al. , 2023 12 regional stress distributions by themselves, and that no other interactions increased the sensitivity of regional stress distribution. All sensitivities and interactions are include d in the Electronic Supplementary Materials (Fig. 3S). Figure 4: Sensitivity of regional stress distribution , ∆𝜎𝑟𝑒𝑔𝑖𝑜𝑛𝑎𝑙 to normalized parameter ranges . Shown are the percentage of variability of ∆𝜎𝑟𝑒𝑔𝑖𝑜𝑛𝑎𝑙 accounted for by a single parameter (R2), and the slope of a linear regression fit (in MPa/parameter range). Parameter ranges are given in Table 1. Blank et al. , 2023 13 Figure 5: Interaction between the nonuniform load distribution, ΔF, and the fiber alignment fac tor, kf, resulted in predictable variations in the regional stress distr ibuti on across the tissue wi dth. 3.3. Regional Axial Stress Predictions are Insensitive to Probabilistic Tissue Factors When predicting the regional axial stress using the regional shear wave speed, we found that no individual parameter accounted for more than 4% of the variability in regional stress prediction errors. We found that interactions between geometric terms ( t, w), load terms ( ΔF, F), and fiber alignment ( kf) accounted for more variability in regional stress prediction errors, but no grouping of variables accounted for more than 8% of the v ariability observed. By themselves, fiber alignment and stiffness terms had no effect on regional stress prediction errors (R2 = 0). All sensitivities and interactions are included in the Electronic Supplementary Materials (Fig. 4-5S, Table 1S) . 4. DISCUSSION The objective s of this study were to (1) determine the relationship between regional shear wave speed and regional axial stress, and ( 2) determine the sensitivity of regional axial stress and Blank et al. , 2023 14 regional shear wave speed measurements to nonuniform load distribution and fiber alignment in the presence of variable tissue geometry and constitutive properties. Our first key finding was that regional shear wave speed tracks regional axial stress in a manner that is well de scribed by a tensioned beam model (Martin et al. 2018) . Our second key finding was that nonuniform fiber alignment acted to distr ibute axial stress across the tissue leading to more uniform regional stresses during nonuniform loading. Our last key finding was that the accuracy of axial stress predictions made using shear wave speed were generally independent of tissue geometry, nonu niform load, fiber material properties, and fiber alignment. Concerning our first objective, we found that regional shear wave speeds are an excellent predictor of the regional axial stress (Fig. 3). Based on a tensioned beam model, s hear wave speed was previously shown to be a good predictor of axial stress in an axially loaded fibrous tissue s (Martin et al. 2018) . We were motivated to consider nonuniform loading scena rios given that ligaments routinely undergo nonuniform deformations due to joint motion . Our model results suggest that shear wave speed may be a good predictor of regional tissue stress under such conditions. The regional w ave speed -stress relationships w ere relatively independent of fiber alignment and fiber linear region stiffness (§3.3.) , particularly at higher loads (J. L. Blank et al. 2021) . We found that regional wave speed -stress relationships were instead dependent on the extent of nonuniform loading in the tissue. More specifically, we observed diminished capacity of regional shear wave speeds to predict regional stress at low overall loads when the nonuniform load distribution was high (Fig. 3 a). This may indicate that near -unloaded or buckled regions of a nonuniformly -loaded tissue can introduce artifacts of higher shear wave propagation speeds on the more highly loaded regions of the tissue. However, the RMSE of predicted regional stress was small (RMSE = 0.5 5 MPa on average). The extent of the Blank et al. , 2023 15 nonuniform loading or stress gradient that can be applied while still accurately resolving regional stress using regional shear wave speed requires further study. Our second objective was to determine the sensitivity of regional axial stress and regional shear wave speed measurements to nonuniform load distribution and fiber alignment in the presence of variable tissue geometry and constitutive properties . First, we found that unaligned fibers distribute d axial stress more uniformly across the tissue during a nonuniform load (Fig. 5). We also found that the interaction between thickness and nonuniform load accounted for much of the variability in regional stress distribution ( R2 = 0.8 2). However, t his is likely due to thinner tissues having a decreased cross - sectional area, and thus a higher axial stress . At higher average axial stresses, regional stresses are more sensitive to nonuniform load ing because the fibers are more likely to be in the linear region of the fiber stress -strain response . More notably , we found that errors in predicting regional axial stress based upon regional shear wave speeds were generally independent of all geometric, boundary, and hyperelastic parameters included in this study . We did observe that more highly aligned fiber s resulted in greater stress distributions (Fig. 5 ) and thus greater average stress prediction errors (Fig. 3 b), which strongly suggests that shear wave speed is highly sensitive to measurement region in tendon -like tissues undergoing nonuniform loading. This should be a consideration in nonuniform load -bearing tendons , like the Achilles tendon (Franz et al. 2015) . In tissues that more effectively distribute stress, like ligament or joint capsule, regional shear wave speed variation may be less of a concern , but experimental follow -ups are needed to confirm . There are four limitations to consider when interpreting the res ults presented in this study. First, we implemented a finite element model using a range of parameters commonly found in soft biological tissues, but we did not consider other mechanical phenomena in these tissues that could Blank et al. , 2023 16 cause nonuniform loading, such as different hierarchical levels of load transmission like interfibrillar shearing and sliding (Szczesny and Elliott 2014; Rigozzi et al. 2013) or sliding between h igher -level tissue constituents such as sub -tendons (Pękala et al. 2017; Ekiert, Tomaszewski, and Mlyniec 2021) . However, we expect the results presented in this study to be representative of common structure - scale evaluations of ligaments and tendons . Second, there are numerous anatomical features of biological soft tissues that can give rise to nonuniform axial stress distribu tions. Here, we only considered linear stress distributions as caused by nonuniform loading, but in biological specimens a non-linear axial stress gradient can be the cause o f boundary conditions like the myotendinous junction (Sharafi et al. 2011) , bone insertions (Butler et al. 1984; Thomopoulos et al. 2003) , or geometric conditions such as changes in cross -section along the tissue length , such as in the Achilles tendon (Magnusson and Kjaer 2003) . These conditions during physiological loading can also give ris e to buckling in unloaded regions of the tissue, such as in the medial collateral ligament during simple knee flexion (John C. Gardiner, Weiss, an d Rosenberg 2001) . However, buckled regions of the tissue likely do not carry positive axial stress and are not useful for tensiometry stress measurements . Thus, we did not consider buckling caused by more extreme non uniform load di stributions (i.e., deep k nee flexion) in this study . Third , we measured shear wave speeds using a transient analysis and did not consider frequency -dependent phase speeds as measured using shear wave dispersion (Helfenstein -Didier et al. 2016; Brum et al. 2014) , which may a lter regional shear wave speeds observed in this model and in vivo. However, our finding that a simple transient analysis allowed for accurate predictions of regional stress suggests that this analysis can be applied despite the possible dispersive behavio r in biological tissues. Finally, shear wave speeds are measured through layers of subcutaneous tissue in vivo , which change s the effective density ( peff) of the tissue of interest, and may alter the interpretation of regional Blank et al. , 2023 17 shear wave speed at the super ficial measurement point (J. L. Blank and Thelen 2023; Sadeghi and Cortes 2020) . Whether regional loading can be ascertained using regional wave speeds in vivo is an area of ongoing study. In summary, based on our computational results, regional shear wave speeds are a strong predictor of regional stress in nonuniformly loaded fibro us soft tissues. We found that the regional shear wave speed -stress relationship was independent of factors that change between tissue types , and that this relationship was generally only sensitive to a combination of the de gree of nonuniform loading and total tissue load. As such, regional shear wave propagation is a promising measure to gauge regional tissue load directly and objectively. ACKNOWLEDGMENT We would like to thank the University of Wisconsin Center for High Throughput Computing for access to their computing resources. FUNDING This material is based upon work supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE -1747503. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. Blank et al. , 2023 18 NOMENCLATURE w Tissue width [mm] t Tissue thickness [mm] ΔF Nonuniform load distribution [N] F Axial load [N] C3 Toe region exponential stress coefficient [M pa] C4 Toe region uncrimping term C5 Linear region fiber modulus [M pa] kf Fiber alignment factor c Shear wave speed [m/s] peff Effective tissue density [kg/m3] μ Tissue shear modulus [M pa] σregional Regional axial stress [M pa] σaverage Average axial stress (F/A) [Mpa] εpred Regional stress prediction error [M pa] A Tissue cross -sectional area [mm2] Blank et al. , 2023 19