1 2 3 4 5 6 7 Comparison of human gastrocnemius forces predicted by Hill-type muscle models and estimated from 8 ultrasound images 9 10 Taylor J.M. Dick1,*, Andrew A. Biewener2 and James M. Wakeling1 11 1Department of Biomedical Physiology and Kinesiology, Simon Fraser University, Burnaby, BC, Canada 12 2Concord Field Station, Harvard University, Bedford, MA, USA 13 * Corresponding author (taylor_dick@sfu.ca) 14 Telephone: 778-782-8445 15 Fax: 778-782-3040 16 17 18 19 20 21 22 23 24 Running Title: Testing gastrocnemii Hill-type models 25 Keywords: B-mode ultrasound, electromyography, motor unit recruitment, musculoskeletal simulation 1 26 Abstract 27 Hill-type models are ubiquitous in the field of biomechanics, providing estimates of a muscle’s 28 force as a function of its activation state and its assumed force-length and force-velocity properties. 29 However, despite their routine use, the accuracy with which Hill-type models predict the forces generated 30 by muscles during submaximal, dynamic tasks remains largely unknown. This study compared human 31 gastrocnemii forces predicted by Hill-type models to the forces estimated from ultrasound-based 32 measures of tendon length changes and stiffness during cycling, over a range of loads and cadences. We 33 tested both a traditional model, with one contractile element, and a differential model, with two 34 contractile elements that accounted for independent contributions of slow and fast muscle fibres. Both 35 models were driven by subject-specific, ultrasound-based measures of fascicle lengths, velocities, and 36 pennation angles and by activation patterns of slow and fast muscle fibres derived from surface 37 electromyographic recordings. The models predicted on average, 54 % the time-varying gastrocnemii 38 forces estimated from the ultrasound-based methods. However, differences between predicted and 39 estimated forces were smaller under low speed-high activation conditions, with models able to predict 40 nearly 80 % of the gastrocnemii force over a complete pedal cycle. Additionally, the predictions from the 41 Hill-type muscle models tested here showed that a similar pattern of force production could be achieved 42 for most conditions with and without accounting for the independent contributions of different muscle 43 fibre types. 44 45 46 47 48 49 50 2 51 Introduction 52 Hill-type muscle models are ubiquitous in biomechanical analyses of human movement and more 53 recently are being applied to animal studies to provide estimates of a muscle’s force as a function of its 54 activation state and assumed force-length and force-velocity properties. For example, Hill-type models 55 have been used, together with reconstructions of musculoskeletal geometry, to analyze the moment56 generating capacities of animals during terrestrial locomotion (O'Neill et al., 2013; Hutchinson et al., 2015) 57 and to estimate the locomotor capabilities of extinct species based on fossil records (Hutchinson and 58 Garcia, 2002). Muscle models have also been used in simulations of human locomotion to infer the 59 functions of individual muscles during walking (Arnold et al., 2013) and running (Hamner et al., 2010) and 60 to uncover important insights into the biomechanical factors that contribute to gait abnormalities (e.g. 61 Peterson et al., 2010; Steele et al., 2010). However, the accuracy of the muscle forces predicted by Hill62 type models, especially during submaximal, dynamic tasks, remains largely unknown. 63 The predictive accuracy of Hill-type models has been investigated in a small number of animal 64 studies in which tendon forces were directly measured. For example, Sandercock and Heckman (1997) 65 and Perreault et al. (2003) examined whether a Hill-type model could predict in situ forces generated by 66 the cat soleus when they imposed length changes corresponding to those measured during locomotion. 67 Their model reproduced soleus forces measured in vivo, using an implanted tendon buckle transducer, to 68 within 10 % of maximal tension during force rise (Sandercock and Heckman, 1997). However, at low motor 69 unit firing rates, differences in the predicted and measured forces were greater than 50 % (Perreault et 70 al, 2003). Our own efforts to predict in situ and in vivo forces generated by the gastrocnemii muscles of 71 goats (Wakeling et al., 2012; Lee et al., 2013) confirmed that Hill-type models are sensitive to assumptions 72 about activation state and the force-velocity properties of the model’s contractile element. For example, 73 we have shown that Hill-type models reproduce the time-varying in vivo forces with an average r2 of 0.40, 74 and errors in force greater than 15 % and 28 % of maximum isometric force for the medial (MG) and lateral 75 (LG) gastrocnemii of goats, respectively, when averaged across the gait cycle and different locomotor 76 speeds. However, errors in force predictions were reduced when models incorporated the independent 77 activations of slow and fast contractile elements, but only during the fastest locomotor speeds (Lee et al., 78 2013). Additionally, Millard and co-workers (2013) have demonstrated that Hill-type models are better 79 able to reproduce muscle forces under maximally activated conditions in comparison to submaximally 80 activated conditions. They compared the ability of three different Hill-type formulations to reproduce in 3 81 vivo forces measured in the cat soleus under maximally active (Krylow and Sandercock, 1997) and 82 submaximally active conditions (Perreault et al., 2003) and found that average errors were 12 % of 83 maximum tension for maximal contractions but increased to more than 17 % for submaximal contractions. 84 Together, these previous studies and others leave many unanswered questions in regards to the ability of 85 Hill-type models to reproduce time-varying muscle force. 86 To gain more insight into the strengths and weaknesses of Hill-type models, independent 87 estimates of time-varying muscle forces, for a wide range of movements, are needed. However, muscle 88 forces often cannot be directly measured. We have recently developed an ultrasound-based approach to 89 indirectly estimate the in vivo forces generated by the human triceps surae muscles during dynamic tasks 90 (Dick et al., 2016). In particular, the gastrocnemii-Achilles tendon (AT) complex is advantageous because 91 of its’ superficial location and short fascicles which allows us to use ultrasound to measure both tendon 92 and fascicle length changes during movement, and couple these measurements with subject-specific 93 estimates of tendon stiffnesses and slack lengths to estimate time-varying forces. This is the first study to 94 test Hill-type models against in vivo estimates of time-varying muscle forces from human subjects at a 95 range of mechanical conditions. 96 Electromyographic (EMG) recordings from the gastrocnemii of cyclists have shown that the 97 activation patterns of different muscle fibre types vary with cadence (Citterio & Agostoni 1984; Wakeling 98 et al., 2006; Wakeling and Horn, 2009). However, while muscles are composed of a mix of slow and fast 99 fibres with different physiological and biomechanical properties, Hill-type models are typically comprised 100 of a single contractile element with force-length and force-velocity properties that have been estimated 101 from in vitro data collected from single fibres under maximally-activated conditions. To represent whole 102 muscle, the fibre’s force is scaled based on the muscle’s level of activation, accounting for the muscle’s 103 physiological cross-sectional area, fibre length, pennation angle, and parallel elasticity. Here we test a 104 traditional Hill-type model, with one contractile element, compared with a differential model, having 105 independent slow and fast contractile elements to account for the contributions of slow and fast muscle 106 fibres. Both models were driven by subject-specific measures of fascicle lengths, velocities, and pennation 107 angles derived from ultrasound images and by activation patterns of slow and fast fibres derived from 108 surface EMG recordings and wavelet based time-frequency decomposition techniques (von Tscharner, 109 2000; Wakeling and Horn, 2009). We hypothesized that the two-element model would better reproduce 110 the ultrasound-based estimates of gastrocnemii force at the higher cadences, where preferential 4 111 recruitment of fast fibres has been reported (Citterio & Agostoni 1984; Wakeling et al., 2006). This is an 112 extension of our previous studies on the gastrocnemii muscles in goats, where we showed that a two113 element model, driven by the independent activation of slow and fast muscle fibres, provided better 114 predictions of time-varying muscle forces than traditional one-element models for some in situ (Wakeling 115 et al., 2012) and in vivo (Lee et al., 2013) conditions. 116 117 Materials and methods 118 We collected comprehensive sets of kinematic, kinetic, EMG, and ultrasound data from 16 119 competitive cyclists (8 male, 8 female; Suppl. Table 1). Subjects pedalled on a stationary cycle ergometer 120 (Indoor Trainer, SRM, Julich, Germany) at cadences between 60-140 r.p.m. and crank loads between 14121 44 N m. In this study, we analyzed data from 8 different combinations of cadence and load. We predicted 122 the time-varying forces generated by each subject’s LG and MG muscles during cycling using a one123 element Hill-type model and a two-element Hill-type model (Fig. 1). The models were driven by fascicle 124 lengths, velocities and pennation angles that we determined experimentally from ultrasound images, and 125 by muscle activations calculated from simultaneous recordings of surface EMG (Fig. 1). We compared each 126 subject’s predicted LG and MG forces to the forces estimated in vivo from ultrasound-based measures of 127 tendon length changes and stiffness (Fig. 1). Subjects gave informed consent and protocols were approved 128 by the Institutional Ethics Review Boards at Simon Fraser University and Harvard University. The 129 experimental protocol has been described and evaluated in a previous study (Dick et al., 2016) so a brief 130 overview is provided here. 131 Muscle models for estimating LG and MG forces 132 We predicted the time-varying forces generated by each subject’s LG and MG muscles during 133 pedalling using a traditional Hill-type model: 134 m = max[̂()̂a(̂f)̂a(̂) + ̂p(̂f)] cos (Eq. 1) 135 In previous efforts to predict gastrocnemius forces in goats, the intrinsic properties of the contractile 136 element have been formulated in several different ways (Wakeling et al., 2012). However, it was found 5 137 that the predicted forces were similar across the different formulations (Lee et al., 2013). In this study, we 138 assigned contractile element properties that were consistent with a ‘homogeneous one-element model’ 139 (Wakeling et al., 2012) that assumes the muscle consisted of fibres with homogenous properties. In this 140 model, the maximum intrinsic speed, often referred to as the maximum unloaded shortening velocity (̂0), 141 was the maximum intrinsic speed of the different fibre types weighted by their fractional cross-sectional 142 areas. The curvature of the force-velocity relationship (α) was an intermediate value between the slow 143 and fast fibres and was assumed to be the same all the fibres within the muscle (Zajac, 1989; Wakeling et 144 al., 2012; Lee et al., 2013) (Fig.1C). 145 We also predicted the subjects’ time-varying LG and MG forces using a ‘differential two-element 146 model’ (Wakeling et al., 2012; Lee et al., 2013) that incorporates independently activated slow and fast 147 contractile elements arranged in parallel (Fig.1C): 148 m = max [[1̂()fast̂a(̂f)̂a,fast(̂) + ̂()sloŵa(̂f)̂a,slow(̂)] + ̂p(̂f)] cos (Eq. 2) 149 In both models, total muscle force m is a function of the time-varying activation level ̂(), the normalized 150 active and passive force-length relationships ̂a(̂f) and ̂p(̂f), respectively, and the normalized force151 velocity relationship ̂a(̂) (Fig. 1C). max represents the subject-specific maximum isometric force of 152 either the LG or the MG: 153 max = ( m ) f,opt 0 (Eq. 3) 154 max was calculated based on the muscles’ scaled volumes m estimated using the regression equations 155 in Handsfield et al. (2014), and optimal fibre lengths f,opt estimated from a subject-specific scaled 156 musculoskeletal model (OpenSim v3.3; Delp et al., 2007; Arnold et al., 2010; Dick et al., 2016), 0 is 157 estimated specific muscle stress (225 kPa; Spector et al., 1980; Roy et al., 1982) (Suppl. Table 1). 1 is a 158 scaling factor that accounts for the lower power of the EMG signals that would be expected from action 159 potentials with higher spectral frequencies (Wakeling et al., 2012). c1 affects the balance between fast 160 and slow fibre contributions and it was given a value of 2 as described below. max scales the fibre (or 161 fascicle) force to whole muscle force and pennation allows the component of fascicle force that is 162 aligned with the tendon to be estimated. 6 163 The normalized active force-length curve (Fig. 1C; Otten, 1987) was given by: 164 ̂a(f) = −(̂f00.6.3−1)2.3 (Eq. 4) 165 The normalized passive force-length curve (Fig. 1C) was given by: 166 ̂p = 2.64̂f2 − 5.30̂f + 2.66 for ̂f > 1 (Eq. 5) 167 ̂p(̂f) = 0 for ̂f ≤ 1 (Eq. 6) 168 where ̂f is fascicle length normalized to fascicle slack length 0,f. 0,f was determined to be the length of 169 the fascicle at tendon slack length. This passive force-length curve is similar to that provided by Millard 170 and co-workers (2013) and is based on a combination of experimental data from chemically-skinned 171 human gastrocnemius fibres (Gollapudi and Lin, 2009) and rabbit whole muscles (Winters et al., 2011). 172 The normalized force-velocity curve (Fig. 1C) was given by: 173 ̂a() = 1+̂̂0 1−̂̂0α for ̂ ≤ 0 (Eq. 7) 174 ̂a() = 1.5 − 0.5 1−̂̂0 1+7̂.506α̂ for ̂ > 0 (Eq. 8) 175 where ̂ is fascicle velocity normalized to 0,f, which is equivalent to strain rate. α describes the curvature 176 of the force-velocity relationship and depends on fibre type (Otten, 1987). A literature survey of 59 species 177 from 88 papers concluded that the α values are 0.18 and 0.29 for slow and fast fibres, respectively 178 (Hodson-Tole and Wakeling, personal communication). ̂0 is the maximum intrinsic speed, and we used 179 values of 5 and 10 −1 for the slow and fast fibres, respectively (Wakeling et al., 2012; Lee et al., 2013). 5 180 and 10 −1 are thought to be within the appropriate range for human muscle (Faulkner et al., 1986; 181 Epstein and Herzog, 1998) and are similar to values reported for mouse, cat, and rat muscle measured at 182 physiological temperatures, which range from 4.8 to 7.3 −1 for slow fibres (Close, 1964; Spector et al., 183 1980; Askew and Marsh, 1998) and from 9.2 to 24.2 −1 for fast fibres (Close, 1964; Close, 1965; Luff, 7 184 1975; Spector et al., 1980). These values were also selected based on previous modelling results where 185 the sensitivity of predicted forces to different values of ̂0 were assessed (Lee et al., 2013). 186 Histological analysis has shown that the human gastrocnemii contain equal cross-sectional areas 187 of slow and fast muscle fibres (Johnson et al., 1973). Therefore, the one-element homogenous model 188 assumed that the active LG and MG have intrinsic properties that represent the average values for α and 189 ̂0 between slow and fast fibres ( =0.235; ̂0=7.5 −1) (Wakeling et al., 2012; Lee et al., 2013). 190 The muscle models (Eq. 1 and Eq. 2) were driven directly by the experimentally derived activation, 191 fascicle length and pennation angle, and the predicted muscle forces were subsequently compared to 192 estimates of the tendon force. The models did not involve the balancing of the fascicle force with tendon 193 force (series elastic element force), as is necessary if the whole muscle-tendon unit is incorporated in the 194 model. As such, our approach (which parallels earlier animal studies: Hodson-Tole and Wakeling, 2009; 195 Wakeling et al. 2012; Lee et al. 2013) circumvents the need to apply additional serial damping components 196 to prevent high-frequency oscillations in the solution (Günther et al., 2007; Millard et al., 2013; Haeufle 197 et al., 2014). 198 Experimental data for driving models 199 Subjects pedalled at eight different combinations of cadence and crank torque: 60 r.p.m. at 44 N 200 m, 80 r.p.m. at 14 N m, 80 r.p.m. at 26 N m, 80 r.p.m. at 35 N m, 80 r.p.m. at 44 N m, 100 r.p.m. at 26 N 201 m, 120 r.p.m. at 13 N m, and 140 r.p.m. at 13 N m, corresponding to crank powers of 275, 115, 220, 290, 202 370, 270, 165, and 200 W, respectively. Each trial lasted 15 seconds. Each block of eight experimental 203 conditions was repeated, in a random order, four times so that we could separately image both the LG 204 and MG muscle bellies and the LG and MG muscle-tendon junctions (MTJ). Subjects rested for five minutes 205 in between blocks of conditions. “Maximum effort” sprint trials (high power and cadence) were collected 206 at the beginning and end of each session in an effort to elicit maximum muscle activity, and we used these 207 data as a reference when normalizing the muscles’ EMG intensities (e.g., Rouffet and Hautier, 2008). We 208 used an optical motion capture system (Certus Optotrak, NDI, Waterloo, Canada) sampling at 100 Hz to 209 track the 3D locations of 32 LED markers placed bilaterally on the lower extremities, the pedals, and on 210 the ultrasound probe. Pedal reaction forces effective and ineffective (normal and radial) to the right and 211 left cranks were recorded at 2000 Hz (Powerforce, Radlabor, Freiburg, Germany). Previously, we have 8 212 verified that the plantarflexion moments estimated from ultrasound-based measures of tendon strain are 213 consistent with the net ankle moments determined from inverse dynamics (Dick et al., 2016). 214 Muscle fascicle lengths and pennation angles 215 A linear-array B-mode ultrasound probe (7 MHz, 60 mm field-of-view; Echoblaster, Telemed, 216 Vilnius, Lithuania) recording at 40 Hz was used to alternately image LG and MG fascicles from the right leg 217 during pedalling. The order in which muscles were imaged was randomized; at the end of each block of 218 conditions, the probe was re-positioned over the other muscle. The probe was secured in a custom-made 219 foam support and positioned to image muscle fascicles that were in the middle of the muscle belly (Suppl. 220 Fig. 1) where the fascicle architecture is homogeneous (Maganaris et al., 1998) and in plane with the 221 scanning image (Kawakami et al., 1993). An ultrasound gel pad (Parker Laboratories, NJ, USA) was placed 222 at the probe-skin interface to enhance image quality and allow muscle bulging. 223 Ultrasound images were manually digitized (ImageJ, NIH, Maryland, USA) to determine each 224 subject’s time-varying fascicle lengths and pennation angles for the eight pedalling conditions. Eight points 225 in each ultrasound image were digitized—three points on each of the superficial and deep aponeuroses 226 and two points on a representative muscle fascicle (Suppl. Fig. 2). Pennation angle was calculated as the 227 mean of the angles made by the fascicle with the deep and superficial aponeuroses. Fascicle lengths and 228 pennation angles from 10 complete crank revolutions were fit using a least-squares minimization of a 229 Fourier series [1 + 1Sin(∅1 + ) + 2Sin(∅2 + 2)], where 1, 1, 2, ∅1, ∅2 are fitted constants and 230 is the crank angle. These fascicle lengths and pennation angles were used as inputs to the muscle 231 models for each muscle, each pedalling condition, and each subject. 232 Muscle activation patterns 233 Bipolar Ag/AgCl surface electrodes (10 mm diameter and 21 mm spacing; Norotrode; Myotronics, 234 Kent, WA, USA) were used to record muscle excitations from the LG, MG, soleus (SOL), and seven other 235 muscles on the contralateral (left) limb to which the ultrasound transducer was placed. EMG signals were 236 preamplified (gain 1000), band-pass filtered (bandwidth 10–500 Hz; Biovision, Wehrheim, Germany), and 237 sampled at 2000 Hz as described elsewhere (e.g., Wakeling and Horn, 2009; Blake and Wakeling, 2014). 238 Electrodes were placed in the mid-region of the muscle belly after the skin had been shaved and cleaned 9 239 with isopropyl alcohol solution. We secured electrodes with stretchable adhesive bandages and tubular 240 net bandages to reduce movement artefacts during pedalling. 241 Surface EMG signal has been shown to contain distinct frequency characteristics that allow for 242 characterization of recruitment between different types of motor units (Wakeling and Rozitis, 2004; Reaz 243 et al., 2006; Wakeling, 2009). We characterized EMG signals from the LG and MG in two different ways: 244 first by total EMG intensity, representing all active motor units within the muscle, and second by the EMG 245 intensity at the low and high frequency bands, to represent the contributions of slow and fast motor units, 246 respectively. EMG signals for the muscles were quantified by their intensities during each crank 247 revolution. The total intensity is a close approximation to the power of the signal and was calculated 248 across a 10–450 Hz frequency band using an EMG-specific wavelet analysis (von Tscharner, 2000). 249 Optimized wavelets were derived, for each muscle, using principal component analysis to identify the 250 major features of the intensity spectra from the 16 subjects (Wakeling and Rozitis, 2004; Von Tscharner 251 and Goepfert, 2006; Hodson-Tole and Wakeling, 2007; Wakeling and Horn, 2009; Lee et al., 2011) (Fig. 2). 252 A least-squares minimization of a wavelet function () (Eq. 9) was used to construct two optimized 253 wavelets to describe the low and high frequency spectra that correspond to the slow and fast motor units: 254 () = ( )c ∙ ∙ (−c+1)∙c ∙ c (Eq. 9) 255 where is a scaling factor that defines the width and shape of the wavelet and c is the central frequency 256 of the wavelet (von Tscharner, 2000) (Table 1). 257 To compare predicted forces derived from the total intensities (one-element model) to the forces 258 derived from the optimized wavelets (two-element model), the time resolutions from the two approaches 259 were defined to be the same; to do this, we set the Gauss filters at the end of the wavelet analysis to a 50 260 ms time resolution for all wavelets, which is within the range of physiological time to peak twitch for 261 skeletal muscle (10 ms and 100 ms, Levin et al., 1999; Spaegele et al., 1999). We used the square-root of 262 the EMG intensity as a measure of muscle excitation , as muscle force is linearly related to the EMG 263 amplitude and not its power (Milner-Brown and Stein, 1975). A lower EMG intensity would be expected 264 from action potentials with higher spectral frequencies (Wakeling et al., 2012). Here we selected a value 265 of c1 = 2 (in Eq. 2) to scale the excitations (based on the ratio of their central frequencies: Table 1) to 266 account for this effect. 10 267 We used a first-order differential equation (transfer function) to estimate the active state of the 268 muscle from the muscle excitation (Zajac, 1989) that is commonly used for modelling human muscle. This 269 choice was partly motivated by a parallel project that will compare the muscle models from this study to 270 muscle force predictions from the OpenSim musculoskeletal simulation environment (Delp et al. 2007) 271 that also uses this transfer function. However, it should be noted that other transfer functions are capable 272 of estimating better activation dynamics (Hatze, 1977; Lee et al. 2011; Rockenfeller et al. 2015) but require 273 more model parameters that are currently unknown for muscle in man. 274 () + [a1ct (act + [1 − act]())] () = (1) act () (Eq. 10) 275 The activation time constant act, specific to slow and fast fibres, and the ratio act of the activation to 276 deactivation time constants were derived using tetanic contraction data in the literature (Table 2). 277 Histological analysis has shown that the human gastrocnemii contain equal cross-sectional areas of slow 278 and fast muscle fibres (Johnson et al., 1973). Therefore, act for the activation transfer function for the 279 “one-element” model was calculated as the average between slow and fast fibres. The ratio of activation 280 to deactivation, act, was independent of fibre-type since both activation and deactivation rates are 281 associated with the motor-unit type (Burke et al., 1973). We used the transfer function to derive the active 282 state for the whole muscle (from the total EMG intensity), as well as for slow and fast motor units. In order 283 to make comparisons between the one- and two-element models, the activation level of the summed 284 slow and fast motor units was scaled in amplitude to equal the total activation during the maximum-effort 285 reference pedalling trials. Activation traces from 10 cycles were combined to compute a mean LG and MG 286 total, slow, and fast muscle activation trace for each pedalling condition. Activation levels for each muscle 287 are presented as normalized values ̂() calculated relative to the maximum activation detected during 288 the reference maximum-effort pedalling trials. 289 Experimental data for evaluating models 290 We estimated in vivo muscle forces from subject-specific tendon length changes and stiffness to 291 evaluate the predictions of our Hill-type models. A linear-array B-mode ultrasound probe was secured 292 over the LG or MG distal muscle-tendon junction. A rigid-body cluster of LED markers (Prager et al., 1998) 293 firmly attached to the ultrasound probe, allowed the 3D coordinates of the LG or MG muscle-tendon 294 junction to be tracked (ImageJ, NIH, Maryland, USA). We calculated the length changes of each tendon 11 295 AT−LG or AT−MG during cycling as the distance from the AT insertion on the calcaneus (identified using 296 an LED marker) to the 3D position of the LG or the MG muscle-tendon junction. We estimated each 297 muscle’s contribution to tendon force, LG and MG (equation shown for MG), during cycling based on the 298 measured length changes AT−LG or AT−MG, the subject-specific values of toe region stiffness SEE−LG,T, 299 SEE−MG,T and linear region stiffness SEE−LG, SEE−MG for either the LG or MG (Suppl. Table 1), and the 300 tendon slack length 0,AT−LG or 0,AT−MG: 301 0, AT−MG ≤ 0,AT−MG MG = {SEE,−MG,T(AT−MG − 0,AT−MG), 0,AT−MG < AT−MG ≤ (0,AT−MG)1.0103 SEE−MG(AT−MG − (0,AT−MG)1.0103) + max,T−MG, (0,AT−MG)1.0103 < AT−MG (Eq. 11) 302 Tendon stiffnesses were estimated from ramped isometric tests (Morrison et al., 2015; Dick et al., 2016), 303 where we quantified the total stiffness of the AT and also determined the proportion of the total AT 304 stiffness contributed by the LG (SEE−LG) and MG (SEE−MG) separately. These muscle-specific stiffnesses 305 were calculated based on ratio of the maximum force-generating capacity max of either the LG or MG to 306 the combined triceps surae maximum force (Suppl. Table 1). Each tendon force-length curve had toe 307 region and a linear region. Linear stiffnesses, SEE−LG and SEE−MG, were subject-specific but due to the 308 difficulty in measuring the toe region properties of the force-length curve in all subjects during the 309 isometric protocol (Dick et al., 2016), the same toe region stiffness, SEE−LG,T and SEE−MG,T, was used 310 for the LG and the MG in all subjects. max,T−MG corresponds to the estimated force at the end of the toe 311 region. Tendon slack lengths for the LG and MG portion of the AT (0,AT−LG and 0,AT−MG) were measured 312 at 310° of the crank cycle, averaged over all cycles; this choice was motivated by published tendon buckle 313 data (Gregor et al., 1987) showing force rise to occur after 310°. The time-varying tendon force from 10 314 cycles was used to compute a mean force trace for each condition. 315 Comparisons of predicted and estimated forces 316 We used the one- and two-element models to predict time-varying LG and MG forces for the 16 317 subjects and 8 pedalling conditions. Additionally, because the forces predicted by Hill-type models are 318 sensitive to the assumed maximum intrinsic speed ν0 and curvature α of the force-velocity relationship 319 (Shue et al., 1995; Perreault et al., 2003; Wakeling et al., 2012; Lee et al., 2013), we performed a sensitivity 320 analysis and ran the models varying the values of α (0.155-0.315) and 0 (4-10 −1) to determine the 321 extent to which the modelled forces were sensitive to these force-velocity properties. We compared the 12 322 predicted forces to the forces estimated from ultrasound-based measures of tendon length changes and 323 stiffness, and characterized differences across the crank cycle using two measures: the coefficient of 324 determination r2 and the root mean square error (RMSE). A general linear model ANOVA was conducted 325 to determine if differences in the r2 and RMSE (dependent variable) existed between models, muscles, 326 cycling conditions, subject (random), and values of α and 0. We used a Chi-square test at each pedalling 327 condition to determine whether the two-element model reproduced the estimated forces in the LG and 328 MG better than the one-element model (higher r2 and lower RMSE) over the entire crank cycle, more 329 often than the one-element model. Additionally, an ANCOVA identified changes in the mean muscle 330 activity for total, fast, and slow motor units with cadence, crank load, and fascicle length ̂f as covariates 331 and subject as a random factor. ̂f was included as a covariate because muscle strain has been shown to 332 affect the frequency content of the EMG signal (Doud and Walsh, 1995). Differences were considered 333 significant at the p<0.05 level. Data are reported as mean ± SE. 334 335 Results 336 The Hill-type models tested here, driven by EMG-derived activations and ultrasound-based 337 measures of fascicle lengths, velocities and pennation angles, were able to capture the general features 338 of the ultrasound-based estimates of force (Fig. 3). On average the models performed with an r2 of 0.54 339 and an RMSE of 13.52 % Fmax across muscles, conditions, and models (Table 3; Suppl. Fig. 3). However, 340 models performed even better at low cadence-high load conditions with r2 values as high as 0.62 and 0.79 341 and errors as low as 8.98 % Fmax and 9.86 % Fmax, for the LG and MG, respectively. A comparison of the 342 modelled forces between the pedalling conditions, between muscles, and between the one- and two343 element models are shown in Table 3 and Suppl. Fig. 3. 344 In both the LG and MG, there were minor differences between the muscle forces predicted by the 345 two-element model when compared to the ultrasound-based estimates of tendon force and those 346 predicted by the traditional one-element model. In particular, when averaged across subjects and 347 conditions, the two-element model predicted forces with r2 values of 0.47±0.09 and 0.63±0.1 for the LG 348 and MG, respectively, compared to 0.46±0.1 and 0.61±0.12 for the one-element model LG and MG, 349 respectively. RMSE errors were similar between the two models. The two-element model did perform 13 350 significantly better than the one-element model at high cadences (100 r.p.m. to 140 r.p.m.) (Suppl. Table 351 2), although the differences were small (Table 3; Suppl. Fig. 3). 352 Consistent with previous literature (Wakeling et al., 2006; Blake and Wakeling, 2014) the 353 conditions tested in our study elicited a range of activation patterns between slow and fast muscle fibres, 354 providing rationale to test the two-element model (Fig. 4). When averaged over the whole pedal cycle, 355 the total activation increased with both crank load and cadence in both the LG and MG across all subjects 356 (p<0.05). Activation of the slow motor units decreased with cadence (LG: p=0.02, MG: p=0.045) but 357 increased with crank load (LG and MG: p<0.01). Fast motor unit activation increased with both cadence 358 and load in the LG and MG (p<0.01). Therefore, increased cadence resulted in a reduction of slow motor 359 unit activity coupled with an increase in the fast motor unit activity. 360 The Hill-type models tested in this study reproduced the estimated forces for the MG more 361 accurately than for the LG. Both the one- and two-element models predicted forces with a higher r2 over 362 the crank cycle for the MG as compared to the LG across all pedalling conditions (p<0.05) (Table 3; Suppl. 363 Fig. 3). The average coefficient of determination r2 across subjects, conditions and models decreased from 364 0.62 (range: 0.45-0.79) for the MG to 0.47 (range: 0.29-0.62) for the LG. However, the average RMSE over 365 the crank cycle was similar between muscles: 13.55 % max (range: 9.86-18.25) for the MG and 13.51 % 366 max (range: 9.98-21.91) for the LG. 367 Time-varying patterns of the predicted and estimated forces, as assessed by the coefficient of 368 determination r2 and RMSE, showed greater differences at high cadences in comparison with low 369 cadences (Table 3; Suppl. Fig. 3); this was a common feature for both the one-element and two-element 370 models. In particular, when averaged across the 2 different models, the RMSE over the crank cycle was 371 11.1±1.2 % Fmax at cadences 60-80 r.p.m. and increased to 17.6±2.1 % Fmax at cadences 100-140 r.p.m. 372 when averaged across muscles, subjects, and conditions. 373 The forces predicted by the Hill-type models tested here were sensitive to the force-velocity 374 curvature α used (Fig. 5). In both the LG and MG, the models predicted force with a higher r2 and lower 375 RMSE (data not shown) when using α =0.235 or α =0.315 for the one-element model, as compared to α 376 =0.155 when averaged across all pedalling conditions (Fig. 5A). The most notable differences in model 377 performance with α were evident at cadences of  100 r.p.m., where models predicted force with a higher 14 378 r2 and lower RMSE when using greater α values (Fig. 5A). In contrast, the predicted forces were not 379 sensitive to the maximum intrinsic speed 0 implemented within the models. These trends were similar 380 for the one-element and two-element models tested (data shown for the one-element model; Fig. 5). 381 382 Discussion 383 Our results reveal that Hill-type muscle models, driven by EMG-derived activations and 384 ultrasound-based measurements of fascicle lengths, velocities and pennation angles, predicted on 385 average, 54 % the gastrocnemii forces generated by human subjects during cycling, as estimated via 386 ultrasound-based estimates of tendon length changes and stiffness. In general the models performed best 387 when the fascicles shortened at lower velocities and the activation levels were high, with models at these 388 conditions able to predict nearly 80 % of the medial gastrocnemii force over a complete pedal cycle. This 389 work provides the first comparison of forces predicted by Hill-type models –driven with in vivo human 390 experimental data under submaximal dynamic tasks –to ultrasound-based estimates of muscle-tendon 391 force. The two-element model with independent slow and fast contractile elements showed no significant 392 improvement for predicting time-varying and maximum muscle forces, when compared with a traditional 393 one-element model for the slower conditions tested, however, small but significant improvements were 394 noted at the fastest cadences. This suggests that other factors play an important role in determining time395 varying patterns and magnitudes of muscle force in addition to the contractile properties of the 396 constituent motor units. Identifying the underlying sources of error in force predictions from Hill-type 397 models remains challenging, particularly when having to rely on non-invasive estimates of muscle-tendon 398 force. However, the approach described here provides an important step and an experimental framework 399 for investigating how Hill-type models may be improved. 400 Differences between muscles and pedalling conditions 401 The performance of the Hill-type models depended on which muscle was being modelled. The 402 human LG and MG vary in architecture (Maganaris et al., 1998), activation patterns (Wakeling et al., 2011), 403 and motor-unit twitch profiles (Vandervoort and McComas, 1983). Specifically, the MG has shorter, more 404 pennate fascicles than the LG (Wakeling et al., 2006; Ward et al., 2009), while the LG has a more complex 405 architecture than the MG, with multiple heads that have different fascicle arrangements (Wolf et al., 15 406 1998). This complexity may have diminished our ability to accurately estimate representative fascicle 407 lengths, velocities, and pennation angles in the LG using ultrasound—and suggests that models are likely 408 more accurate for muscles with homogeneous architecture. Consequently, these differences may help to 409 explain the variation in model performance shown here, particularly the finding that models reproduced 410 the estimated forces more closely across the pedal cycle for the MG than the LG, which is consistent with 411 the results from the gastrocnemii of goats (Lee et al., 2013). 412 The muscle models were compared to estimates of tendon force that were derived from 413 ultrasound measures of tendon length and subject-specific tendon stiffness (methods described in Dick et 414 al., 2016). The tendon stiffnesses for the LG and MG portions of the AT were previously determined on an 415 isometric frame in which the knee was held at 130 ˚ and the subjects performed maximal effort 416 plantarflexions (Morrison et al., 2015; Dick et al., 2016). However it is possible that at different knee 417 angles, where the relative lengths of the LG, MG, and SOL are affected differently due to the biarticular 418 nature of LG and MG versus the uniarticular nature of SOL, we may get a different value of stiffness and 419 thus a different value of force. During cycling, in the middle of the pedal upstroke (270 ˚) the AT is 420 unloaded and muscle activation is zero; consequently, both estimated tendon force and predicted muscle 421 force will be zero. Towards the top of the pedal cycle, the AT starts to become loaded; however, as the 422 knee is still relatively flexed the contribution of the gastrocnemii to the combined AT force is likely less 423 relative to SOL than occurred during our isometric measures of force and tendon strain. This may lead to 424 an over-estimation of the tendon force during this phase, although the forces predicted by the muscle 425 models are independent of tendon stretch and are thus unaffected by this. However, when the knee angle 426 reaches 130 ̊ at a crank angle of 90 ̊, this matches the isometric test conditions and coincides with the 427 period when the muscle forces are the greatest. Additionally, the ankle plantarflexors (LG, MG, and SOL) 428 have a similar activity during the low cadence-high load conditions (Dick et al., 2016), which most closely 429 matches the isometric conditions for which tendon stiffnesses were determined. There was also a greater 430 resolution of ultrasound frames per pedal cycle at the lowest cadences due to the finite sampling rate of 431 the ultrasound scanner. Thus, we would expect the best match between the modelled muscle forces, and 432 the estimated tendon stiffnesses to occur for the low cadence and high load conditions. Indeed, the r2 433 reached 0.79 with a corresponding RMSE of 10 % Fmax for these conditions in the MG (Table 3), and this is 434 likely the limit to the accuracy that can be achieved using these non-invasive experimental techniques. 16 435 The accuracy of the models tested here is partly limited by the measures of fascicle length, 436 pennation and tendon length as well as estimates of fascicle and tendon slack lengths. We have previously 437 shown that small increases in tendon slack length delayed the onset of AT force during cycling and 438 decreased the magnitude of the predicted forces (Dick et al., 2016). Additionally, fascicle slack length is 439 set at the time when the tendon is at its slack length. It is possible that the accuracy of these models could 440 be further improved if the fascicle and tendon slack lengths were independently determined by 441 optimization (e.g., Gerus et al., 2012; Gerus et al., 2015), rather than being experimentally guided, 442 however this was not the purpose of this study. Additionally, the overall accuracy of the muscle model is 443 scaled by Fmax (equation 3), that was determined from both an optimal fibre-length and an estimated 444 muscle volume that were derived from regression equations from other studies (Delp et al., 2007; 445 Handsfield et al. 2014). 446 Comparison to Hill-type models in the literature 447 Previous studies have assessed the accuracy of Hill-type models against direct measures of force 448 from tendon buckle recordings in animals (Sandercock and Heckman, 1997; Perreault et al., 2003; 449 Wakeling et al., 2012; Lee et al., 2013; Millard et al., 2013; Kim et al., 2015) or against measures of heat 450 and work (Umberger et al., 2003; Lichtwark and Wilson, 2005). In particular, modelling in situ forces 451 yielded higher r2 values (Sandercock and Heckman, 1997; Perreault et al., 2003; Wakeling et al., 2012) 452 than in vivo forces (Lee et al., 2013), which likely relates to the more controlled contractions studied during 453 in situ experiments in comparison to in vivo experiments. The errors achieved in this study for the slower 454 cadences (RMSE of 10 % Fmax with r2 of 0.79) match the best in vivo models (RMSE of 10 %; Lee et al. 2013) 455 that have been achieved in animal studies where the data were collected from more invasive techniques 456 (fine-wire EMG for activation, sonomicrometry for fascicle lengths, and tendon buckles for force). 457 Differences between the forces predicted by the one- and two-element models in this study were 458 smaller than the differences predicted by similar one- and two-element models in the goat gastrocnemii 459 (Lee et al., 2013). However, both the goat study and this current study did find the only significant 460 improvements in the two-element model to occur at the fastest muscle strain-rates. Surprisingly, the 461 greater benefits from the two-element model occurred in the goat study that did not show such fast 462 strains rates as seen here. However, the differences in model performance may reflect the additional 463 challenges of extracting EMG information from surface EMG signals in man, or from the effect whereby 17 464 contractile differences between muscle fibre-types may be mitigated in larger muscles that are 465 submaximally activated (Holt et al. 2014; Ross et al. 2016). 466 Intrinsic properties and the two-element model 467 Previous studies have shown that the errors in forces predicted by Hill-type models are often 468 related to the assumed maximum intrinsic speed 0 and curvature α of the force-velocity relationship 469 (Shue et al., 1995; Perreault et al., 2003; Wakeling et al., 2012; Lee et al., 2013). In contrast to our previous 470 study (Lee et al., 2013), predicted forces of the models here were not sensitive to maximum intrinsic speed 471 0 (Fig. 5). However, predicted forces were sensitive, in a cadence-specific manner, to the curvature α of 472 the force-velocity relationship. Most notably, models reproduced the estimated forces more closely 473 (higher r2 and lower RMSE over the pedal cycle) using low α values at low cadences, but models performed 474 better using higher α values at high cadences (Fig. 5A). For a given shortening velocity, a higher α would 475 result in a larger force. As mentioned previously, slow fibres are characterized by a higher curvature (lower 476 α) whereas fast fibres are characterized by a lower curvature (higher α) (Otten, 1987; Umberger et al., 477 2003). These results suggest that the forces predicted by traditional Hill-type models may be improved by 478 incorporating a force-velocity curvature α that varies appropriately with type of muscle fibres recruited 479 for the specific mechanical demands of the task— a larger α for high speed tasks where faster fibres are 480 recruited, and a smaller α for low speed tasks where slower fibres are recruited. Umberger and colleagues 481 (2003) have previously alluded to this, and have shown that Hill-type models with force-velocity 482 parameters, to include 0 and α, that are matched to the proportions of slow and fast muscle fibres, 483 generate reasonable predictions of whole muscle energetics during isolated muscle contractions, single 484 joint motion, and whole body movement. 485 Consistent with previous cycling studies (Wakeling et al., 2006; Wakeling and Horn, 2009; Blake 486 and Wakeling, 2014), higher cadences elicited differential recruitment patterns that displayed an increase 487 in the activation of the faster motor units coupled with a decrease in the activation of the slower motor 488 units (Fig. 4). However, the differences between predicted forces for the one- and two-element models 489 at even these highest cadences were minor (average 3 %). This could be attributed to the fact that (i) the 490 one-element model was not sensitive to ̂0, (ii) the differences in recruitment across the range of loads 491 and cadences tested here may not have been large enough to elicit substantial differences in model 492 performance, or (iii) the forces from human-sized muscle at submaximal activations may be less sensitive 18 493 to differences in fibre-type recruitment than previously expected. To date, our understanding of whole 494 muscle behaviour has largely come from studies in which muscles are maximally activated, and from 495 single-fibre data or data from smaller muscles (< 10 g). However, recent evidence suggests that muscle 496 behaviour also depends on the muscle’s size (Ross and Wakeling, 2016) and level of activation (Rack and 497 Westbury, 1969; Josephson and Edman 1988; Rassier et al., 1999; Holt and Azizi, 2014; Holt et al., 2014; 498 Holt and Azizi, 2016; Ross and Wakeling, 2016). Understanding the role of different fibre-types within 499 large mixed muscles during voluntary submaximal contractions is clearly an avenue that warrants further 500 investigation. 501 Conclusions 502 In this study, we compared gastrocnemii forces predicted by Hill-type models to the forces 503 estimated from ultrasound-based measures of tendon length change and stiffness. Similar forces were 504 predicted for the traditional one-element model and the two-element model that accounted for the 505 independent contributions of different muscle fibre types across most of the speeds tested, and the 506 models performed best when fascicle velocities were low and muscle activation was high. What additional 507 factors should we consider in order to yield more accurate muscle models? Future work should aim to 508 additionally consider the mechanical effects of a muscles’ physical properties (e.g., mass, viscosity), as 509 well as the potential contributions of elastic tissue (e.g., aponeurosis) and surrounding structures on 510 muscle contractile function. This will likely be critical for improving the reliability of muscle models as well 511 as for understanding the mechanisms underlying whole-muscle function during dynamic sub-maximal 512 contractions—which we currently know little about. Benchmark experimental data sets, like the data 513 we’ve analyzed here, are necessary and useful for testing and improving models. 514 515 Acknowledgements 516 The authors thank Dr Allison Arnold for insightful discussions. 517 518 Competing Interests 519 No competing interests declared. 520 19 521 Author contributions 522 T.D, A.B., and J.W. designed the experiments. T.D. and J.W. conducted the experiments. T.D. analyzed 523 the data and wrote the first draft of the manuscript. J.W. and A.B. provided guidance throughout the 524 study, assisted with writing, and approved the manuscript. 525 526 Funding 527 We gratefully acknowledge funding from NIH Grant R01 2R01AR055648 and an NSERC Graduate Research 528 Fellowship to T.J.M. Dick. 529 530 References 531 1. Arnold, E.M., Ward, S.R., Lieber, R.L. & Delp, S. L. (2010). A model of the lower limb for analysis of 532 human movement. Ann. Biomed. Eng., 38(2), 269-279. 533 2. Arnold, E. M., Hamner, S. R., Seth, A., Millard, M., & Delp, S. L. (2013). How muscle fiber lengths and 534 velocities affect muscle force generation as humans walk and run at different speeds. J. Exp. Biol., 535 216(11), 2150-2160. 536 3. Askew, G.N. & Marsh, R.L. (1998). Optimal shortening velocity (V/Vmax) of skeletal muscle during 537 cyclical contractions: length-force effects and velocity-dependent activation and deactivation. J. 538 Exp. Biol., 201: 1527-1540. 539 4. Blake, O. M., & Wakeling, J. M. (2014). Early deactivation of slower muscle fibres at high movement 540 frequencies. J. Exp. Biol., 217(19), 3528-3534. 541 5. Burke, R. E., Levine, D. N., Tsairis, P. & Zajac, F. E. (1973). Physiological types and histochemical 542 profiles in motor units of the cat gastrocnemius. J. Physiol., 234, 723-748. 543 6. Citterio, G. & Agostoni, E. (1984). Selective activation of quadriceps muscle fibers according to 544 bicycling rate. J. Appl. Physiol. 57, 371-379. 545 7. Close, R.I. (1964). Dynamic properties of fast and slow skeletal muscles of the rat during 546 development. J. Physiol. 173: 74-95. 14. 547 8. Close, R.I. (1965). Force: velocity properties of mouse muscles. Nature 206: 718-719. 20 548 9. Delp, S. L., Anderson, F. C., Arnold, A. S., Loan, P., Habib, A., John, C. T., ... & Thelen, D. G. (2007). 549 OpenSim: open-source software to create and analyze dynamic simulations of movement. IEEE 550 Trans. Biomed. Eng., 54(11), 1940-1950. 551 10. Dick, T.J.M., Arnold, A.S., Wakeling, J.M. (2016). Quantifying Achilles tendon force in vivo from 552 ultrasound images. J. Biomech. 49(14), 3200-320. 553 11. Doud, J. R. & Walsh, J. M. (1995). Muscle fatigue and muscle length interaction: effect on the EMG 554 frequency components. Electromyogr. Clin. Neurophysiol., 35, 331-339. 555 12. Epstein, M. & Herzog, W. (1998). Theoretical Models of Skeletal Muscles: Biological and 556 Mathematical Considerations (New York: John Wiley and Sons). 557 13. Faulkner, J. A., Claflin, D. R., & McCully, K. K. (1986). Power output of fast and slow fibers from 558 human skeletal muscles. Human Muscle Power, 81-94. 559 14. Gerus, P., Rao, G., & Berton, E. (2012). Subject-specific tendon-aponeurosis definition in Hill-type 560 model predicts higher muscle forces in dynamic tasks. PloS one, 7(8), e44406. 561 15. Gerus, P., Rao, G., & Berton, E. (2015). Ultrasound-based subject-specific parameters improve 562 fascicle behaviour estimation in Hill-type muscle model. Comput Methods Biomech Biomed Engin., 563 18(2), 116-123. 564 16. Gollapudi, S. K. & Lin, D. C. (2009). Experimental determination of sarcomere force–length 565 relationship in type-I human skeletal muscle fibers. J. Biomech., 42, 2011-2016. 566 17. Gregor, R. J., Komi, P. V., & Järvinen, M. (1987). Achilles tendon forces during cycling. Int. J. Sports. 567 Med., 8, 9-14. 568 18. Günther, M., Schmitt, S., & Wank, V. (2007). High-frequency oscillations as a consequence of 569 neglected serial damping in Hill-type muscle models. Biol. Cybern., 97(1), 63-79. 570 19. Haeufle, D.F.B., Günther M., Bayer, A., Schmitt, S. (2014). Hill-Type Muscle Model with Serial 571 Damping and Eccentric Force-Velocity Relation. J. Biomech., 47(6), 1531–1536. 572 20. Hamner, S.R., Seth, A., Delp, S.L. (2010). Muscle contributions to propulsion and support during 573 running. J. Biomech., 43(14):2709-16. 574 21. Handsfield, G. G., Meyer, C. H., Hart, J. M., Abel, M. F., & Blemker, S. S. (2014). Relationships of 35 575 lower limb muscles to height and body mass quantified using MRI. J. Biomech., 47(3), 631-638. 21 576 22. Hatze, H. (1977). A myocybernetic control model of skeletal muscle. Biol. Cybern., 25(2), 103-119. 577 23. Hodson-Tole, E.F., Wakeling, J.M. (2007). Variations in motor unit recruitment patterns occur within 578 and between muscles in the running rat (Rattus norvegicus). J. Exp. Biol., 210:2333–45. 579 24. Hodson-Tole, E. F., & Wakeling, J. M. (2009). Motor unit recruitment for dynamic tasks: current 580 understanding and future directions. J. Comp. Physiol. B, 179(1), 57-66. 581 25. Holt, N. C., & Azizi, E. (2014). What drives activation-dependent shifts in the force–length curve?. 582 Biol. Lett., 10(9), 20140651. 583 26. Holt, N. C., Wakeling, J.M., & Biewener, A.A. (2014). The effect of fast and slow motor unit activation 584 on whole-muscle mechanical performance: the size principle may not pose a mechanical paradox. 585 Proc. R. Soc. B, 281(1783), 20140002. 586 27. Holt, N. C., & Azizi, E. (2016). The effect of activation level on muscle function during locomotion: 587 are optimal lengths and velocities always used?. Proc. R. Soc. B., (Vol. 283, No. 1823, p. 20152832). 588 28. Hutchinson, J. R., & Garcia, M. (2002). Tyrannosaurus was not a fast runner. Nature, 415(6875), 1018589 1021. 590 29. Hutchinson, J. R., Rankin, J. W., Rubenson, J., Rosenbluth, K. H., Siston, R. A., & Delp, S. L. (2015). 591 Musculoskeletal modelling of an ostrich (Struthio camelus) pelvic limb: influence of limb orientation 592 on muscular capacity during locomotion. PeerJ, 3, e1001. 593 30. Johnson, M.A., Polgar, J., Weightman, D., Appleton, D. (1973). Data on the distribution of fibre types 594 in thirty-six human muscles. An autopsy study. J. Neurolog. Sci., 18:111– 129. 595 31. Josephson, R. K., & Edman, K. A. P. (1988). The consequences of fibre heterogeneity on the force‐ 596 velocity relation of skeletal muscle. Acta Physiol. Scand., 132(3), 341-352. 597 32. Kawakami, Y., Abe, T., Fukunaga, T. (1993). Muscle-fiber pennation angles are greater in 598 hypertrophied than in normal muscles. J. Appl. Physiol., 74:2740–2744. 599 33. Kernell, D., Eerbeek, O., & Verhey, B. A. (1983). Relation between isometric force and stimulus rate 600 in cat's hindlimb motor units of different twitch contraction time. Exp. Brain Res., 50(2-3), 220-227. 601 34. Kim, H., Sandercock, T. G., & Heckman, C. J. (2015). An action potential-driven model of soleus 602 muscle activation dynamics for locomotor-like movements. J. Neural Eng., 12(4), 046025. 22 603 35. Krylow, A.M. & Sandercock, T.G. (1997). Dynamic Force Responses of Muscle Involving Eccentric 604 Contraction. J. Biomech., 30 (1), 27–33. 605 36. Lee, S.S.M., de Boef Miara, M., Arnold-Rife, A., Biewener, A.A., Wakeling, J.M. (2011). EMG analysis 606 tuned for determining the timing and level of activation in different motor units. . J. Electromyogr. 607 Kinesiol., 21:557–565. 608 37. Lee, S. S., de Boef Miara, M., Arnold, A. S., Biewener, A. A., & Wakeling, J. M. (2011). EMG analysis 609 tuned for determining the timing and level of activation in different motor units. J. Electromyogr. 610 Kinesiol., 21(4), 557-565. 611 38. Lee, S. S., Arnold, A. S., de Boef Miara, M., Biewener, A. A., & Wakeling, J. M. (2013). Accuracy of 612 gastrocnemius muscles forces in walking and running goats predicted by one-element and two613 element Hill-type models. J. Biomech., 46(13), 2288-2295. 614 39. Levin, O., Mizrahi, J., Isakov, E. (1999). Transcutaneous FES of paralyzed quadriceps: is knee torque 615 affected by unintended activation of the hamstrings. J. Electromyogr. Kinesiol., 10:47–58 616 40. Lichtwark, G. A., & Wilson, A. M. (2005). A modified Hill muscle model that predicts muscle power 617 output and efficiency during sinusoidal length changes. J. Exp. Biol., 208(15), 2831-2843. 618 41. Luff, A.R. (1975). Dynamic properties of fast and slow skeletal muscles in the cat and rat following 619 cross-reinnervation. J. Physiol., 248: 83-96. 620 42. Maganaris, C.N, Baltzopoulos, V., Sargeant A.J. (1998). In vivo measurements of the triceps surae 621 complex architecture in man: implications for muscle function. J. Physiol., 512:603–614. 622 43. Millard, M., Uchida, T., Seth, A., & Delp, S. L. (2013). Flexing computational muscle: modeling and 623 simulation of musculotendon dynamics. J. Biomed. Eng., 135(2), 021005. 624 44. Milner‐Brown, H. S., & Stein, R. B. (1975). The relation between the surface electromyogram and 625 muscular force. J. Physiol., 246(3), 549-569. 626 45. Mörl, F., Siebert, T., Schmitt, S., Blickhan, R., & Guenther, M. (2012). Electro-mechanical delay in 627 Hill-type muscle models. J. Mech. Med. Biol., 12(05), 1250085. 628 46. Morrison, S. M., Dick, T. J., & Wakeling, J. M. (2015). Structural and mechanical properties of the 629 human Achilles tendon: Sex and strength effects. J. Biomech., 48(12):3530-3. 630 doi:10.1016/j.jbiomech.2015.06.009 23 631 47. O'Neill, M. C., Lee, L. F., Larson, S. G., Demes, B., Stern, J. T., & Umberger, B. R. (2013). A three632 dimensional musculoskeletal model of the chimpanzee (Pan troglodytes) pelvis and hind limb. J. 633 Exp. Biol., 216(19), 3709-3723. 634 48. Otten, E. (1987). A myocybernetic model of the jaw system of the rat. J. Neurosci. Meth., 21: 287635 302 636 49. Perreault, E.J., Heckman, C.J., Sandercock, T.G. (2003). Hill muscle model errors during movement 637 are greatest within the physiologically relevant range of motor unit firing rates. J. Biomech., 36:211– 638 218. 639 50. Peterson, C. L., Hall, A. L., Kautz, S. A., & Neptune, R. R. (2010). Pre-swing deficits in forward 640 propulsion, swing initiation and power generation by individual muscles during hemiparetic walking. 641 J. Biomech., 43(12), 2348-2355. 642 51. Prager, R. W., Rohling, R. N., Gee, A. H., & Berman, L. (1998). Rapid calibration for 3-D freehand 643 ultrasound. Ultrasound Med. Biol., 24(6), 855-869. 644 52. Rack, P.M.H., Westbury DR. (1969). The effects of length and stimulus rate on tension in the 645 isometric cat soleus muscle. J. Physiol., 204, 443–460. 646 53. Rassier, DE, MacIntosh BR, Herzog W. (1999). Length dependence of active force production in 647 skeletal muscle. J. Appl. Physiol. 86, 1445–1457 648 54. Reaz, M. B. I., Hussain, M. S., & Mohd-Yasin, F. (2006). Techniques of EMG signal analysis: detection, 649 processing, classification and applications. Biol. Proced. Online, 8(1), 11-35. 650 55. Rockenfeller, R., Günther, M., Schmitt, S., Götz, T. (2015). Comparative Sensitivity Analysis of 651 Muscle Activation Dynamics. Comput Math Methods Med., 2015, 585409, 652 doi:10.1155/2015/585409 653 56. Ross, S.A., & Wakeling, J.M. (2016). Muscle shortening velocity depends on tissue inertia and level 654 of activation during submaximal contractions. Biol. Lett. 12 20151041; DOI: 10.1098/rsbl.2015.1041. 655 57. Rouffet, D.M., Hautier, C.A. (2008). EMG normalization to study muscle activation in cycling. J. 656 Electromyogr. Kinesiol., 18, 866-878. 657 58. Roy, R. R., Meadows, I. D., Baldwin, K. M., & Edgerton, V. R. (1982). Functional significance of 658 compensatory overloaded rat fast muscle. J. Appl. Physiol., 52(2), 473-478. 24 659 59. Sandercock, T.G. & Heckman, C.J. (1997). Force from cat soleus muscle during imposed locomotor660 like movements: experimental data versus Hill-type model predictions. J. Neurophysiol., 77:1538– 661 1552. 662 60. Shue, G., Crago, E., Chizeck, H.J. (1995). Muscle-joint models incorporating activation dynamics, 663 moment-angle, and moment-velocity properties. IEEE Trans. Biomed. Eng., 42:212–223. 664 61. Spaegele, T., Kistner, A., Gollhofer, A. (1999). Modelling simulation and optimization of human 665 vertical jump. J. Biomech., 32:521–30. 666 62. Spector, S.A., P.F. Gardiner, R.F. Zernicke, R.R. Roy, & V.R. Edgerton. (1980) Muscle architecture and 667 force-velocity characteristics of cat soleus and medial gastrocnemius: implications for motor 668 control. J. Neurophysiol., 44: 951-960. 669 63. Steele, K. M., Seth, A., Hicks, J. L., Schwartz, M. S., & Delp, S. L. (2010). Muscle contributions to 670 support and progression during single-limb stance in crouch gait. J. Biomech., 43(11), 2099-2105. 671 64. Umberger, B. R., Gerritsen, K. G. & Martin, P. E. (2003). A model of human muscle energy 672 expenditure. Comput. Methods Biomech. Biomed. Engin. 6, 99-111. 673 65. Vandervoort, A.A., McComas, A.J. (1983). A comparison of the contractile properties of the human 674 gastrocnemius and soleus muscles. Eur. J. Appl. Physiol., 51:435–440. 675 66. Von Tscharner, V. (2000). Intensity analysis in time-frequency space of surface myoelectric signals 676 by wavelets of specified resolution. J. Electromyogr. Kinesiol., 10:433–445. 677 67. Von Tscharner, V., & Goepfert, B. (2006). Estimation of the interplay between groups of fast and 678 slow muscle fibers of the tibialis anterior and gastrocnemius muscle while running. J. Electromyogr. 679 Kinesiol., 16(2), 188-197. 680 68. Wakeling, J.M., Rozitis, A.I. (2004). Spectral properties of myoelectric signals from different motor 681 units in the leg extensor muscles. J. Exp. Biol., 207(Pt 14):2519–28. 682 69. Wakeling, J.M., Uehli, K., Rozitis, A.I. (2006). Muscle fibre recruitment can respond to the mechanics 683 of the muscle contraction. Interface, 3:533–544. 684 70. Wakeling J.M. (2009). Patterns of motor recruitment can be determined using surface EMG. J. 685 Electromyogr. Kinesiol., 19(2):199–207 25 686 71. Wakeling, J. M. & Horn, T. (2009). Neuromechanics of muscle synergies during cycling. J. 687 Neurophysiol., 101, 843-854. 688 72. Wakeling, J. M., Blake, O. M., Wong, I., Rana, M., & Lee, S. S. (2011). Movement mechanics as a 689 determinate of muscle structure, recruitment and coordination. Phil. Trans. R. Soc. B., 366(1570), 690 1554-1564. 691 73. Wakeling, J.M., Lee, S.S.M., de Boef Miara, M., Arnold, A.S., Biewener, A.A. (2012). A muscle’s force 692 depends on the recruitment patterns of its fibers. Ann. Biomed. Eng., 40:1708–1720. 693 74. Ward, S. R., Eng, C. M., Smallwood, L. H., & Lieber, R. L. (2009). Are current measurements of lower 694 extremity muscle architecture accurate?. Clin. Orthop. Relat. Res., 467(4), 1074-1082. 695 75. Winters, T. M., Takahashi, M., Lieber, R. L., & Ward, S. R. (2011). Whole Muscle Length-Tension 696 Relationships Are Accurately Modeled as Scaled Sarcomeres in Rabbit Hindlimb Muscles. J. 697 Biomech., 44, 109–115. 698 76. Wolf, S. L., Ammerman, J., & Jann, B. (1998). Organization of responses in human lateral 699 gastrocnemius muscle to specified body perturbations. J. Electromyogr. Kinesiol., 8(1), 11-21. 700 77. Zajac, F.E. (1989). Muscle and tendon: properties, models, scaling, and application to biomechanics 701 and motor control. Crit. Rev. Biomed. Eng., 17:359–411. 702 703 704 705 706 707 708 709 710 711 712 713 714 26 715 Tables 716 717 718 719 Table 1. Centre frequencies and scale factors of optimized wavelets. 720 Muscle fc low [Hz] slow fc high [Hz] shigh LG 72.87 0.085 156.69 0.078 MG 62.32 0.078 146.95 0.051 721 722 Centre frequencies fc low and fc high and scale factors slow and shigh that characterize the low and high 723 frequency components from the EMG intensity spectra for the muscle-specific optimized wavelets (Eq. 724 9) for the lateral (LG) and medial (MG) gastrocnemii. 725 726 727 728 729 Table 2. Constants for the excitation-activation transfer functions to determine muscle active state from 730 EMG intensities 731 732 Constant Fast Excitation Slow Excitation Total Excitation act [ms] 25 45 35 act 0.6 0.6 0.6 733 act is the activation time constant and act is the ratio of the activation to deactivation time. 734 To derive these constants, we manually traced (ImageJ, NIH, Maryland, USA) tetanic curves for cat slow 735 and fast motor units (Kernell et al., 1983) and optimized the parameters act and act to reproduce the 736 force traces for a matched stimulation duration. 737 738 739 740 741 27 Table 3. Average lateral (LG) and medial gastrocnemii (MG) coefficient of determination, r2 and RMSE between model predictions and ultrasound-based estimates of force, using the one-element and the two-element models across pedalling conditions. LG MG Cadence (r.p.m.) crank load (N m) r2 r2 RMSE RMSE r2 r2 RMSE one-element two-element (% Fmax) (% Fmax) one-element two-element (% Fmax) one-element two-element one-element RMSE (% Fmax) two-element 60 44 0.54 ± 0.07 0.52 ± 0.07 8.98 ± 0.86 9.38 ± 0.95 0.79 ± 0.03 0.78 ± 0.03 10.17 ± 1.60 80 14 0.46 ± 0.07 0.47 ± 0.07 10.98 ± 0.83 11.32 ± 1.05 0.64 ± 0.05 0.65 ± 0.05 10.60 ± 0.99 80 26 0.61 ± 0.03 0.62 ± 0.04 10.27 ± 1.47 10.45 ± 1.40 0.66 ± 0.05 0.66 ± 0.05 10.83 ± 1.20 80 35 0.46 ± 0.07 0.46 ± 0.07 10.88 ± 1.20 11.46 ± 1.09 0.68 ± 0.05 0.68 ± 0.05 11.96 ± 1.49 80 44 0.55 ± 0.08 0.55 ± 0.08 13.05 ± 2.33 13.78 ± 1.95 0.71 ± 0.05 0.71 ± 0.05 12.39 ± 1.86 100 26 0.41 ± 0.08 0.42 ± 0.09 15.47 ± 2.44 14.90 ± 2.36 0.50 ± 0.07 0.51 ± 0.07 18.25 ± 2.03 120 13 0.38 ± 0.08 0.40 ± 0.09 16.49 ± 1.99 15.78 ± 1.73 0.45 ± 0.05 0.48 ± 0.05 16.93 ± 2.09 140 13 0.29 ± 0.06 0.32 ± 0.05 21.91 ± 3.50 21.22 ± 3.36 0.46 ± 0.08 0.49 ± 0.08 17.53 ± 2.25 Average 0.46 0.47 13.50 13.51 0.61 0.63 13.65 Data are presented as mean ± SE (n=16); r2 and RMSE were calculated for the time-varying forces across one complete pedal revolution. 9.86 ± 1.26 10.77 ± 0.89 10.67 ± 1.04 11.36 ± 1.40 12.09 ± 1.59 18.26 ± 1.82 17.45 ± 2.01 17.43 ± 2.25 13.48 28 Table 4. Model and Equation Parameters. Parameter Definition ̂() Normalized activation Scaling factor, represents the lower EMG 1 excitations that would be expected from action potentials with higher spectral frequencies ̂() Normalized excitation LG Lateral gastrocnemius force MG Medial gastrocnemius force ̂a(f) ̂a() m max ̂p(f) SEE−LG,T SEE−MG,T SEE−LG SEE−MG AT−LG AT−MG 0,AT−LG 0,AT−MG f 0,f ̂ Normalized active force-length relation Normalized force-velocity relation Muscle Force Maximal isometric force Normalized passive force-length relation LG toe region tendon stiffness MG toe region tendon stiffness LG linear region tendon stiffness MG linear region tendon stiffness LG tendon length MG tendon length LG tendon slack length MG tendon slack length Fascicle length Fascicle slack length Normalized fascicle velocity ̂0 Maximum intrinsic speed m Muscle volume α Curvature of force-velocity relation Pennation angle 0 Maximum isometric stress of muscle fibres Source Derived from measured EMG (see Eq. 10) Calculated, used a value of 2 Calculated from measured EMG Calculated from tendon length, stiffness, and slack length (see Eq. 1) Calculated from tendon length, stiffness, and slack length (see Eq. 11) Estimated from literature (Lee et al., 2013; see Eq. 4) Literature (see Eq. 7 and 8) Calculated from model (Eq. 1) Estimated (see Eq. 3) Estimated from literature (Millard, 2013; Eq. 5 and 6) Measured in Dick et al., 2016 Measured in Dick et al., 2016 Measured in Dick et al., 2016 Measured in Dick et al., 2016 Measured using B-mode ultrasound Measured using B-mode ultrasound Estimated from AT−LG at 310 ̊ of pedal cycle Estimated from AT−MG at 310 ̊ of pedal cycle Measured using B-mode ultrasound Determined from f at 0,AT Derived from f Estimated from literature (5 and 10 −1 for slow and fast fibres) Estimated using the regression equations of Handsfield et al. (2014) Estimated from literature (0.18 and 0.29 for slow and fast muscle fibres; Wakeling et al., 2012; Lee et al., 2013) Measured using B-mode ultrasound Estimated from literature (225 kPa; Spector et al., 1980; Roy et al., 1982) 29 Figure Legends Figure 1. Approach for comparing LG and MG forces during cycling estimated from ultrasound-based measures of tendon length changes and stiffness and predicted from Hill-type muscle models. During the cycling protocol (A), subjects pedalled on a stationary bike while we measured tendon lengths from tracked ultrasound images of the LG and MG muscle-tendon junctions, fascicle lengths and pennation angles from ultrasound images of the LG and MG muscle bellies, and activation patterns derived from surface EMG. A trigger from the ultrasound system was used to synchronize all data. Forces were estimated from the tracked ultrasound images of tendon length changes during cycling and tendon stiffnesses measured for each subject in an isometric protocol (Dick et al., 2016). The experimentally determined fascicle lengths f, pennation angles β, and normalized muscle activations ̂() for total (black), slow (red), and fast (blue) motor units were used as inputs for the muscle models (B). We tested a traditional one-element Hill-type muscle model and additionally a two-element model that accounted for the independent contributions of slow and fast muscles fibres (C). Estimated forces from the ultrasound-based approach were compared to the forces predicted from the models for the LG and MG. CE: contractile element, PEE: parallel elastic element, β: pennation. 30 Figure 2. Myoelectric intensity spectra reconstructed from the pooled frequency spectra (thin lines) and optimized wavelets (thick lines) for human LG (A) and MG (B). Low frequency spectra are shown in red and high frequency spectra in blue. Figure 3. Time-varying force profiles for the LG and MG at 60 r.p.m at 44 N M (A) and 100 r.p.m. at 26 N m for one representative subject. Ultrasound-based estimates of force are represented in black, predicted forces from the one-element model in grey, and from the two-element model in red. Muscle forces m are normalized to the maximum isometric force max of either the LG or the MG. 31 Figure 4. Differential recruitment of slow and fast muscle fibres with pedalling condition. Raw EMG (grey) and total (black), slow (red), and fast (blue) intensity and activation traces for one representative subject LG (A) and MG (B) pedalling at 60 r.p.m. at 44 N m (left panel) and 140 r.p.m. at 13 N m (right panel). Vertical dashed lines show timing of pedal top-dead-centre. These conditions have been chosen to highlight the differences in recruitment between slow and fast muscle fibres that can occur across the range of mechanical conditions tested here. 32 Figure 5. Sensitivity of one-element model predictions to force-velocity parameters: curvature of forcevelocity relation α (A), and maximum intrinsic speed 0 (B) in the LG and MG. Data points shown represent the mean ± SE across all subjects (n=16). Different values of α and 0 are shown using different shades of grey. 33 Supplementary Files: Supplementary Table 1. Characteristics of the 16 cyclists tested. Subject Sex Age (y) Height (cm) Body Mass (kg) LG Max Isometric Force (N)* MG Max Isometric Force (N)* LG tendon Stiffness (N mm-1)† MG tendon Stiffness (N mm-1)† 1 M 29 186 90.85 674 1179 34.5 63.2 2 M 28 167 68.00 588 1028 22.0 38.9 3 M 33 175 68.00 546 956 36.2 50.0 4 M 32 183 63.00 488 854 22.5 39.0 5 M 37 173 71.80 599 1047 39.2 55.0 6 M 30 179 70.30 555 972 24.5 41.8 7 M 47 183 86.50 661 1157 42.7 59.1 8 M 31 183 82.80 676 1184 26.3 41.2 9 F 30 167 65.80 584 1022 34.5 73.6 10 F 43 174 64.40 552 966 28.9 55.4 11 F 24 171 65.80 571 999 27.8 52.0 12 F 26 167 58.00 500 874 22.1 47.2 13 F 22 160 51.30 474 830 25.3 43.4 14 F 23 168 68.00 608 1065 29.5 58.8 15 F 20 165 67.30 622 1089 26.9 42.4 16 F 32 167 59.40 541 948 38.2 59.2 * Estimates of the muscles’ maximum isometric force-generating capacity based on the muscles’ volumes (Handsfield et al., 2014) and optimal fibre lengths (Arnold et al., 2010), assuming a maximum isometric muscle stress 0 of 225 kPa. (Spector et al., 1980; Roy et al., 1982). †Linear region LG and MG AT stiffnesses were estimated as the proportion of the total AT stiffness contributed by each muscle based on ratio of the maximum force-generating capacity max of either the LG or MG to the combined triceps surae maximum force. 34 Supplementary Table 2. Chi-square test results for comparison of r2 and RMSE between one-element and two-element Hill-type muscle models. Coefficient of determination, r2 RMSE Cadence (r.p.m.) Load (N m) two-element (% total) p-value two-element (% total) p-value 60 44 46.9 0.73 37.5 0.15 80 14 59.4 0.29 28.1 0.01 80 26 62.5 0.16 43.8 0.47 80 35 65.6 0.07 43.8 0.47 80 44 59.4 0.29 46.9 0.72 100 26 71.9 0.01* 62.5 0.15 120 13 68.8 0.03* 68.8 0.03* 140 13 71.9 0.01* 68.8 0.03* Values in % total column indicate the number of times for the combined MG and LG results across individual subjects and pedalling conditions, as a %, that the two-element model performed better (higher r2 and lower RMSE) (n=16). Differences were considered significant at the p<0.05 level. *indicates p<0.05 Supplementary Figure 1. Experimental set up displaying B-mode ultrasound probe positioned over the LG muscle belly to image muscle fascicles and LED motion capture markers to record kinematics (A) and surface EMG electrodes placed on contralateral limb to record simultaneous muscle excitations (B). 35 Supplementary Figure 2. B-mode ultrasound images of the MG muscle belly (A) and MG muscle-tendon junction (B). Red dots indicate digitized points that were used to determine time-varying fascicle lengths and pennation angles from (A) and tendon lengths from (B). 36 Supplementary Figure 3. Model predictions differ between muscles, pedalling conditions, and model type. Comparison of the one-element (grey) and two-element (black) Hill-type models using r2 (A) and RMSE (B) across the pedal cycle for the LG and MG. Data points shown represent the mean ± SE across all subjects at each pedalling condition. 37