Simulated work loops predict maximal human cycling power

Several species, including fish, birds and lizards, perform some maximal locomotor activities with coordination patterns that maximize muscle power (Askew and Marsh, 2002; Askew et al., 2001; Curtin et al., 2005; Franklin and Johnston, 1997; James and Johnston, 1998; Syme and Shadwick, 2002; Wakeling and Johnston, 1998). That is, muscle power for a complete shortening–lengthening cycle during voluntary movement is at or near the maximum possible for that muscle, even when a large parameter space is evaluated using in situ or in vitro work loops. For example, previous authors have reported that muscle power is maximized during escape responses (Curtin et al., 2005; Franklin and Johnston, 1997; James and Johnston, 1998; Wakeling and Johnston, 1998) and steady-state swimming in fish (Syme and Shadwick, 2002), and flight take off in quail (Askew et al., 2001). Because power for a shortening–lengthening cycle arises from complex interactions of force–length, force–velocity and activation/deactivation characteristics (Josephson, 1999), these findings suggest that an animal's movement patterns develop in concert with muscle characteristics so as to maximize muscle power.

Most investigations in which in vivo voluntary movements have been compared with muscle contractions measured through in situ work loops have focused on studying movements performed dominantly by one or two muscles (e.g. Biewener and Corning, 2001). This approach has allowed scientists to evaluate important functional movements while instrumenting and dissecting only the few dominant muscle(s). However, this approach is problematic for studying many movements, particularly locomotor movements that involve multiple muscles (including biarticular muscles) spanning multiple joints. Studying such complex movements in situ is difficult because of the surgical complexity of instrumenting all of the relevant muscles. Consequently, complex locomotor movements have not been studied by comparing in vivo and in situ work loops, and the extent to which these locomotor activities are performed with maximized muscle power remains unknown.

Understanding whether muscle power is maximized during complex mammalian movements, and human movements in particular, is important for studying basic aspects of motor control. Notably, such understanding could clarify why some biarticular muscles appear to perform contradictory actions (Lombard's paradox: Andrews, 1987; Gregor et al., 1985). For example, the biceps femoris long head is anatomically positioned to both extend the hip and flex the knee, but is active during whole-leg extension; thus, this muscle appears to produce the desired action (extension) at the hip, but a counterproductive action (flexion) at the knee. Understanding the role of biarticular muscles could provide unique insight into voluntary control of whole-limb movement. Gaining such insight by performing experiments using in vivo and in situ techniques is not feasible for human muscles, but mathematical modeling could facilitate similar comparisons. Indeed, mathematical muscle models (e.g. Millard et al., 2013; Thelen, 2003; Winters, 1995) have been used to study how individual muscle actions contribute to complex activities such as walking (e.g. Anderson and Pandy, 2003; Buchanan et al., 2004; Piazza, 2006; Steele et al., 2010; Thelen and Anderson, 2006; Zajac et al., 2002), running (e.g. Dorn et al., 2012; Hamner et al., 2010; Lloyd and Besier, 2003) and cycling (e.g. Rankin and Neptune, 2008; van Soest and Casius, 2000; Yoshihuku and Herzog, 1990). Within the context of maximized power, a muscle model could be subjected to any specified length trajectory and stimulation onset and offset timing could be set to maximize work for a complete shortening–lengthening cycle. That is, a muscle model could be used to form a simulated work loop with realistic length trajectory and this could then be compared with experimental data recorded during maximal-effort human movement.

One human locomotor action that might be performed with maximized muscle power is maximal cycling. Indeed, we previously reported that overall maximal cycling power, measured at the level of the cranks, exhibited characteristics similar to power produced during maximized in situ work loops (Martin, 2007, 2000). However, to what extent power is maximized at the level of the joints and muscles during cycling remains an open question. Therefore, the aim of this investigation was to determine whether humans maximize muscle power during maximal voluntary cycling within the constraints imposed by the cycling action. To accomplish this aim, we developed simulations of work loops for the leg muscles using a mathematical muscle model (Thelen, 2003) with cycling-specific length trajectories. We compared the work-loop simulation results with experimentally measured power produced by humans performing maximal cycling. We specifically examined power production at the level of the joints and muscles by testing three hypotheses. First, given our previous work demonstrating similar characteristics between maximal cycling power at the level of the cranks and power produced during work loops, we hypothesized that the net function of the leg muscles crossing the hip and knee would exhibit similar power production to that observed during maximal cycling by human cyclists. Second, given that the ankle's primary purpose may not be to maximize power generation, but rather to transfer power delivered by the hip and knee to the pedal (Zajac et al., 2002), we hypothesized that the experimental and modeled ankle power would not agree as closely as power at the hip and knee. Finally, we hypothesized that biarticular muscles might produce joint power that differed substantially from muscle power, thus providing novel insight into biarticular muscle function (e.g. Lombard's paradox).

