Muscle force interacts with stature to influence functionally related polar second moments of area in the lower limb among adult women

Abstract Objectives We sought to determine the relationships between muscle size, function, and polar second moments of area (J) at the midshaft femur, proximal tibia, and midshaft tibia. Materials and Methods We used peripheral quantitative computed tomography to quantify right femoral and tibial J and soft tissue cross‐sectional areas, and force plate mechanography to quantify peak power output and maximum force of the right limb, among athletic women and control subjects. Results Lower limb bone J exhibited strong relationships with estimated force but not power between both groups. Among controls, the strongest relationships between force and J were found at the midshaft femur. Among athletes, these relationships shifted to the tibia, regardless of body size, likely reflecting functional strain related to the major knee extensors and ankle plantarflexors. Together, muscle force and stature explained as much as 82 and 48% of the variance in lower limb bone J among controls and athletes, respectively. Discussion Results highlight the importance of considering relevant muscle function variables (e.g., force and lever arm lengths) when interpreting behavioral signatures from skeletal remains. Future work to improve the estimation of muscle force from skeletal remains, and incorporate it with lever arm length into analyses, is warranted. Results also suggest that, in doing so, functional relationships between a given section location and musculature should be considered.

and known behavior (Macintosh, Pinhasi, & Stock, 2017;Macintosh & Stock, 2019;Niinimäki et al., 2019;Shaw & Stock, 2009a;Shaw & Stock, 2009b). Observed differences in lower limb bone J relative to known and inferred mobility derive from the responsiveness of bone surfaces to prevailing strain conditions (Garn, 1972;Gosman, Stout, & Larsen, 2011). Known as functional adaptation (Ruff, Holt, & Trinkaus, 2006), this process involves bone deposition and resorption relative to localized strain (Ruff, 2008;Seeman, 2008) and can alter the cross-sectional distribution of bone tissue. The most mechanically relevant forces acting on the lower limbs, bending and torsion (Ruff & Hayes, 1983), produce stresses about a neutral cross-sectional axis (bending) or centroid (torsion), whose magnitudes are zero at these neutral locations and increase proportionally as distance from them increases (Frankel & Nordin, 1980;Ruff & Hayes, 1983). This means that the addition of bone to the periosteal surface maximizes diaphyseal J, as it is bone fibers in the periosteal region that must resist the majority of stresses during loading (Jepsen, 2009). By quantifying variation in limb bone diaphyseal J within and among populations in the past, bioarchaeological analyses can infer relative amounts of bending/torsion during life based on the bone's functional response, particularly at the periosteal surface.

