Over the last two decades, it has become clear that skeletal muscle exhibits complex nonuniform strain distributions during contraction (Azizi and Deslauriers, 2014, Pappas et al., 2002), which ultimately influences muscle force generation. This advancement in understanding emerged from measurements of tissue-level deformations in muscle using in vivo techniques, including cine MRI (Fiorentino et al., 2013, Hodgson et al., 2006, Pappas et al., 2002, Zhong et al., 2008), sonomicrometry (Gillis and Biewener, 2002, Soman et al., 2005), and ultrasound (Bennett et al., 2014, Leitner et al., 2020). While the overall conclusion that strain distributions are nonuniform is consistent across the literature, the precise patterns of deformation vary substantially across studies because local strain distributions are impacted by a wide variety of factors, such as muscle architecture, neuromuscular activation, muscle tissue properties, internal aponeurosis morphology, and more (Azizi and Deslauriers, 2014, Rehorn and Blemker, 2010). It is challenging to quantify how these factors influences muscle deformation through experiments alone, due to the difficulty in experimentally isolating the influence of each factor.

Computational modeling provides an excellent framework for capturing and understanding complex muscle deformation and contraction, isolating the influence of biomechanical factors, performing what-if experiments, and explaining cause-and-effect relationships. Indeed, finite element muscle modeling studies have shown that complex 3D muscle morphology (Blemker and Delp, 2005, Blemker et al., 2005, Fiorentino and Blemker, 2014, Fiorentino et al., 2014, Knaus et al., 2022, Rahemi et al., 2015, Röhrle and Pullan, 2007), aponeuroses (Chi et al., 2010, Fiorentino and Blemker, 2014, Fiorentino et al., 2014, Knaus et al., 2022, Rehorn and Blemker, 2010), and interactions with neighboring tissues and muscles (Maas et al., 2003) all contribute to strain distributions. Due to the nonlinear nature of the passive and active force–length characteristics of muscle, varying degrees of muscle activation and length change also lead to varied degrees of strain heterogeneity. For example, conditions of elevated activated lengthening, such as during sprinting, lead to greater levels of strain heterogeneity and thus higher localized strains, which provides an explanation for increased susceptibility of the hamstrings to strain injury (Fiorentino et al., 2014, Fiorentino and Blemker, 2014). While these simulations account for the nonlinear active and passive force–length characteristics of muscle fibers, they do not incorporate the force–velocity property of muscle (Hill, 1938). It stands to reason that the force–velocity property could have a balancing effect on force transmission between in-series fiber segments, as individual sarcomeres lengthen and shorten individually to compensate for their neighbors (Morgan, 1990, Morgan and Proske, 2006). Therefore, it is plausible that strain distributions predicted by finite-element simulations would be impacted by the force–velocity property.

Implementation of the force–velocity property in constitutive models of skeletal muscle tissue has been successful (Ehret et al., 2011, Marcucci et al., 2017, Ross et al., 2021); however, these formulations remain reliant on specialized software packages that limits dissemination to other researchers. Furthermore, it remains unclear how much the force–velocity relationship impacts strain distributions predicted by finite element models. Taking advantage of recent developments in free and open-source finite element modeling, such as FEBio (Maas et al., 2012) with its robust plugin system (Maas et al., 2018), can provide additional avenues to use and share specialized formulations of constitutive models.

The goals of this study were to: 1) develop, implement, and numerically validate a force–velocity constitutive model of muscle within FEBio that can be easily shared and implemented in both existing and new finite element models, and 2) use the force–velocity constitutive model to determine if, and under what conditions, the force–velocity property impacts tissue-level strain distributions. To explore this question, we created models with both simplified and complex geometries, simulated during both simplified and physiological conditions, including a model of the biceps femoris long head undergoing eccentric activation during sprinting.