To determine whether humans maximize muscle power during maximal voluntary cycling, we compared joint power values measured during a maximal cycling activity with those derived from simulated work loops of muscles having parasagittal action in the lower extremity (Fig. 1). Experimental cycling data, including limb kinematics and pedal reaction forces, were collected during maximal isokinetic cycling at a pedaling rate of 120 rev min⁻¹, a pedaling rate generally associated with maximum power (Martin et al., 1997). The Thelen muscle model (Thelen, 2003) was used to simulate cyclic contractions of lower limb muscles with each constrained to the length trajectory imposed by the experimental cycling kinematics.

Previously reported kinematics (joint angle and angular velocity) and kinetics (net joint moment and power) during maximal cycling (Martin and Brown, 2009) were used for this investigation. To briefly summarize the experimental study, 13 highly trained cyclists (1 female, 12 males, mean±s.d. mass 74.8±6.5 kg) performed maximal isokinetic cycling trials at 120 rev min⁻¹ for one 30 s trial. For this investigation, we used only data from the first complete cycle for each subject, which represents a non-fatigued state at a constant cycling velocity. During each trial, pedal reaction force, pedal and crank angle, and limb segment position were recorded at 240 Hz. Specifically, pedal reaction force was recorded from the right pedal using two 3-component piezoelectric force transducers (Kistler 9251, Kistler USA, Amherst, NY, USA). Pedal and crank angle were recorded using digital encoders (S5S-1024-IB, US Digital, Vancouver, WA, USA) attached to the right pedal and crank. Limb segment position, defined as the position of the hip, knee, ankle and fifth metatarsal head, was derived from measurements from an instrumented spatial linkage system (Martin et al., 2007). The limb segment positions were used to calculate ankle, knee and hip joint angle and joint angular velocity (termed ‘experimental joint angle’ and ‘experimental joint angular velocity’ in Fig. 1). Parasagittal plane net joint moment values (termed ‘experimental joint moment’ in Fig. 1) at the ankle, knee and hip were determined using inverse dynamic techniques (Elftman, 1939). Joint power (termed ‘experimental joint power’ in Fig. 1) was calculated as the product of net joint moment and joint angular velocity. Net power (termed ‘experimental leg power’ in Fig. 1) was calculated as the sum of hip, knee and ankle joint power. Joint angle, joint angular velocity and joint power for the hip, knee and ankle from these experimental data are presented graphically.

To estimate maximal muscle power during cycling, we created 38 work-loop simulations, one for each lower extremity muscle with parasagittal plane actions (flexion and extension). Importantly, the work-loop simulations represent the muscular work generated during one complete pedal revolution. The kinematics of one pedal revolution matched the mean hip, knee and ankle joint angles measured across all cycling participants.