| Behavioral signals in lower limb bone crosssectional geometry
Functional adaptation is influenced by a variety of loading-related parameters such as weight-bearing, body size and shape, ground reaction forces, muscle contractions, and the magnitude and frequency of loading (Niinimäki et al., 2019;Robling, 2009). Loading also modulates the normal sequence of bone deposition and resorption that is ultimately under biochemical regulation, involving reproductive hormones, growth hormones, and others (Gosman et al., 2011;Libanati, Baylink, Lois-Wenzel, Srinvasan, & Mohan, 1999;Parfitt, 2002;Tanner, 1989). This means that variation in hormonal status, as well as in the dietary-derived nutrients required for normal hormone secretion (Garn, Guzmán, & Wagner, 1969;Garn, Rohmann, Béhar, Viteri, & Guzmán, 1964;Kenney, McCoy, & Williams, 1994), can alter crosssectional bone distribution in ways unrelated to loading. Because many of these factors are simply unknown among past individuals, inferences about behavior based on bioarchaeological analyses of skeletal variation tend to remain quite broad. When inferred behavior differs distinctly, for example among foragers relative to farmers or industrialized populations, the former tend to exhibit higher sizestandardized mean midshaft femoral and tibial J, attributed to relatively higher terrestrial mobility (Bridges, Blitz, & Solano, 2000;Ruff et al., 1984). These inferences are supported by consistent and similar differences among living athletes relative to controls, where behavior is known (Macintosh & Stock, 2019;Niinimäki et al., 2019;Shaw & Stock, 2009a;Shaw & Stock, 2009b). However, our ability to infer more nuanced detail about the loading environment that has shaped patterns in limb bone structural variation in the past is limited by an inability to reconstruct more specifically the parameters exerting force on the bones during locomotor behaviors.
One of these parameters is muscular contraction, which accounts for a substantial amount of the force applied to bone in normal daily movement (Burr, 1997;Robling, 2009). So closely related are muscle and bone that the biomechanical relationship between the two is typically viewed as a functional "muscle-bone" unit (Fricke & Schoenau, 2007;Schoenau & Frost, 2002). As the lower limb impacts the ground and pushes back off of it during locomotor movements, ground reaction forces (GRFs) are exerted on the limb bones, reflecting the impact loading produced by accelerating and decelerating the body as it contacts the ground (Vainionpää et al., 2007). More specifically, GRFs reflect the activity of the skeletal muscles as they contract eccentrically during the braking phase to slow the body's dropping center of mass and contract concentrically to raise it again during the push-off phase (Ishikawa & Komi, 2008). As a result, variation in activities performed during life can induce substantial variation in the loading to which bones must functionally adapt.
The performance of particularly high-impact activities can exert substantially higher GRFs on the lower limb bones during ground contact than is typical of normal daily movement. For example, the loading associated with some sports, such as triple jump, can exert peak impact GRFs of over 15x body weight (Perttunen, Kyröläinen, Komi, & Heinonen, 2000). These high impacts have been associated with substantial bone functional adaptation in the lower limb, leading to increases of over 20% in midshaft tibial bone strength indices (BSIs) among triple jump athletes relative to controls (Heinonen, Sievänen, Kyröläinen, Perttunen, & Kannus, 2001). In contrast, sport-specific loading regimes involving lower-impact activities, for example power-lifting or rowing, may not necessarily provide sufficient loading stimulus to enhance limb bone CSG properties among athletes relative to controls, despite substantial muscular effort (Macintosh & Stock, 2019;Niinimäki et al., 2017Niinimäki et al., , 2019Nikander, Sievänen, Uusi-Rasi, Heinonen, & Kannus, 2006). In some instances, muscle activity explains more variance in bone strength parameters than sport-related impact loading. For example, among 133 premenopausal female national-level athletes and controls, estimated muscle performance-related joint moments explained 42% more of the variance in BSI at the tibial midshaft than did sportrelated impact loading modality (Nikander et al., 2006). This is supported by Schipilow, Macdonald, Liphardt, Kan, and Boyd (2013), who demonstrated that muscle strength accounted for 13% more of the variance in distal tibial failure load (N), as quantified from finite element analysis, than did sport-related impact loading relative to controls. In the proximal femur, when the effect of gluteal muscle size, strength, and body size are removed, sport-related ground impact does still exert a significant effect at midshaft among females, but only in the presence of very high ground impacts performed at a high rate (Niinimäki et al., 2019). Given this complexity, interpreting the loading environment experienced by limb bones during life from skeletal remains alone can be challenging, and would benefit from the development of more rigorous methods of estimating muscle activity in the past.
One means of estimating local muscle force among living humans is by quantifying muscle cross-sectional areas (MCSAs), as these are roughly proportional to the force that muscle can produce among both males and females (Maughan, Watson, & Weir, 1983b). MCSAs have demonstrated significant explanatory power for cross-sectional geometry, bone areas, and bone strength throughout the lifespan (Heinonen, Mckay, Whittall, Forster, & Khan, 2001;Macdonald, Kontulainen, Petit, Janssen, & McKay, 2006;Rantalainen, Nikander, Kukuljan, & Daly, 2013). However, the estimation of MCSAs from skeletal remains alone has proven difficult, and existing attempts to do so have had mixed results. Though Slizewski, Onau, Shaw, and Harvati (2013) found that radial cortical bone area predicted forearm MCSA among living humans with percent SEs of between 11 and 15%, accuracy varied highly by age and sex. Among young males, Shaw (2010) found that muscle areas were generally poor predictors of tibial cross-sectional bone area and J, with relationships largely reflecting covariation with body size.
These studies assessed relationships between MCSAs and the immediately adjacent cortical bone obtained from the same section location through the limb. However, cortical bone parameters may be more influenced by functionally-related MCSAs located elsewhere in the limb, rather than immediately adjacent to the crosssection of interest. For example, tibial midshaft BSI is more strongly associated with muscle area in the thigh than the calf among older men (>65 years of age) (Rantalainen et al., 2013), while the relative influence of gluteal muscle loading on proximal femoral geometry has also been shown to vary by cross-sectional location (Niinimäki et al., 2019). As a result, the estimation of muscle from bone in the past may prove more successful if we can better account for functional muscle-bone relationships within the limb. An improved understanding of these relationships may also help explain some of the regional variation in CSG properties within the lower limb that is consistently documented among archaeological and living populations.
Tibial diaphyseal structural properties often appear to reflect inferred loading more clearly than femoral diaphyseal properties (Davies & Stock, 2014;Macintosh, Pinhasi, & Stock, 2014;Macintosh & Stock, 2019;Stock, 2006), which is currently attributed to a wider range of influences, such as body size and breadth, acting on the femoral diaphysis (Davies & Stock, 2014). However, the contribution of functional relationships between tibial CSG and the powerful knee extensors of the proximal limb segment could hold explanatory power for these regional patterns in the limb, at least in part.
Though a muscle's cross-sectional area is a good proxy for its ability to produce force, muscle strength per unit of cross-sectional area exhibits considerable variation among subjects, attributable to variation in fiber type composition (Kanehisa, Ikegawa, & Fukunaga, 1994;Maughan et al., 1983b;Maughan, Watson, & Weir, 1983a). Thus, more direct measures of functional strength may be more accurate.
Vertical GRFs can be measured directly from subjects while jumping on a force plate, and these GRFs can be used to estimate maximum muscle force and power produced during whole-limb movements (Anliker, Rawer, Boutellier, & Toigo, 2011;Binkley & Specker, 2008;Rantalainen, Heinonen, Komi, & Linnamo, 2008). The highest peak voluntary GRFs on the force plate are produced by single-leg hopping (Veilleux & Rauch, 2010), providing a good estimate of maximum voluntary muscle force (Fmax) of the lower limb ankle plantarflexors (Anliker et al., 2011). Lower limb Fmax derived from single-leg hopping has demonstrated stronger relationships with tibial bone mass than does calf MCSA in both sexes (Anliker et al., 2011), and has been associated with increased periosteal expansion at the midshaft tibia (Hardcastle et al., 2014). As maximum and polar second moments of area are maximized by periosteal expansion (Ruff, 2008;Ruff & Hayes, 1983), we might expect them to exhibit a strong relationship with force plate-derived Fmax. In support of this, maximum second moments of area (Imax; mm 2 ) at the midshaft tibia have demonstrated significant positive correlations with maximum voluntary force among young adult women (Rantalainen et al., 2008).
Mechanography also allows for the estimation of peak power output (Pmax; W), derived from force estimates combined with velocity during the upward phase of a bilateral countermovement jump. This peak power output is thought to reflect primarily the actions of the hip and knee extensors in addition to those of the ankle plantarflexors (Hardcastle et al., 2014), and Pmax has shown significant associations with tibial compressive and bending strength indices among both sexes (Binkley & Specker, 2008;Hardcastle et al., 2014;Rantalainen et al., 2010). However, of relevance to maximizing polar second moments of area, the relationship between Pmax and cortical bone strength may be reflecting reduced endosteal resorption rather than enhanced periosteal expansion (Hardcastle et al., 2014). This suggests that Pmax may not exhibit a strong relationship with J specifically, though this remains to be determined. Currently, the extent to which limb bone CSG properties reflect whole-limb muscle force or power among past populations is unknown.
1.2 | The functional interaction of body size, composition, and cross-sectional geometry Given the close relationship between muscle and bone, bone mass and total body lean mass are strongly correlated throughout life, regardless of sex or reproductive status (Capozza, Cointry, Cure-Ramírez, Ferretti, & Cure-Cure, 2004;Ferretti et al., 1998). However, the importance of muscle activity in shaping lower limb bone strength parameters is such that some of the behavioral signature being inferred from skeletal variation, even among otherwise similarly-sized individuals, may be attributable to variation in body composition, as well as age-related changes in it. Body size variation is known to impact polar second moments of area through its influence on diaphyseal diameters, the amount of weight borne by the lower limb bones, the extent of mechanical stimulation on them during ground impact, the amount of contractile muscle tissue, and the lever arms about which that muscle tissue exerts bending and joint moments during movement (Niinimäki et al., 2019;Ruff, 2000;Ruff, 2018). As a result, many CSG properties scale with body mass (Parfitt, 1994;Seeman et al., 1996), and polar second moments of area must be standardized to bone length (a proxy for lever arm length) and estimated body mass prior to comparison among individuals and populations (Ruff, 2000;Ruff, 2008;Ruff, 2018). In the lower limb, this enables variation in CSG properties to be attributed more directly to functional adaptation in response to behavior while controlling for the influence of variation in body size (Ruff, 2000).
However, controlling for the effect of lever arm length as a component of body size standardization makes it difficult to explore the contribution of muscle activity to the behavioral signal that we are detecting from skeletal remains. Our potential to explain variation in midshaft diaphyseal bone strength parameters would be improved if we could better estimate muscle force from bone, and incorporate force estimates alongside bone length into statistical analyses. For example, estimated joint moments alone (proxied from stature and maximum muscle force) explain almost 50% of the variance in midshaft tibial BSI among young female athletes and controls, and 42% more variance than ground impact loading (Nikander et al., 2006).
Further, maximum torques directly measured during knee extension are highly correlated with MCSA of the active muscle group multiplied by stature as a proxy for lever arm length among both males and females (Schantz, Randall-Fox, Hutchison, Tydén, & Astrand, 1983). At present, we do not know how much of the variation in lower limb J, within or between limbs, is explained by muscle activity or by variation in relative lean mass. To do so requires the development of better methods of estimating, and then incorporating, relevant muscle function variables (e.g., force and lever arm lengths) into our statistical comparisons. The first step to doing so is building a more accurate picture of the relationships between muscle size, function, and commonly used CSG properties and section locations throughout the weight-bearing lower limb as a whole.