The inputs to the work-loop simulations were muscle–tendon length, velocity and moment arm trajectory (Fig. 1), which were estimated from a musculoskeletal model of the lower extremity (OpenSim, 3DGaitModel2392; Delp et al., 2007). Specifically, the musculoskeletal model, including muscle–tendon parameters, was scaled to match the mean segment lengths across all participants in the cycling experiment, and the experimentally measured hip, knee and ankle joint angles were input into the scaled model. Given that the experimental data only measured parasagittal plane motion (flexion/extension at the hip, knee and ankle), all other degrees-of-freedom (e.g. abduction/adduction and internal/external rotation) were held constant in a neutral position, with the exception of pelvic tilt. Pelvic tilt, which describes the position of the trunk relative to the thigh and will influence simulated hip joint angle and muscle length, was estimated by matching model and experimental kinematics and found to be −3 deg, indicating a slight forward lean of the trunk. Based on a scaled musculoskeletal model and the experimentally prescribed kinematics, muscle–tendon length and moment arm trajectories were calculated as a function of crank angle for 38 muscles with parasagittal actions at the hip, knee and ankle in the right limb (muscle names and abbreviations are summarized in Table 1). To present the kinematic muscle data, the following parameters have been summarized: maximum and minimum muscle–tendon length (relative to resting muscle length, where resting muscle length equals the tendon slack length plus the product of the optimal fiber length and the cosine of the pennation angle), moment arms, crank angles representing muscle shortening, and shortening velocities.

To perform the work-loop simulations, the muscle–tendon length, velocity and moment arm trajectories were input into a mathematical muscle model and the onset and offset of muscle stimulations were optimized to maximize power generation. For the mathematical muscle model, we specifically used the mathematical description provided by C. T. John (http://simtk-confluence.stanford.edu:8080/download/attachments/2624181/CompleteDescriptionOfTheThelen2003MuscleModel.pdf; https://simtk-confluence.stanford.edu:8443/display/OpenSim24/Muscle+Model+Theory+and+Publications) to develop custom-written code (Microsoft Excel 2013) of the Thelen (2003) muscle model. This muscle model includes differential equations describing the activation and deactivation dynamics that occur during muscle contraction. To derive muscle force for a given level of muscle activation, forward integration is required. To avoid the computational instability often associated with numerical integration, we used small time steps (0.042 ms, which is equivalent to a sampling frequency of 24 kHz). This provided stability for all 38 muscles when initial conditions were set within the passive lengthening phase. Given that the sampling rate of our input data (muscle–tendon length, velocity and moment arm trajectories) matched the 240 Hz sampling rate of our experimental data, we used a fourth-order Fourier series to resample the data at the required 24 kHz. This order for the Fourier series approximations agreed well with raw muscle–tendon length [mean±s.d. root mean square (RMS) error 0.02±0.01% of mean length] and moment arm (RMS error 0.06±0.2% of mean moment arm) trajectories. Each muscle was simulated individually in order to incorporate muscle-specific definitions of maximum isometric force, force–velocity shape, pennation angle, optimal fiber length and tendon slack length into the mathematical muscle model. All muscle parameters were defined to match those in the scaled musculoskeletal model. For all muscles, activation and deactivation time constants were defined as 10 and 40 ms (Winters and Stark, 1985), respectively. To maximize muscle power in the simulation, we optimized the stimulation onset and offset timing of each muscle. Specifically, onset and offset timing were selected to maximize net work and average power for complete shortening–lengthening cycles for each muscle. This is common practice in work-loop experiments and we sought to replicate that using our mathematical muscle model.

The outputs of the work-loop simulations were muscle force, muscle power, muscle joint moment, muscle joint power, net joint power and net leg power (Fig. 1). Muscle–tendon force (termed ‘simulated muscle force’ in Fig. 1) was directly derived from the mathematical muscle model based on each individual muscle's force–length and force–velocity characteristics. Muscle power (termed ‘simulated muscle power’ in Fig. 1) was calculated as the product of absolute muscle–tendon force and muscle–tendon velocity. Muscle joint moment (termed ‘simulated muscle joint moment’ in Fig. 1) for each muscle was calculated as the product of the muscle force and muscle–tendon moment arm; muscle joint moments were separately calculated at each joint crossed by a given muscle. Muscle joint power (termed ‘simulated muscle joint power’ in Fig. 1) was calculated from muscle joint moment and the joint angular velocity (from experimental data). Net joint power (termed ‘simulated joint power’ in Fig. 1) was calculated at each joint as the sum of the muscle joint power values at that joint. Net leg power (termed ‘summed muscle power’ in Fig. 1) was calculated as the sum of the power produced by all 38 muscles across all joints. Muscle stimulation onset and offset, peak and average force, net, positive and negative work, and peak and average muscle and joint power are reported to characterize these simulation results.

To test our first two hypotheses that humans perform maximal cycling with maximized muscle power at the hip and knee (hypothesis 1) but not at the ankle (hypothesis 2), we performed intraclass correlation and Pearson's correlation analyses of simulated versus experimentally measured power values throughout the pedaling cycle. Comparisons included summed muscle power versus experimental leg power, as well as simulated joint power versus experimental joint power. Intraclass correlation provides a quantitative assessment of the agreement between each set of measures, while Pearson's correlation provides a measure of similarity of shape without regard to amplitude. To test our third hypothesis, we explored biarticular muscle function by comparing simulated muscle and joint power. Specifically, we used intraclass correlations to compare simulated muscle power with (a) the simulated joint power for each joint spanned by the uniarticular muscle, (b) the simulated joint power for each joint spanned by the biarticular muscle and (c) the sum of those two joint power values. We expected that power at either joint would not agree with muscle power, but that the sum of the power at the two joints would match that of the muscle. Further, we expected that power at both joints would exhibit substantial negative power, while the muscle would actually produce very little negative power. This would underscore the importance of considering both joints spanned by biarticular muscles.

Experimental joint angle and angular velocity values exhibited clear extension and flexion phases within each crank cycle (Fig. 2). The knee exhibited the greatest range of motion (Fig. 2A) and angular velocity (Fig. 2B), followed by the hip and ankle. The hip (357–184 deg, negative angular velocity) and knee (339–166 deg, positive angular velocity) were in extension for 187 deg of crank rotation, whereas the ankle was in extension for 210 deg of crank rotation (51–261 deg, negative angular velocity). The hip and ankle joints produced substantial power (Fig. 2C) during extension (448 and 141 W, respectively), with minimal power during flexion (20 and −15 W, respectively), whereas the knee joint produced substantial power in both extension (215 W) and flexion (188 W).

To provide example traces of muscle–tendon length, velocity and moment arms for a representative uniarticular (vastus lateralis, VL) and biarticular (biceps femoris long head, BFLH) muscle, we plotted those values against crank angle (Fig. 3). The uniarticular VL exhibited a clear shortening/lengthening pattern, whereas the biarticular BFLH remained nearly isometric for approximately 25% of the cycle (Fig. 3A). VL reached a peak shortening velocity (negative value) of ∼2.5 fiber lengths s⁻¹, whereas peak shortening velocity of BFLH was ∼1.0 fiber length s⁻¹ (Fig. 3B). The biarticular BFLH exhibited shortening during portions of the crank cycle involving both hip extension (<184 deg) and knee flexion (>166 deg). Moment arms for VL and BFLH exhibited large variation across the cycle and moment arms for BFLH were substantially different at the proximal and distal joints (Fig. 3C).

For the entire muscle set, maximum and minimum muscle–tendon length was 94±11% and 86±13% (Table 1) of resting length, respectively. Muscles began and ended shortening at a wide variety of crank angles depending on the joint(s) spanned and the primary action (Table 1). Maximum and minimum values of muscle tendon moment arms were 40±20 and 26±18 mm, respectively (Table 1). Peak and average shortening velocities were 1.59±1 and 0.89±0.52 fiber lengths s⁻¹, respectively (Table 1).

To illustrate force and power production characteristics, we plotted those measures against crank angle for VL and BFLH (Fig. 4A,B). Active force production began slightly before muscle shortening and continued into the lengthening phase, with a peak value closely following the onset of shortening when velocity was small. Peak power occurred near the midpoint of shortening when shortening velocity was near its peak (compare peak velocity in Fig. 3B with peak power in Fig. 4B). We also plotted force against muscle length to form a modeled work loop (Fig. 4C). These modeled work loops displayed the features of in situ work loops: the data progress counter-clockwise, the area under the top (concentric) portion of the trace represents positive work and the area under the lower trace represents negative work. For the entire muscle set, mean (±s.d.) values for muscle stimulation onset and offset that produced maximum work and power occurred at 17±6 ms (i.e. 1.7 times the activation time constant) prior to the beginning of muscle shortening and 49±8 ms (1.2 times the deactivation time constant) prior to the end of shortening, respectively. Average force during shortening and lengthening was 53±15% and 7±3% of isometric force, respectively (Table 2). These concentric and eccentric forces produced 8.8±8.6 J of positive work, 1.0±1.2 J of negative work and 7.8±7.4 J of net work (Table 2). Peak and average power was 56±51 and 16±15 W, respectively (Table 2).

Representative power values for VL and BFLH demonstrate characteristics of joint power produced by uniarticular and biarticular muscles (Fig. 5). Muscle and joint power were nearly identical for VL (Fig. 5B) as were muscle power and the sum of joint power at the hip and knee for BFLH (Fig. 5A). However, hip and knee joint power differed dramatically from BFLH muscle power. Peak and average muscle power of BFLH was, respectively, 49 and 12 W. In contrast, peak and average joint power of BFLH was, respectively, 150 and 21 W at the hip, and 76 and −8 W at the knee (Fig. 5A). These large differences between power produced by the muscle and power delivered to the joints support our hypothesis that simulations can be used to elucidate how biarticular muscles function at their proximal and distal joint.

Peak and average joint power produced by uniarticular muscles was closely related to muscle power [peak: r²=0.999, intraclass correlation coefficient (ICC)=0.9993 (ICC confidence limits: 0.9985–0.9997), average: r²>0.999, ICC=0.9995 (0.9989–0.9998); Table 2] with minor differences (peak: 1.1±1.7 W and average: 0.2±0.4 W) arising from estimations of muscle–tendon moment arms. Peak and average joint power produced by biarticular muscles was not closely related to their respective muscle power [peak: r²=0.24, ICC=0.29 (−0.12–0.65)], average: r²>0.08, ICC=0.22 (−0.16–0.61); Table 2]. However, peak [r²=0.988, ICC=0.999 (0.993–0.999)] and average [r²=0.998, ICC=0.993 (0.964–0.999)] muscle power agreed quite well with the sum of the joint power from both joints spanned by biarticular muscles.