| Objectives
In the current study, we seek to determine what relationship estimated lower limb muscle force and power (proxied from thigh and calf MCSA and whole-limb force plate mechanography) have with polar second moments of area quantified at commonly used section locations in bioarchaeology: the midshaft femur, proximal tibia, and midshaft tibia. We assess these relationships among a group of female athletes and control subjects from whom several key relationships have already been demonstrated. First, significant enhancement of midshaft femoral and tibial J relative to controls was only documented among sports involving ground impact loading (soccer and endurance running); high muscle magnitudes alone (rowing) were not sufficient to elicit significant change in midshaft lower limb bone J relative to controls (Macintosh & Stock, 2019). Second, midshaft femoral and tibial J were most strongly correlated with total body lean mass (r ≤ .72), followed by body mass (r ≤ .60), and fat mass (r ≤ .27), as quantified from bioelectrical impedance (Pomeroy, Macintosh, Wells, Cole, & Stock, 2018). Third, midshaft tibial J was a stronger predictor of total body lean mass than of body mass (adjusted R 2 of 0.51 vs. 0.36, respectively), particularly when bone length was included in regression models (adjusted R 2 of 0.6 vs. 0.42, respectively) (Pomeroy et al., 2018).
Among this mixed sample of females, we assess relationships between muscle size, function, and cross-sectional J among two groups split by physical activity history: recreationally-active control subjects with no history of sport-related impact loading or intensive muscular activity, and active athletes with significant current and past histories of competitive sport. Among the athlete group, sport-related loading characteristics varied (see Section 2) but all were combined into one group in order to more accurately reflect the variable behaviors and activity histories of individuals from fossil, prehistoric, or archaeological populations. Our three main research questions, and associated hypotheses, are as follows: 1 To what extent do estimates of regional and/or whole-limb muscle force and power influence femoral and tibial polar second moments of area?
We hypothesize that estimates of maximum force and peak power output will be significantly predictive of polar second moments of area at all section locations, but that peak power output may be less strongly related to lower limb bone J than maximum force.
2 How does the relationship between J and MCSAs pattern regionally within the lower limb?
Reflecting functional relationships, we hypothesize that proximal and midshaft tibial J will be more strongly correlated with thigh MCSA than with calf MCSA and that midshaft femoral J will be more strongly correlated with thigh MCSA than with calf MCSA.
We further hypothesize that overall relationships between femoral J and MCSAs may be weaker than is the case with tibial J, due to our inability to account for the influence of functionally-related gluteal musculature on the midshaft femur.
3 Are estimated force (MCSAs and Fmax) and/or power variables, in combination with a proxy for lever arm length (stature), better independent predictors of variance in femoral and/or tibial J than body mass and stature alone?
We hypothesize that regression models including a proxied measure of muscle force or power will explain more variance in J than models incorporating only stature and body mass.