To illustrate the individual muscle contributions to net joint power, we plotted the power produced by each muscle at the hip, knee and ankle (Fig. 6), as well as the net power of all the muscles spanning the joint (Fig. 7). Note that the net powers at the hip and knee are substantially influenced by negative joint powers produced by biarticular muscles as previously shown for BFLH.

Experimental leg power (534 W) was closely approximated by the sum of all modeled muscle power values [589 W, r²=0.91, ICC=0.91 (0.86–0.94); Fig. 7A]. Experimental joint power, calculated at 1 deg increments of crank angle, was also closely approximated by modeled joint power for the hip [Fig. 7B; r²=0.94, ICC=0.92 (0.79–0.96)] and knee [Fig. 7C; r²=0.90, ICC=0.95 (0.94–0.96)], but not the ankle [Fig. 7D, r²=0.89, ICC=0.74 (−0.09–0.92)]. These results provide strong support for our hypotheses that voluntary maximal cycling is performed with maximized muscle power at the hip and knee but less so at the ankle. When muscle power was averaged over the complete crank cycle, modeled joint power underestimated experimental joint power at the hip (190 versus 248 W), agreed well with experimental joint power at the knee (217 versus 208 W) and substantially overestimated experimental joint power at the ankle (179 versus 78 W).

Cycling, like many human locomotor activities, involves coordinated extension and flexion of the hip, knee and ankle, which are powered by uniarticular and biarticular muscles. Coordination strategies for controlling activation of these muscles may involve optimizing force direction, power transfer and/or power production (Zajac et al., 2002). In this study, we demonstrated that simulations which maximized power of muscles that cross the hip and knee closely approximated joint power measured experimentally during maximal voluntary cycling. This finding supports our hypothesis that humans maximize muscle power during voluntary maximal cycling, as do birds, fish and lizards during some maximal activities (Askew and Marsh, 2002; Askew et al., 2001; Curtin et al., 2005; Franklin and Johnston, 1997; James and Johnston, 1998; Syme and Shadwick, 2002; Wakeling and Johnston, 1998). Importantly, the optimization method implemented in this study only altered the timing of muscle activation and deactivation to maximize muscle power. Thus, these results imply that human motor control patterns optimize the timing of activation and deactivation to maximize power for complete shortening/lengthening contraction cycles. In contrast to those for the hip and knee, the experimental data were only modestly approximated by the simulations that maximized power of muscles that cross the ankle. This supports the notion that the primary function of the ankle muscles during maximal cycling at 120 rev min⁻¹ is energy transfer rather than energy production (Zajac et al., 2002). Finally, as discussed further below, our modeling provided novel insight by demonstrating that biarticular muscle power differed substantially from individual joint power.