| Participants
A total of 102 women took part in the study, and all provided written informed consent prior to participation. All participants were healthy, pre-menopausal women aged 19-34 years, primarily of European descent. They included a range of current and past physical activity levels, from recreationally-active women through to women with elite sporting histories. Participants had neither recent soft tissue injuries that affected muscle performance nor medical conditions or medication use known to affect bone metabolism. All participants completed a questionnaire assessing their past and current recreational and competitive physical activity histories, and a variety of nutritional and hormonal factors known to affect bone, including age at menarche, family history of osteoporosis, active avoidance of dairy products, past and current estrogen-based hormonal contraceptive use, age at first use of estrogen-based hormonal contraceptives, and any history of menstrual irregularity, defined as any history of amenorrhea or oligomenorrhea (fewer than 10 periods in 12 months) not clearly related to contraceptive use (e.g., intrauterine devices).
The control subjects (N = 34) had never taken part in competitive sport and had no current or past participation of more than 3 hr per week of intensive weight-bearing sport, but some were active at a recreational level in yoga, cycling, pilates, swimming, and other nonweight-bearing exercise. The athlete group (N = 68) consisted of women with current and previous participation in sport as high as national-team level, often participating at a competitive level in multiple different sports throughout their lives. Athletes were recruited via their current participation in one of three sports (rowing, endurance running, soccer), but for the purposes of this study were not screened by previous or concurrent loading history, so represent women with a range of impact loading histories. The sports most frequently reported in addition to rowing, endurance running (middle-distance through ultramarathon) and soccer were: triathlon, swimming, field hockey, rugby, cricket, netball, basketball, volleyball, tennis, lacrosse, gymnastics, and equestrian. No athlete had any significant injury within the past year that rendered them inactive for over 1 month.

| Anthropometry
Stature was measured to the nearest 0.01 cm using a SECA 274 stadiometer, and was used as a proxy for lever arm length for the purposes of this study. Body mass was recorded to the nearest 0.1 kg with a SECA electronic scale. Femoral length and maximum tibial length were obtained from participants using sliding calipers according to the methods in International Standards for Anthropometric Assessment (International Society for the Advancement of Kinanthropometry, 2001). For the femur, maximum length was measured as the distance between the proximal border of the greater trochanter and the distal margin of the lateral condyle, and does not reflect the true maximum length of the bone. For the tibia, maximum length was measured as the distance between the proximal medial border of the tibial plateau and the distal margin of the medial malleolus.

| Bone, muscle, and fat cross-sectional parameters
All cross-sectional bone images were collected using peripheral quantitative computed tomography (pQCT; XCT-3000; Stratec Medizintechnik GmbH, Pforzheim, Germany) at the University of Cambridge. Crosssections were obtained from the midshaft femur and tibia, at 50% of maximum length from the distal end, and from the proximal tibia, at 66% of maximum length from the distal end (see Figure 1). However, because it was not possible to obtain true maximum length from the femur among living participants, the "midshaft" femoral slice used here is slightly distal to true midshaft, though is still taken close to diaphyseal transverse cross-sectional minima. Muscle cross-sectional areas were used as proxies for localized muscle force that are not reliant on participant skill or motivation (Frank, Lorbergs, Chilibeck, Farthing, & Kontulainen, 2010). Both muscle and subcutaneous fat cross-sectional areas (MCSA and FCSA; mm 2 ) were quantified at the midshaft femur (thigh MCSA and FCSA) and proximal tibial scans (calf MCSA and FCSA), using Macro analyses in the pQCT manufacturer software (XCT, version 6.2.0), with Contour Mode 1 (−100 mg/cm 3 ), Peel Mode 2 (40 mg/cm 3 ), Cort mode 1 (710 mg/cm 3 ), and an F03F05F05 muscle filter applied. Cross-sectional pQCT images were then imported into ImageJ (http://rsbweb.nih.gov/ij/). The polar second moment of area (J; mm 4 ), a measure of twice average bending and torsional rigidity, was quantified using the bone image analysis plugin BoneJ (Doube et al., 2010), with the "Optimize Threshold" function used to isolate cortical bone from surrounding tissue.

| Whole-limb performance parameters
Whole lower limb maximum force and power were evaluated using a Leonardo Mechanograph Ground Reaction Force Platform (NOVOTEC Medical GmbH, Pforzheim, Germany) and manufacturer's software F I G U R E 1 pQCT slices at the two section locations from which soft tissue parameters were quantified, demonstrating the range of cross-sectional variation in muscle and subcutaneous fat in the sample (Leonardo Mechanography version 4.4). This force plate is split into two sections that enable measurement of the left and right legs separately. The Leonardo Mechanograph system measures force applied to the plate over time, but also measures acceleration, allowing it to quantify vertical velocity and ultimately power output from these force and velocity measurements (Dionyssiotis, Galanos, Michas, Trovas, & Lyritis, 2009). The force platform was calibrated prior to each assessment, and subjects performed all jumps in socked or bare feet.
Maximum voluntary force (Fmax; kN) of the right leg was obtained from the lift-off phase during multiple single-leg hopping (M1LH) on the ball of the foot. Force was acquired from single-leg hopping because this movement pattern produces the highest peak GRFs on the force plate (Veilleux & Rauch, 2010), so is the best proxy for maximum voluntary force exerted by the ankle plantarflexors (Anliker et al., 2011). Ground reaction force curves during multiple single-leg hopping of an athlete and control subject, each most representative of their group's mean Fmax, are presented in Figure 2. The multiple single-leg hopping procedure begins with the subject in an upright standing position, after which they lift the non-testing leg off the plate, and begin hopping on the testing leg while keeping the knee stiff. Subjects were instructed to stand still on both legs until hearing a beep from the computer, at which point they were to perform six to eight consecutive single-leg hops on the ball of the foot, keeping their heel off of the plate at all times, with the objective of minimizing contact time with the plate. Upon completion of each set of hops, the subjects were asked to stand still on both legs again, until hearing two beeps, indicating completion of the test.
Each subject was then allowed to warm up and perform approximately 10-15 submaximal hops to practice. For testing, each subject performed three sets of six to ten consecutive single-leg hops on the ball of the foot. Any hop in which the heel contacted the force plate was excluded; of remaining hops, the one producing the highest Fmax relative to their body weight (indicated as the detected maximum in Figure 2) was identified, and the raw Fmax used for analyses.
Peak power output (Pmax; kW) was obtained for the right leg during the lift-off phase of a single two-legged countermovement jump (S1LJ). The countermovement jump provides insight into an individual's neuromuscular capabilities while remaining resistant to participant fatigue and requiring minimal habituation prior to testing (Cormie, McBride, & McCaulley, 2009;McMahon, Suchomel, Lake, & Comfort, 2018). Power, force, and velocity curves during the countermovement jump of an athlete and control subject, each most representative of their group's mean Pmax, are presented in Figure 3. The jump F I G U R E 2 Ground reaction force curves during multiple single-leg hopping (M1LH) of an (a) athlete and (b) control subject, each most representative of their group's mean Fmax. Detected Maximum tracing (in red) is indicative solely of the timing of peak force and does not reflect force amplitude begins with the subject in an upright standing position, followed by an initial lowering of the center of mass as the subject squats (the countermovement), during which energy is stored in their elastic soft tissues, and then the immediate driving of the body vertically off of the force plate (Dionyssiotis et al., 2009;McMahon et al., 2018). Subjects were instructed to stand as still as possible on the plate until hearing a beep from the computer, at which point they were to rapidly squat to their preferred depth and immediately jump vertically, as high and fast as they could, using both legs as well as freely swinging their arms to elevate their head and chest as high as possible above the ground. Each participant was also asked to soften their landing by absorbing some of the impact through flexing their hips, knees, and ankles, and then to stand upright and still until hearing two beeps, signifying completion of that attempt. Each participant was allowed to warm up and familiarize herself with the force plate prior to testing, after which each subject performed three attempts. The jump in which the subject obtained the highest total peak power output (produced by both legs together) relative to their body weight during the lift-off phase was identified, and raw Pmax from the right leg specifically was isolated from the right half of the force plate and used for analyses.

| Statistical analyses
All data were tested for normality by assessing standardized skewness (skew/SE), with values between −3.29 and 3.29 considered normally distributed. Variables that were not normally distributed were log10 transformed prior to analyses. Independent samples t tests were used to assess differences between athletes and control subjects in 14 continuous parameters. Chi-squared or Fisher's Exact tests were used to compare five categorical parameters. Bivariate correlations were used to assess the relationships between limb bone polar second moments of area at each section location and stature, body mass, muscle areas, and performance variables, and partial correlations were used to do so after controlling for stature and body mass. Bivariate correlations were also used to assess the relationship between muscle areas and body size parameters. Age did not significantly correlate with any muscle or bone variables in our sample of participants, so was not included in regression models.
Hierarchical linear regression models were constructed to assess the variance in muscle area (MCSA) and whole-limb performance variables (Pmax, Fmax) accounted for in J at different locations in the F I G U R E 3 Power, force, and velocity curves during the countermovement jump (S2LJ) of an (a) athlete and (b) control subject, each most representative of their group's mean Pmax. Each subject's maximum force and velocity are also provided. Units for the y-axis are as follows: power, kW; force, kN; velocity, m/s. Each subject's take-off and landing point are indicated; all Pmax values were obtained during the upward phase of the jump, just prior to take-off. Maximum power during take-off was obtained at 3.43 s in (a), and 3.55 s in (b) limb. Hierarchical models for J at the midshaft femur, proximal tibia, and midshaft tibia were constructed, and each was adjusted for baseline stature and body mass. Subsequently, four models were built for each section location, with either thigh MCSA, calf MCSA, Pmax, or Fmax entered. For each model, multicollinearity was assessed based on tolerance, variance inflation factors, and variance proportions; no regression models violated the assumptions of multicollinearity. This procedure was performed separately among athletes and control subjects. Statistical analyses were performed using SPSS Version 24, with significance set at p < .05.

| RESULTS
Descriptive statistics for all variables, and results of independent samples t tests between athletes and controls, are provided in Table 1.
The athlete group was, on average, significantly taller (p = .018) than the control group and had significantly larger absolute thigh and calf MCSA (all p < .001) but significantly smaller absolute thigh and calf FCSA (p = .001). Despite these differences in body composition, there were no significant differences in body mass between the two groups.
Thus, our two groups were similar in body mass but differed significantly in the proportion of that mass that is capable of generating performance-related bending and joint moments on limb bone diaphyses. Images of pQCT slices at the midshaft femur and proximal tibia demonstrating the range of variation in soft tissue parameters in the sample are provided in Figure 1. Given their substantially larger stature and absolute muscle areas, athletes were able to generate significantly higher raw Fmax and Pmax (p < .001) than controls. Athletes also had significantly higher raw values of J throughout the limb (all p < .001). There were no significant differences between the athlete and control group in any questionnaire-based self-reported dietary or hormonal parameters.

| Body size
The results of bivariate correlations are provided in Table 2   relationship and its patterning throughout the limb differed between athletes and controls. Among control subjects, stature and body mass were more highly correlated with lower limb J (r = .572-.754) than among athletes (r = .524-.580), and these relationships were strongest at midshaft locations, particularly in the femur (stature: r = .754; body mass: r = .648; all p < .001), and weakest at the proximal tibia (stature: r = .629; body mass: r = .612; all p < .001). In contrast, among athletes, relationships between body size variables and limb bone J were consistently strongest at the proximal tibia (stature: r = .580; body mass: r = .575; all p < .001), and weaker at the midshaft locations.

| Soft tissue areas
The extent of the relationship between body size variables, MCSA, and FCSA is assessed in Table 3 The results of partial correlations assessing relationships between muscle areas and limb bone J after controlling for stature and body mass are provided in Table 2 and depicted in Figure 5.
Among controls, only thigh MCSA remained significantly positively correlated with J, with relationships highest at the midshaft femur (r = .570, p < .01) and lowest at the proximal tibia (r = .368, p < .05); FCSAs no longer exhibited any significant relationship with limb bone J. Among athletes, calf MCSA remained significantly positively T A B L E 2 Bivariate and partial correlations assessing relationships between regional lower limb bone J and both regional soft tissue areas and dynamic performance variables (N = 102)  F I G U R E 4 Bivariate correlations between regional polar second moments of area (J) and both regional soft tissue areas and dynamic performance variables Pearson's correlations assessing relationships between regional MCSAs and body size variables

| Explanatory power of MCSAs and whole-limb performance variables: Control subjects
Among control subjects, hierarchical linear regression models revealed that baseline adjustment for body size (stature and body mass) explained 70.6% of the variance in midshaft femoral J (Table 4), 52.4% of the variance in proximal tibial J (Table 5), and 49.6% of the variance in midshaft tibial J (Table 6)

| DISCUSSION
This study aimed to determine what relationship estimated lower limb muscle force and power (proxied from thigh and calf MCSA and whole-limb force plate mechanography) have with J quantified at commonly used section locations in bioarchaeology. Polar second moments of area at the midshaft femur, proximal tibia, and midshaft tibia were all strongly correlated with muscle size and function F I G U R E 5 Partial correlations between regional polar second moments of area (J) and both regional soft tissue areas and dynamic performance variables variables. However, these relationships between bone and muscle were patterned differently within the limb of athletes and controls, and reflected body size to differing extents. A summary of the main findings relevant to each research question is found in Table 7.

| Muscle force production is significantly associated with functionally related lower limb bone polar second moments of area
Our research sought to investigate the extent to which estimates of regional and/or whole-limb muscle force and power influence femoral and tibial polar second moments of area, and the regional patterning of these relationships within the lower limb. Results found that, among both female athletes and controls, force estimates (both Fmax and MCSAs) consistently exhibited significant relationships with limb bone J even when controlling for body size. Whole-limb Fmax exhibited stronger relationships with J in the tibia relative to the femur among both athletes and controls, suggesting that at least some of the variance in tibial J across the bioarchaeological literature may be reflecting the functional influence of muscle force-generating capacity. The single-leg hopping from which Fmax was quantified exposes the tibia in particular to maximal forces exerted by the main ankle plantarflexors, as well as the knee extensors to a lesser extent (Rantalainen et al., 2009;Sumnik et al., 2013;Veilleux & Rauch, 2010). The main ankle plantarflexors of the calf run parallel to the tibia, so during normal bipedal locomotion, they exert mainly compressive force on the diaphysis, particularly distally (Rittweger et al., 2000;Wehner, Claes, & Simon, 2009). However, as compressive load interacts with tibial diaphyseal curvature in more proximal section locations (Bertram & Biewener, 1988;Biewener, 1991;Capozza et al., 2010;Garcia & da Silva, 2004;Wehner et al., 2009), bending moments on the tibia increase. This may contribute to Fmax being a stronger predictor of J at the proximal rather than midshaft tibia.
The interpretation of lower limb bone J clearly requires the consideration of functional relationships with musculature. Among athletes, calf MCSA consistently maintained stronger relationships with J more distally, at the midshaft location (unadjusted r = .442; adjusted r = .423; p < .01), where it explained 4% more of the variance in J than it did at the proximal tibia (12% vs. 8.1%), after accounting for stature. Among control subjects, thigh MCSA was a significant independent predictor of J more distally, at both the proximal and midshaft tibial J, explaining an additional 6% and 6.6% its variance, T A B L E 4 Linear regression modeling investigating the explanatory power of muscle areas and dynamic performance variables for midshaft femoral J (N = 102) respectively, while calf MCSA was not. This pattern is similar to that documented by Rantalainen et al. (2013) among older adult males, whereby density-weighted polar cross-sectional moment of inertia (mg/cm) at the tibial midshaft were better predicted by thigh MCSA than calf MCSA (Rantalainen et al., 2013). Results of the current study may explain why earlier attempts to estimate muscle area from skeletal remains have had poor results: limb bone J is not necessarily strongly correlated with adjacent muscle areas from the same section location, but rather functionally related areas more proximally in the limb.
Despite being quantified in part from force measurements, Pmax consistently exhibited much weaker relationships with lower limb bone J than did Fmax, and never significantly explained any of the var- tories, may also be adding an additional source of variation that is known to affect peak power output: muscle conditioning (Hardcastle et al., 2014). For example, strength training has been shown to influence peak power output in the countermovement jump (Cormie et al., 2009). Certain types of training in particular, for example plyometrics (Holcolm, Lander, Rutland, & Wilson, 1996) or resistance training (McBride, Triplett-McBride, Davie, & Newton, 2002), improve the capability to generate not just force but also acceleration, both crucial to maximizing power output. The variation in athletic history of the participants in this study introduces variation in experience and training that may be influencing their ability to generate power. In contrast, a variety of factors influence Fmax production, including kinematics, elastic properties of the Achilles tendon, and the biomechanical properties of the ankle and foot joints, all of which are less influenced by training and are relatively consistent across the lifespan (Hardcastle et al., 2014).
Interestingly, results of the current study also found that regional FCSAs were not related to lower limb bone J among control subjects, but they were among athletes, where both thigh and calf MCSAs were significantly negatively correlated with J throughout the lower limb.
This likely reflects variation in the magnitude and duration of sport participation and the predominant ground impact loading characteristics among the athlete sample, rather than any direct functional T A B L E 5 Linear regression modeling investigating the explanatory power of muscle areas and dynamic performance variables for proximal tibial J (N = 102) relationship between fat and bone. Athletes who had participated intensively for many years in sport that exerts high ground impact loading on the lower limb have high polar second moments of area (Macintosh & Stock, 2019), and likely also have low fat mass relative to athletes with less intensive sporting histories.

| Considering the interaction between muscle force and stature when interpreting cross-sectional geometry
Our research also sought to investigate whether or not estimated force (MCSAs and Fmax) and/or power variables, in combination with a proxy for lever arm length (in this case, stature) are better independent predictors of variance in femoral and/or tibial J than stature and body mass alone. Stature and body mass together explained a substantial portion of the variance in J among control subjects, decreasing progressively from 71% at the midshaft femur to 53% at the midshaft tibia. This pattern supports bioarchaeological work documenting stronger correspondence between body size/breadth and lower limb bone CSG properties in the femur than tibia. However, this pattern was not documented among female athletes, where stature and body mass together exerted their strongest influence on J at the proximal tibia, explaining 37% of its variance. The major knee extensors and ankle plantarflexors make up a large component of body mass, and that musculature is functionally related to proximal tibial CSG. Thus, it may be that body mass itself is not particularly influential for shaping J among athletic women, but rather the muscular component of that mass (Pomeroy et al., 2018).

| Implications for behavioral reconstruction
There are many components shaping the behavioral signal that we attempt to infer from variation in lower limb bone structural properties. The two groups of women in our study, athletes and control subjects, did not differ in bone-related nutritional or hormonal status variables derived from questionnaire data, including age at menarche, family history of osteoporosis, dairy avoidance, past or current hormonal contraceptive use, age at first hormonal contraceptive use, and history of menstrual irregularity. Further, both athletes and controls were relatively homogeneous in age and ethnicity, none had a self-reported history of eating disorder, all were premenopausal, and none had had children. The major differences between these two groups of women were related to loading characteristics, predominantly ground impact loading and muscle activity. In the absence of any significant lifetime ground impact loading, as was the case with the control group, from 63.7 to 81.6% of the variance in femoral and tibial J was accounted for just by a proxy for lever arm length (stature) and muscle force (Fmax), while body mass was not significant.
This result suggests that variation in diaphyseal bone strength parameters may be better explained by the incorporation of a proxy for lever arm length (bone length in the case of skeletal remains) and improved estimates of muscle force from bone. Even among women with significant histories of impact loading associated with changes in their midshaft femoral and tibial J (Macintosh & Stock, 2019), 39.8 to 41.2% of the variance in femoral and tibial J was still accounted for just by a proxy for lever arm length (stature) and muscle force (Fmax). Results also demonstrate that the estimation of muscle force from limb bone cross-sections should consider functional relationships between muscle and bone. For example, though midshaft femoral J exhibited stronger correlations with thigh MCSA obtained from the same section location than with less functionally related calf musculature distally in the limb, both proximal and midshaft tibial J exhibited stronger correlations with MCSAs from more proximal section locations. Further, though the midshaft tibial section location is utilized more commonly than the proximal 66% section location in behavioral inferences, the latter exhibits particularly strong relationships between muscle and bone among active individuals.

| CONCLUSION
Lower limb bone polar second moments of area exhibited strong relationships with muscle force production, both estimated from mechanography and proxied by muscle areas, but not with power output. Muscle force, in combination with a proxy for lever arm length (in this case, stature), explained as much as 82% of the variance in lower limb bone J, and muscle force consistently replaced body mass as an independent predictor of variance in every instance. Results highlight the potential of better methods of estimating relevant muscle function variables (e.g., force and lever arm lengths) and their incorporation into statistical analyses for understanding behavioral signals inferred from skeletal remains. Potential functional relationships between musculature and limb bone diaphyseal J were documented, with particularly strong relationships between muscle force estimates and tibial J among athletes, and femoral J among controls.
Future work on the interaction of musculature, size/proportions, and impact loading as interacting mechanisms shaping bone functional adaptation across a variety of contexts and locations in both the upper and lower limbs of both sexes is needed.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available on request from the corresponding author. The data are not publicly available due to privacy or ethical restrictions.