The experimental data utilized in this investigation were obtained from competitive cyclists and, thus, might represent a highly skilled power production technique. However, findings from two previous investigations suggest that trained cyclists perform similarly to non-cyclists. First, Mornieux and colleagues (2008) demonstrated that cyclists and non-cyclists produced nearly identical pedal forces with two types of pedals as well as with and without visual feedback of power within the cycle. Second, Martin and colleagues (2000) previously reported that non-cyclists produced power equal to or slightly greater than that produced by trained racing cyclists with just 2 days of four rehearsal trials (3–4 s each). Taken together, these findings suggest that cycling provides a window through which to observe basic aspects of neuromuscular function and motor control. Thus, we believe that our results represent a global finding that innate extension and flexion patterns are capable of executing muscle stimulation patterns that maximize power within the context of a complete shortening/lengthening contraction cycle.

Maximizing power for a complete shortening/lengthening cycling requires a compromise of stimulating the muscle long enough (e.g. throughout a large portion of muscle shortening) to produce substantial positive power while ending stimulation early enough so as to prevent excessive eccentric work (Caiozzo and Baldwin, 1997). Negative work during lengthening averaged −12% of the work done during shortening, demonstrating the complex trade-offs of positive and negative work associated with stimulation timing to maximize muscle work and power. Our simulated stimulation patterns achieved this balance with onsets beginning an average of 17±6 ms prior to the beginning of shortening and offsets beginning an average of 49±8 ms prior to the end of shortening. With the model's exponential activation time constant of 10 ms, muscles were 82% activated as they began to shorten and thus produced near-maximum force. In the final 17 ms of lengthening, the muscle was nearly isometric, and only 5% of the net negative work resulted from this activation strategy. The majority of negative work occurred during lengthening after deactivation. With the model's deactivation time constant of 40 ms, deactivation at 49 ms prior to lengthening meant that muscles were 29% activated when they began lengthening and thus could produce substantial antagonistic force. Further, residual activation after stimulation offset caused muscles to be activated at >1% for 135 ms of the lengthening phase or 27% of the cycle. This residual activation produced 95% of the negative work. Thus, the overwhelming majority of negative work from our simulations was due to the lack of complete relaxation during lengthening, as has been described previously (e.g. Josephson, 1985).

Importantly, our stimulation onset and offset timing values agree reasonably well with previously reported electromyography (EMG) data with respect to cycle crank angle (Fig. 8). Specifically, Dorel et al. (2012) reported surface EMG data for 11 muscles of 15 trained cyclists during maximal cycling at 100 rev min⁻¹. We compared the optimized onset and offset timing from our simulation with their data (digitized values from fig. 5 in Dorel et al., 2012). Note that this is not an ideal comparison because simulated onset and offset timing represent muscle stimulation (neural command stimulating muscle), while recorded surface EMG represents muscle activation (muscle contraction already past a given threshold). Correlations demonstrate that EMG onsets agreed reasonably well (r²=0.94) with our simulated muscle onset timing for gluteus maximus (Gmax), tensor fasciae latae (TFL), rectus femoris (RF), vastus lateralis (VL), lateral gastrocnemius (LG), medial gastrocnemius (MG) and soleus (Sol), but differed substantially for semimembranosus (SM), biceps femoris long head (BFLH), vastus medialis (VM) and tibialis anterior (TA). EMG and simulated muscle offsets also agreed reasonably well (r²=0.82, r²=0.88 without TFL) for all muscles except TFL and Sol. Differences in these values for muscles that span the ankle (TA and Sol) likely reflect their role as stabilizers responsible for power transfer rather than direct power producers during maximal cycling. Indeed, examination of Figs 6E and 7D suggests that Sol power late in its shortening phase (crank angles of 183–255 deg) accounts for almost all of the differences in ankle power for crank angles greater than 180 deg. Maximized Sol power during that portion of the cycle may have imparted negative power to the pedal, reducing overall power. Differences in the values for biarticular muscles (SM and BFLH) may indicate that our simulations did not fully describe the actions of these muscles. For example, EMG of those cyclists indicated that they activated the biarticular SM and BFLH well before the onset of muscle shortening (Table 1). Consequently, these muscles may have produced large near-isometric force that delivered opposing moments at the hip and knee but no net power; such moments might suggest that these muscles perform a power transfer or kinematic role that was not clearly evident in our simulations. Alternatively, these biarticular hip extensors may simply have been activated synergistically with other hip extensors as a single muscle group. The difference in offset of TFL may reflect its thigh-abduction role, which could act to stabilize the pelvis, although frontal plane actions were not included in our model. Indeed, without TFL, the coefficient of determination for offset increased to r²=0.88. The difference in onset for VM is more difficult to explain but could be due to the lower pedaling rate adopted by Dorel et al. (2012) of 100 rev min⁻¹ because different pedaling rates may require different kinematic strategies (McDaniel et al., 2014).

Our data confirmed our hypothesis that biarticular muscles produce joint power that differed substantially from muscle power, thus providing novel insight into biarticular muscle function (Tables 1 and 2, and Figs 5 and 6). BFLH provides a compelling example. At a crank angle of 90 deg, muscle power of BFLH was 7 W but that small power manifested as a hip joint power of 129 W and a knee joint power of −122 W (Fig. 5A). These contrasting muscle, hip and knee power values occurred because the muscle shortening velocity was small, facilitating high muscle tension (∼500 N), while the hip and knee joints had substantial joint angular velocities of 362 and 258 deg s⁻¹, respectively. This combination of high force crossing moving joints produced these large power values, even while the muscle was nearly isometric and therefore producing almost no muscle power. When considered over the entire shortening/lengthening contraction cycle, BFLH produced 11.9 W of muscle power, 21 W of hip joint power and −8.1 W of knee joint power (Table 2). These examples of instantaneous and average power demonstrate the importance of considering effects at proximal and distal joints simultaneously in order to properly interpret biarticular muscle function. Similar results can be seen for the combined effects of other biarticular hip extensors/knee flexors (Fig. 6); hip joint power reaches its highest value at a crank angle of 134 deg when the knee is producing substantial negative power (Fig. 7). This negative knee joint power is due to the combined effects of BFLH, LG, MG, SM and ST, all of which are producing positive power at the hip and ankle while at the same time producing negative power at the knee. Because our simulations represent maximized muscle power production, this negative power was not the result of poor coordination but rather an inevitable consequence of biarticular muscle function. This type of insight regarding coordinated multi-joint human activity can, within the constraints of current technology, only be obtained through simulations, demonstrating that simulations provide a valuable approach for examining biarticular muscle function.

Our experimental biomechanics data were collected using one cycle crank length (170 mm) and one cycle frequency (120 rev min⁻¹), and thus do not encompass a large parameter space of frequency or muscle excursion amplitude. However, Barratt and colleagues (2011) have previously demonstrated that crank length does not influence joint specific power production during maximal cycling and thus a single, standardized crank length can be used to produce data that are broadly representative of joint power production. In other work, McDaniel and colleagues (2014) have reported that each joint action exhibits an individual power–pedaling rate relationship. However, hip extension, knee extension and knee flexion, the three main power-producing actions, were at or near their maximum at 120 rev min⁻¹. Therefore, we believe that using this single pedaling rate is justified for this investigation.

Our modeling approach has several limitations that must be discussed. First, the muscle parameters provided in OpenSim represent a 50th percentile male whereas our experimental cycling biomechanics data were recorded in a group of trained cyclists. While we scaled muscle and tendon lengths to account for the segment lengths of our subjects, we used the default values for each muscle's cross-sectional area and isometric force as we did not have data necessary to scale these parameters based on the cyclists' anatomy. One example where additional model scaling might have been beneficial is the hip joint power, where our cyclists outperformed the model. It is possible and even likely that these cyclists, as a result of their training, had larger hip extensor muscles than 50th percentile values, and that scaling of those muscles could have improved our model prediction. Second, we prescribed the kinematics to the model and thus our modeling solution is not a true forward solution for maximized power. Rather, we maximized muscle power within the constraints set by the cyclists during maximal cycling. Thus, our approach is similar to that of those who have compared in vivo muscle power with in situ muscle power by experimentally imposing in vivo strain patterns onto muscles during work loops (Askew and Marsh, 2002; Askew et al., 2001; Curtin et al., 2005; Franklin and Johnston, 1997; James and Johnston, 1998; Syme and Shadwick, 2002; Wakeling and Johnston, 1998). Third, we scaled the length of muscle fibers and tendons according to the segment lengths of experimental data subjects prior to obtaining muscle–tendon length trajectories in OpenSim. This approach indicated that some of the muscles functioned at lengths well below resting length [e.g. sartorius (Sar), iliacus (IL), pectineus (Pect)]. These lengths may not be realistic because muscles are known to adapt in length to chronic activity (Ullrich et al., 2009). In addition, individual maximized muscle power, averaged for the cycle, ranged from 1 to 65 W, with 10 muscles producing less than 3 W. If some of these smaller muscles did not voluntarily produce maximized power, their contribution may have been too small to substantially influence the summed power at the hip or knee. Consequently, the excellent agreement of simulated and experimental joint power values strongly suggests that humans maximize power during maximal cycling but does not guarantee that power of each and every muscle is necessarily maximized. Finally, the muscle model we used did not include history-dependent effects which are known to influence force production (e.g. McDaniel et al. 2010; Powers et al., 2014). Despite these limitations, our results demonstrate remarkable agreement of modeled maximized muscle power with voluntary maximal cycling.

In summary, we simulated the maximum power that each muscle could produce within the kinematic constraints of human cycling. The combined power of those simulations agreed very well with experimental joint power for the hip and knee, but less well with joint power for the ankle. We interpret these results to support our hypotheses that humans maximize muscular power for complete shortening/lengthening cycles of hip and knee muscles, but that ankle muscles must act primarily to transfer power from the ankle to the pedal. Thus, for muscles spanning the hip and knee, humans join fish, birds and lizards in their ability to maximize muscular power. Additionally, our simulations provide novel insight into the disparate joint power values produced by biarticular muscles at their proximal and distal joints, where individual joint power can appear to be much greater than actual muscle power. Future applications for this simulation technique may include predicting maximal capability of humans in various clinical and exercise scenarios such as traumatic muscle damage, amputation, tendon transfer surgery, peripheral muscle fatigue and adaptations to training.