Estimating Lower Limb Skeletal Loading

Osteoporosis, accidents and subsequent bone fractures cause suffering on an individual level, as well as an economical burden to the society (Ortiz-Luna et al., 2009; Stevens & Olson, 2000). It has been estimated that, in Finland alone, between 30,000 to 40,000 osteoporosis-related fractures occur annually and that 400,000 Finnish people have osteoporosis (Duodecim, 2008). There are a few potential ways of preventing bone fracture, i.e. strengthening bones and/or preventing falls (Ortiz-Luna et al., 2009; Stevens & Olson, 2000). In order to withstand prevalent loading without breaking; while remaining relatively light in weight to allow for locomotion, bones have the ability to adapt their structure to functional loading (Frost, 2000; 2003; Sievanen, 2005). It has been demonstrated that physical activity affects the weight bearing skeleton more than the non-weight bearing one (Mikkola et al., 2008), and it may therefore be argued that, the skeleton is loaded mainly by locomotory actions that impart strains on bones. Bones are loaded in daily activities by muscles accelerating and decelerating body segments and resisting the pull of gravity (Burr et al., 1996). Since falling is the single most significant bone fracture risk factor (Jarvinen et al., 2008) and up to 90% of fractures are caused by falls (Cummings & Melton, 2002; Stevens & Olson, 2000; Wagner et al., 2009), exercise can be viewed as a potential intervention for fracture prevention. Exercise seemingly has a potential of both reducing the fall rate and also increasing bone strength. In agreement, exercise interventions have been shown to successfully decrease the fall rate (Kemmler et al., 2010; Korpelainen et al., 2006), to strengthen the bones and to decrease the fracture rate (Korpelainen et al., 2006; Sinaki et al., 2002). It is quite obvious, that in some cases even the strongest bone could not withstand the loading associated with falling. Therefore, it is reasonable to question whether strengthening bones makes a significant contribution in preventing fractures. A prospective study has shown that people with higher calcaneal bone mineral density had lower fracture rate even while the fall rates between the groups were relatively similar (Cheng et al., 1997). The differences between the groups in bone mineral density were relatively large (~10%). In agreement, prospective studies have shown that increasing the amount of bone mineral with drugs by about 10% is effective in decreasing fracture rates (Cummings et al., 2009). On the other hand, the increments in bone mineral amount associated with year long exercise interventions are relatively modest, i.e. in the order of 1 2% (Nikander et al., 2010). There are examples of fractures for which this kind of relatively modest bone strength increase does play an important role. One of the most convincing examples is vertebral fractures. It has been shown in a prospective study that vertebral bone mineral density (represents 11


Introduction
Osteoporosis, accidents and subsequent bone fractures cause suffering on an individual level, as well as an economical burden to the society (Ortiz-Luna et al., 2009;Stevens & Olson, 2000). It has been estimated that, in Finland alone, between 30,000 to 40,000 osteoporosis-related fractures occur annually and that 400,000 Finnish people have osteoporosis (Duodecim, 2008).
There are a few potential ways of preventing bone fracture, i.e. strengthening bones and/or preventing falls (Ortiz-Luna et al., 2009;Stevens & Olson, 2000). In order to withstand prevalent loading without breaking; while remaining relatively light in weight to allow for locomotion, bones have the ability to adapt their structure to functional loading (Frost, 2000;2003;Sievänen, 2005). It has been demonstrated that physical activity affects the weight bearing skeleton more than the non-weight bearing one , and it may therefore be argued that, the skeleton is loaded mainly by locomotory actions that impart strains on bones. Bones are loaded in daily activities by muscles accelerating and decelerating body segments and resisting the pull of gravity (Burr et al., 1996). Since falling is the single most significant bone fracture risk factor (Järvinen et al., 2008) and up to 90% of fractures are caused by falls (Cummings & Melton, 2002;Stevens & Olson, 2000;Wagner et al., 2009), exercise can be viewed as a potential intervention for fracture prevention. Exercise seemingly has a potential of both reducing the fall rate and also increasing bone strength. In agreement, exercise interventions have been shown to successfully decrease the fall rate (Kemmler et al., 2010;Korpelainen et al., 2006), to strengthen the bones and to decrease the fracture rate Sinaki et al., 2002). It is quite obvious, that in some cases even the strongest bone could not withstand the loading associated with falling. Therefore, it is reasonable to question whether strengthening bones makes a significant contribution in preventing fractures. A prospective study has shown that people with higher calcaneal bone mineral density had lower fracture rate even while the fall rates between the groups were relatively similar (Cheng et al., 1997). The differences between the groups in bone mineral density were relatively large (~10%). In agreement, prospective studies have shown that increasing the amount of bone mineral with drugs by about 10% is effective in decreasing fracture rates (Cummings et al., 2009). On the other hand, the increments in bone mineral amount associated with year long exercise interventions are relatively modest, i.e. in the order of 1 -2% . There are examples of fractures for which this kind of relatively modest bone strength increase does play an important role. One of the most convincing examples is vertebral fractures. It has been shown in a prospective study that vertebral bone mineral density (represents 11 www.intechopen.com bone strength indirectly) may be increased with back strengthening exercise program and that the higher vertebral bone mineral density is also associated with lower incidence of vertebral fractures (Sinaki et al., 2002). Furthermore, exercise induced bone gains may be more appropriately situated in order to maintain bone strength compared to drug induced bone gains. In case of drug therapy, the bone mineral amount increment is distributed rather evenly through the bone 1 , while physical activity associated bone increments are situated mechanically advantageously (Leppänen et al., 2010;Ma et al., 2009;Rubin & Lanyon, 1984). In pure compressive or pure tensile loading, the geometry of the bone makes no difference to bone strength (excluding slender structures such as the ones found in trabeculae, which are susceptible to buckling). The only important cross-section parameter is its area. In bending and torsion, however, the geometry of the bone makes a big difference. As a rule of thumb, the further away from the bending/rotation axis the material is situated, the more significant contribution it makes to the strength of the bone (further details given under section 2.1) (Fig.  1). Coincidentally, the toughest part of bone is the external shield, which is also subjected to the highest loads. Since exercise interventions are effective and feasible, how does one optimize the osteogenic loading of the intervention? Strains (Fig. 4) are one of the most important stimulants to bone adaptation Lanyon (1987); Turner (1998), and therefore designing osteogenic interventions could benefit from the knowledge of bone strains at different cross-sections in a wide range of exercises. General rules for bone adaptation are relatively well established, i.e.; bones respond increasingly with increasing 1) strain rate, 2) strain magnitude and 3) number of loading cycles (Lanyon, 1987;Turner, 1998;Whalen et al., 1988). In addition, loading direction affects the response, i.e.; unusual loading direction causes a larger response (Rubin & Lanyon, 1984). In bone adaptation, also the interval between loadings plays an important role. Bone responsiveness to loading decreases relatively rapidly, and dividing the loading bout into smaller bouts separated by several hours is more osteogenic than doing the exercise all in one bout (Umemura et al., 2002). If bone strains are known, the osteogenicity of loading may be 1 At the moment of writing this chapter, we were unaware of a published article we could cite to confirm this statement. However, the FREEDOM study (Cummings et al., 2009) did include computed tomography measurements, the results of which were presented at least at the ECCEO-IOF11 congress held March 23-26 2011 in Valencia, Spain. Although not mentioned in the presentation, one of the authors of the present chapter was in the audience and queried the presenter (prof. Genant) afterwards whether the bone geometry had changed as a result of the intervention. As one might have expected for an anti-resorptive drug, the geometry had not changed according to the presenter.  (Turner, 1998;Turner & Robling, 2003)) Eq. 1.
Where ǫ is the magnitude of mechanical strain, δǫ δt denotes the mechanical strain rate in [Hz], i is the number of loading cycles, and t is the interval between loading expressed in hours.
(a) Osteogenicity as a function of repetitions.
(b) Sensitivity as a function of loading interval. Fig. 2. Osteogenicity of loading is linearly related to strain magnitude and rate, while loading cycle dependency and interval between loadings dependency are non-linear.
However, for human applications, bone strains are generally not known. Furthermore, measuring bone strains in vivo is invasive and limited to superficial bone sites (Hoshaw et al., 1997;Lanyon et al., 1975;Milgrom et al., 2000), thus not feasible in practice. Consequently, indirect estimates have to be used in practical applications. Since, intuitively there is an association between ground reaction force and bone strains, it has been suggested that osteogenic index may be calculated from reaction forces (Turner & Robling, 2003;Whalen et al., 1988). Subsequently, reaction force based osteogenic index has been shown to be able to successfully differentiate two exercise regimes producing differing bone responses from each other in a prospective human study (von Stengel et al., 2007;. Although, using reaction force as a substitute for bone strains is tempting there are some shortcomings of this approach which will be demonstrated later under heading 3. There are at least a few examples of osteogenic exercise interventions reported, where bone response has been observed on one bone of the lower limb (and spine), while other lower limb bones have not responded (e.g. Vainionpää et al. showed increase in femoral neck bone strength following a high impact exercise intervention in post menopausal women, while tibial mid-shaft bone strength was unaffected by the same intervention (Vainionpää, 2007;Vainionpää et al., 2005). While at first glance it may appear that the aforementioned is unexpected, the phenomenon may be explained by acknowledging that bone response is site specific. High impact exercises are characterized as high impact by the nature of causing high reaction forces. High reaction forces occur in exercises in which the lower limb is mainly compressed (e.g. above 10 times body weight has been measured under one foot in triple jump 2nd jump eccentric phase (Perttunen et al., 2000)), while the leg is straight. This kind of posture enables high reaction forces, since the inertia of the body is efficiently transmitted through hip, knee and ankle joint without the joints giving. If force is not transmitted through the joint centres, the joints give and reaction forces above 6 times body weight are rarely recorded (Ishikawa et al., 2005) even under two feet. In such situation the kinetic energy of the body is damped by muscles and tendons. Bone loading condition (compression or bending) has probably a major effect on bone strains in long bone mid-shaft (Heinonen et al., 1996). Consider a 41 cm long hollow 245 Estimating Lower Limb Skeletal Loading www.intechopen.com cylinder with an inner radius of 6 mm and an outer radius of 12 mm (corresponds to the dimensions of a rather typical tibia). Should the other end of this cylinder be rigidly fixed and either an axial force (compression) or a perpendicular force (bending moment) with the same magnitude be applied on the other end of the cylinder, the latter would result in 50x greater maximal stress (or strain) at the outer surface of the mid-shaft (furhter details on how stress may be estimated in section 2) (Rantalainen et al., 2010). A force during a high impact exercise, which is compressive for tibia will most likely cause bending in femoral neck and therefore it is plausible that a difference may be observed between femoral neck and tibial mid-shaft bone adaptation (or lack of) following the same high impact exercise intervention.

Estimating lower limb loading with reaction forces
Several indirect methods have been used in estimating skeletal loading ranging from qualitative analyses to energy expenditure measurements. However, as per the title of the book, only mechanical loading estimates will be considered in the present text 2 . One relatively popular method to estimate loading, and thus strains, caused by locomotory actions on the bones is to examine the ground reaction forces or accelerations of the center of mass registered during these actions Jämsä et al. (2006); Turner & Robling (2003); Vainionpää (2007); Vainionpää et al. (2006). The advantage of measuring reaction forces or accelerations of the center of mass is the relative simplicity of the measurement. However, it is rarely brought into question if estimating skeletal loading from ground reaction forces is reasonable. Joint angles as well as muscle activity have a great influence on the loading of different bones and should be considered. In addition, diverse bone geometry and mineral content also have a significant influence on bone strains. An equal amount of force applied to different bodies will lead to different bone strains if the mineral amount and/or geometry of the loaded bones differ. The next evident step is thus to include the moment arms and joint angles into the calculations in inverse dynamics. This step, however, comes with an increase in the complexity of measurements as kinematics (i.e. movements of body parts in space) need to be recorded in addition to the aforementioned reaction forces. Furthermore the point of force application and direction of force needs to also be determined (detailed description of inverse dynamics is outside of the scope of the present chapter and detailed explanations can be found from other biomechanical textbooks, e.g. (Robertson, 2004)). However, even with inverse dynamics, the whole loading environment is not determined, since only external forces have been considered thus far. Muscles cause the movement of body parts in relation to each other and have relatively short moment arms compared to external forces. Therefore, the muscle (i.e. internal) forces are much higher than the measured reaction forces and neglecting these internal forces affects the loading estimate. This last step of estimating the loading environment is by far the most complex and includes arguably the biggest sources of error. Many muscles cross a given joint causing the same movement of the joint (e.g. soleus and gastrocnemii muscles plantar flex the ankle) making calculating the internal forces an overdetermined (redundant) problem, which means that there are infinite number of possible solutions of muscle forces 3 .

Stress and strain
It needs to be mentioned that loading, forces, stresses and strains are used somewhat interchangeably in the present text, while generally, forces, stresses and strains should be differentiated. Therefore, the remainder of the paragraph tries to explain what are the differences between forces, stresses and strains and why they should not be mixed. The relationship between the applied load (=force) and relative deformation (Fig. 4) caused by the load can be plotted as a load-deformation curve. The curve may be divided into two parts, the linear elastic region and the non-linear plastic region. The stiffness of the structure is the slope of the elastic region, yield is the load at transition from the linear to elastic region and toughness is the area under curve up to yield or failure ( Fig. 3) (Currey, 2001;Martin & Burr, 1989;Turner & Burr, 1993). The load-deformation curve may be converted to stress-strain curve by accounting for the geometry of the structure and the loading situation. The strain by definition is the deformation divided by the original dimension. If the loading force is compressive (or tensile), then the stress is simply force divided by the cross-sectional area of the structure. In bending and torsion, the stress varies throughout a cross-section. The slope of the stress-strain curve in turn is the Young's modulus, i.e. stiffness, of the material. The strength of the structure is defined as the load at which the structure either yields (yield strength) or breaks (breaking strength, which equals ultimate strength in bone). As was the case with stiffness, the strength may be reported as a material property (yield or breaking stress) or structural property (yield or breaking load). Because of the anisotropic structure of bone material, Young's modulus and breaking stress are different in longitudinal and transverse directions. Furthermore, bone ultimate stress is different in tension, compression and shear. While looking at different types of loading: bending, compression and shear, it becomes clear that the geometry of bone is of utmost importance in bending and shear loading (Currey, 2001;Turner & Burr, 1993).

Example of how geometry affects bone strength
Bending stress in a bone cross-section in a given loading configuration is determined by the second moment of area (also known as cross-sectional moment of inertia) of the cross-section as shown in (Fig. 1 A). Stress (indicated by the shaded area) is distributed unevenly throughout the cross-section, being by definition zero at the neutral axis and increasing its magnitude linearly while moving away from the neutral axis. The stress (σ) at a given point within the cross-section is determined by multiplying the bending moment (M) with the distance coordinate (y) from the neutral axis (dashed line) and dividing the product with the second moment of area (I) of the cross-section in relation to the neutral axis (i.e. σ = M I y). This leads to the conclusion that the bending stress will be positive (tensile) on one side of the neutral axis and negative (compressive) on the other. In case of symmetrical cross-sections the stress magnitudes on both extremes will be equal. In case of a real bone it is rarely the case. Let's simplify the bone mid-shaft cross-section to a hollow oval with external radii a e , b e and internal a i , b i as depicted in (Fig. 1 B). The area of an oval is equal to πab, where a and b are the radii of the oval and thus the area of an oval ring is given by π(a e b e − a i b i ). If the area is increased by 2% in a way observed in association to exercise interventions, only b e will be increased as shown in figure (Fig. 1 C). Therefore the following relation holds: and can be transformed to compute the corresponding new external radius, b e−new , with the means of following equation: The cross-sectional second moment of area of a solid ellipse,I, is equal to: and thus for an elliptic ring, the second moment of area will be: Lets assume a i = 0.5, b i = 0.5 and a e = 1.0, b e = 1.5, which if the units were in cm, corresponds to young adult men tibial mid-shaft. In this particular case the cross-sectional inertia prior to the intervention would have been 2.60cm 4 and after 2.74cm 4 , which is a 5% increase. Accordingly the change would have lead to a 3% decrease in maximal stress. For further details on bone mechanics, please refer to pertinent literature such as prof. Currey's book (Currey, 2002).

Reaction forces and inverse dynamics in estimating skeletal loading
Lower limb skeletal loading (bone stresses) may be estimated with a combination of kinetic and kinematic measurements (Anderson et al., 1996). Kinetic measurements are typically the measurement of vertical and horizontal ground reaction forces e.g. during locomotion, whereas the kinematic measurements (displacements, velocities and accelerations of segments and respective angular equivalents of joints) are derived from the motion captured from the subject. The measured kinematics is combined with the measured kinetics in inverse dynamics to calculate joint moments. Judging by the apparent similarity of the patterns of ankle moment and tibial mid-shaft maximal principal strain in walking (Fig. 5), simply looking at joint moments appears to give a reasonable estimate of tibial mid-shaft bone strains. However, it needs to be remembered that tibia is also affected by the knee joint moment and the patterns of the knee joint moments differ quite markedly from the patterns observed at the ankle joint. Fig. 4. Strain describes the relative deformation of a structure. If the rectangle (representing a solid structure) on the left is subjected to a tension force as indicated by the arrows on the narrower rectangle (representing the same structure under tension), the rectangle will deform. If the structure is incompressible, the width will change roughly half of the relative change in height but to opposite direction (remember that any given real life structure is three dimensional and volume is maintained if the structure is incompressible, i.e. Poisson's ratio = 0.5), i.e. the rectangle becomes longer but narrower. If the force was compressive, the rectangle would become shorter but thicker. In addition to being subjected to both tensile and compressive strain at the same point, due to being incompressible, the opposite forms of strain also cause shear strain at the same point. Unlike the example given, bone is compressible to some extent (cortical bone Poisson's ratio is somewhere between 0.2 to 0.4). Nevertheless, any given loaded point in bone will nevertheless always develop all forms of strains.
Furthermore, if one were to estimate the magnitude of bone stresses and strains, just calculating based on external forces, (i.e. joint moments), and neglect the internal forces (i.e. muscle forces), the stresses and strains could potentially be greatly miss-characterized (an example is given under section 4). Again, using tibial mid-shaft as an example consider standing on the balls of the feet. If the fact that the moments are caused by muscles were neglected, one would include only ankle joint moment and the weight of the body above tibial mid-shaft in calculating the tibial mid-shaft stresses and strains. However, since the ankle joint moment is in fact produced by the ankle extensors, the compressive and bending forces cause by the muscle activity should be taken into consideration instead of joint moments. Due to the ratio of contact point to ankle joint moment arm and ankle extensor muscles' insertion to ankle joint moment arm, the force required from the ankle extensor musculature is in the order of two to three times body weight (Finni et al., 1998;Ishikawa et al., 2005;Komi et al., 1992). Furthermore, bending moments and compressive forces load the bone in a very different manner and neglecting the muscles as the cause of the joint moment has the potential of greatly influencing the result (Fig. 6). Knowing bone stresses and strains in a single point is useful only to a limited extent and for most practical applications, the stresses and strains need to be estimated across the bone cross-section and at several points through the length of the Fig. 5. Schematic illustration of tibial mid-shaft anteromedial aspect maximal principal strain (Lanyon et al., 1975) (upper pane) and ankle sagittal plane joint moment (Silder et al., 2008) during a walking cycle.
bone. Therefore, it quickly becomes evident that building models for individual applications, solving the muscle forces (which are redundant and therefore need to be estimated with some kind of optimization criteria) and calculating the stresses and strains is not feasible and a more general modelling approach needs to be chosen. Flexible multibody dynamics methods can be used as such a generalized modelling method.

Flexible multibody dynamics
Flexibility can be accounted for in a number of ways in multibody applications. The linear theory of elastodynamics can be considered a traditional approach to account for flexibility. This approach relies on the assumption of small deformations in the flexible bodies. Thus, rigid body simulation is decoupled with the deformation computation. The rigid body simulation is performed to obtain the external, as well as internal, forces acting upon each of the bodies. These forces are later imposed on the finite element model of the body for which deformations, stresses and strains can be obtained (Erdman & Sandor, 1972;Lowen & Chassapis, 1986;Lowen & Jandrasits, 1972;Winfrey, 1971). Standard multibody and finite element solvers can be used for this method, which is a great benefit. Additionally, a considerable increase in computational speed can be achieved if the force application points are known a priori and do not vary over the duration of the simulation. In such cases, the use of linear finite element analysis is reduced to a single computation of the full model. Strains, stresses and deformations can then be computed for each time step as post-processing. Conversely, the flexibility of the bodies does not affect the multibody simulation behaviour, which is the main disadvantage, especially in case of considerable deformations. Fig. 6. Freebody diagrams for estimating tibial mid-shaft bone loading while standing on one leg on the ball of the foot. A) Forces accounted for when estimating tibial loading with just the ground reaction force, i.e. ground reaction force (GRF) and its counter force (weight of the subject). B) Estimating tibial loading after inverse dynamics. Compressive loading is accounted for by the GRF and bending is accounted for by the ankle joint moment. It is quite easy to imagine a situation, in which tibial loading is not proportional to ground reaction forces. Consider the situation in pane A and then imagine one standing with the ankle fully plantar flexed. In both cases the GRFs are the same, while ankle joint moment is changed (if one were a ballerina, the ankle joint moment could actually be made to drop to zero). C) Tibial loading is estimated from joint contact forces comprising reaction forces and muscular forces.
Lumped mass formulation is another method that can be used to describe mechanical flexibility (Hughes, 1979;Hughes & Winget, 1980;Huston, 1981;1991;Kamman & Huston, 2001;Raman-Nair & Baddour, 2003). In this formulation, a flexible body is replaced by a set of point masses connected via springs. Using a sufficient amount of springs and masses allows for reasonably accurate mass distribution, within inhomogeneous bodies as well. Similar to the linear theory of elastodynamics, the lumped mass approach does not require that any changes be made in the standard rigid multibody solver. In contrast to the rigid body representation, the performance decays with the rise of discretization precision. Therefore, this method can only be practically applied to beam structures. The third major method for introducing flexibility into the multibody formulation is the floating frame of reference (Shabana, 2005). Formulation relies on coordinate partitioning so that one set of coordinates is used to describe the flexible body's reference frame in the global coordinate system and another set of coordinates is used to describe the deformation of the body in the local frame of reference. Originally, the deformation of the body was described in a similar fashion to the finite element method, resulting in a remarkable increase in regards to computational effort. However, Shabana and Wehage (Shabana, 1983) have developed a solution to this problem by replacing the full finite element models of flexible structures with deformation modes description. This allowed for a reduction in the amount of deformation coordinates (in the range of thousands) to a reasonable amount of modal coordinates, making this method an effective compromise between accuracy and computational effort.

Modelling results obtained with flexible multibody dynamics in human locomotion
There is a relative paucity of dynamic bone loading modelling results in human locomotion in the present literature. One notable exception is the results from Anderson et al. (1996), who modelled femoral neck maximal stresses in several gym exercises. However, Anderson et al. (1996) reported bone stresses from only one location (i.e. the femoral neck). The floating frame of reference method with modal coordinates has been successfully applied in modelling tibial and femoral bone strains during walking Nazer et al., 2011;Nazer, Klodowski, Rantalainen, Heinonen, Sievänen & Mikkola, 2008;. When looking at femoral neck strain estimates and tibial mid-shaft strain estimates it may be seen that the timing of peak strains differ between these two bone sites. Remembering from the introduction, that literature reports femoral neck adaptation while no adaptation is seen at tibial mid-shaft following an intervention, it certainly appears that modelling based approaches of estimating skeletal loading may prove useful in the future.

Numerical example
To illustrate the differences between presented modelling techniques a simple numerical example will be presented. In the example axial stresses and strains at the mid-shaft cross-section of the tibia of a subject standing on the balls of his feet will be estimated. To make the example more realistic, the dimensions of the subject were taken from the LifeMod 4 database and correspond to a male whose height is 180cm, weight 80kg and age 24 years. In the example the subject extends his ankle joints by 20 • , and tilts his body forward by 10 • .I n this condition the center of mass will be slightly ahead of the foot contact, which will add extra torque at joints to make the example more interesting. Only a right half of a lower body will be considered, meaning: pelvis, right upper leg, right lower leg and a right foot. To make simulation models comparable with hand calculations, which are based solely on ground reaction forces, location of the cross-section at which stresses and strains are computed is preserved together with it's geometry described by a polygon. All models share also the material model for the bone which for simplicity is isotropic and homogeneous, with Young's modulus of 17GPa (Dong & Guo, 2004) and Poisson's ratio equal to 0.3 (Reilly & Burstein, 1975).

Approach based on ground reaction force (model 1)
Let's start the analysis with the simplest 2D (saggital plane) model, which would be the most intuitive first choice (Fig. 7). In the model the mass is concentrated in a point at the hip joint, represented by body gravity force, F g . This location is close to the real center of mass for the human. At the same location stabilization force F s is applied to constrain the model from falling. At the foot-ground contact point, ground reaction force is applied, for simplicity divided into normal, GRF n , and tangential components, GRF f which correspond to static friction. Knowing the mass of the subject, m = 80kg, we can assume for symmetric case, that the load carried by each leg is equal, thus: Fig. 7. Model diagram with tibia cross-section at which the stresses and strains will be computed with the use of ground reaction forces only where g is the gravity constant equal to 9.81066 N kg . The unknown stabilization, F s , and friction forces, GRF f , as well as vertical component of the ground reaction force, GRF n can be solved from the equilibrium equations: where, A x , A y , D x , D y are x and y coordinates of hip joint and foot-ground contact point correspondingly. Substituting to equation 6 numerical values F g can be obtained: Knowing the value of F g from equation 10, value of GRF n can be obtained from equation 8, yielding GRF n = 392.43N. Substituting numerical values to equation 9, F s can be solved as: From equation 7 value of GRF f can be expressed as: GRF f = 19.28N. Knowing all the external forces applied to the model in the equilibrium, internal forces and moments can be computed. The cross-section at which the stresses and strains will be computed is located at point P, as marked in figure 8. Solving of internal normal force, F n , transverse force, F t , and bending moment, M will be based on free body diagram depicted in figure 8. Equilibrium equations for the bottom part of tibia and a foot can be expressed in the cross-section coordinate systemx −ȳ for simplicity: Note that coordinates of points P and D are expressed in the global coordinate system due to the fact, that forces GRF f and GRF n act along global coordinate system axes. The resulting moments from those forces around point P are then computed as the force multiplied the moment arm, which is the distance between contact point and point P and is measured along axis perpendicular to the force direction. In this example we assume that moment has a positive sign if it results in clockwise rotation. In order to solve the equilibrium equations location of point P is required. Point P lies at the neutral axis of the bone cross-section. Thus at this point we will compute the cross-sectional parameters: second moment of area along bending axis, centroid location, and area. Bone cross-section is taken from the generic bone model generated by LifeMod. It is described as a polygon (Fig. 7) with vertices coordinates given in table 1. Area of a polygon can be obtained using the following equation: where n is the number of points, point n + 1 refers to point 1 as the polygon must be closed. Note that the equation 15 holds when the points are ordered in clockwise direction, for the anticlockwise order of points the minus sign has to be added. After computations the area of the bone cross-section is A = 393.07mm 2 . Centroid y coordinate,P x , can be obtained from: Substituting numerical values to equation 16 yieldsP x = −1.84mm. Global coordinates of point P can be obtained from the geometrical dependencies depicted in figure 8B. Note that P x is a coordinate of the point P with respect to point P' not a distance (absolute value of a coordinate), thus sign convention in the equations 17 and 18.
P y = C y + |CP |·cos(α) −P x · sin(α) Substituting numerical values to equations 17, 18 yields: Second moment of area of the bone cross-section aboutz axis can be defined as: where: Substituting numerical values to equations 21 and 22 yields: Iz = 15008mm 4 . Finally, having all cross-sectional parameters, and most importantly coordinates of point P in global coordinate system, equilibrium equations 12, 13 and 14, can be solved for: bending moment, M, normal force, F n and tangential force, F t as follows: Each point within the cross-section of the bone is subjected to stresses caused by the normal force, F n , and bending moment, M. In mechanics tensile stresses are denoted by positive sign, consequently compressive strains are marked by the negative sign. Compressive axial stresses caused by the normal force are uniform within the cross-section and can be calculated using the following relationship: Stresses caused by the bending moment vary within the cross-section with respect to the coordinatex. The relationship between point location and bending stresses is: Total axial strains are sum of stresses caused by the normal force and stresses caused by bending: Substituting numerical values to equation 26 results in: Strains are associated with stresses with a linear relationship called the Hook's law: Table 2 contains axial stresses computed at the vertices of bone's cross-section using equations 26, 27, and 28. In the last column of the table total strains are presented computed using the Hooke's law (equation 30). It can be seen that the posterior part of the bone is subjected to tensile strains while the anterior part is in compression. Bending contribution in the loading is dominating, which can be seen from analysis of table 2.

Inverse dynamics simulation (model 2)
Inverse dynamics simulation is similar to the calculations presented in the previous section. The main difference between the two approaches is the mass distribution within the model and complexity of the mechanical system. Moreover the inverse dynamics model, presented in this section, is three dimensional. This allows to account for bending also in the frontal (coronal) (and transverse) plane(s). In model 2 foot is attached to the ground by a revolute joint. Foot, lower leg and pelvis are interconnected via spherical joints, and pelvis at the center of mass is attached to the ground by planar joint in the sagittal plane. Stabilization force is applied at the center of mass of the body by a constraint restricting translation along the  Table 2. Axial stresses resulting from normal force, σ Fn , bending moment, σ M , total stresses, σ, and total strains, ǫ, at the vertices of the polygon defining the cross-section of the bone global x axis. In order to maintain the desired body position, ankle, knee and hip joints cannot remain unconstrained, thus at each of this joints torques are applied. In the presented case this corresponds to fixing the body segments together rigidly as only static case is discussed. However in the general case, where motion would be studied, those torques would vary over time in order to reproduce subject's motion. The model used in inverse dynamics simulation is presented in figure 9. Models used in simulation can be more complex comparing to the simple mechanical calculations. For this reason natural mass distribution along with inertia properties for each body segment can be achieved. Reasonable guess of geometry and mass distribution can be obtained from LifeMod package. According to the subject specification the foot weight 2.14kg, lower leg mass is 4.19kg, upper leg's mass is 7.47kg and the upper body plus pelvis weight 52.4kg. As only half of the model is considered, the upper body and pelvis mass have to be reduced to 26.2kg to maintain symmetry conditions. Stresses and strains calculations could be performed in analogical way like in the model 1 using the multibody model just to obtain the cross-sectional bending moments and normal force. However use of flexible multibody dynamics allows to compute the stresses and strains directly at any nodal location within the flexible bone. To benefit from the flexible multibody approach, the rigid bone at which stresses and strains are to be computed needs to be replaced with it's flexible version. To obtain the most accurate results bone geometry should be reconstructed from medical images to correspond to the subject, however for this example just the simple shell geometry created by LifeMod will be used. After exporting the geometry from Adams 5 , it can be imported to any finite element software which is capable of generating modal neutral file. For this example ANSYS 12 will be used.
In the finite element software the initial shell geometry is converted to solid body which can be discretized with solid elements. As in the model 1, it is assumed that bone is homogeneous, thus the whole volume of the bone is meshed with 5mm in size four-node tetrahedral linear solid elements. Each of the elements is assigned Young modulus of 17GPa and Poisson's ratio of 0.3. In order to perform modal analysis rigid beams are added at the distal ends of tibia connecting surface nodes with knee and ankle joint centres respectively. Craig-Bampton (Craig & Bampton, 1968) modal reduction allows to represent large finite element models with the use of several deformation modes, thus remarkable reducing the computational effort in Fig. 9. Illustration of the three-dimensional model used in the inverse dynamics simulation. The model consists of four segments: foot, lower and upper leg, and a pelvis.
the simulation software. Knee and ankle joint center location nodes are used as boundary nodes in the modal analysis. For the current example 42 deformation modes will be computed.
In the simulation however only 36 modes are usable, due to the fact that first six deformation modes represent rigid body motion, which is accounted for already in multibody formulation.  Table 3. Total axial stresses, σ, and total strains, ǫ, at the vertices of the polygon defining the cross-section of the bone obtained from inverse dynamics model

Forward dynamics simulation (model 3)
Forward dynamics simulation with the use of flexible multibody dynamics is based on the model presented in inverse dynamics example. However in the model 3, torque actuators at the joints are removed, and in turn muscle actuators are added. Muscle actuators are represented by straight line of action between muscle insertion points. The contraction of each muscle is controlled by the proportional-integral-derivative (PID) controller, which aims at replicating muscle contraction patterns obtained in the inverse dynamics. The model is actuated by 16 muscles: tibialis anterior, soleus, gastrocnemius medialis, gastrocnemius lateralis, rectus femoris, vastus lateralis, vastus medialis, biceps femoris (two actuators), illiacus, semitendinosus, gluteus medialis (two actuators), gluteus maximus (two actuators), adductor magnus. Muscle placement is depicted in figure 10. Adductor magnus is not visible in the picture as it is on the inferior site of the femur. Forward dynamics simulation requires significant amount of input parameters, among them muscle contraction patterns, that are used as a reference signal for muscle controllers. The easiest way to obtain muscle contraction splines is the inverse dynamics simulation. During the inverse dynamics model is driven by the motion captured from a subject. In the example case, there is no motion, so during the inverse dynamics the model is simply fixed in desired orientation. After gathering muscle contraction splines, which in this case are just initial muscle length values, the forward dynamics simulation can be performed. External constrains are removed, meaning that the joints are free to move if no muscle force is produced. However support hinge at the foot, planar constraint at the pelvis and translation along x axis constrain at the center of mass are preserved. Simulated ground reaction forces are: GRF n = 392.22N and GRF f = 20.57N for the vertical and frictional components respectively. Stresses and strains at the mid-shaft cross-section of the tibia are presented in table 4.

Summary of the numerical example
It is worth noting that the loading pattern obtained from model 3 is completely different to the ones obtained using models 1 and 2. Model 3 predicts that the posterior part of tibia is subjected to compression while the anterior part is in tension. The magnitude of strains and stresses is around 7.62 times smaller for model 3 compared to model 2. The differences can be explained by analysing differences in mechanics of model 2 and 3. The first one behaves like a rigid structure, where sagittal ankle torque is carried as bending moment by the tibia. Model 1 behaves the same way. Comparing the bending stresses at the considered cross-section to normal stresses one can notice that the bending stresses are up to 19 times higher than the normal stresses. In contrast, in model 3 tibia, foot and muscles create a mechanism. Due to the fact that ankle joint is almost frictionless it carries almost no torque at all, and the contact force is transmitted through Achilles tendon, and through ankle joint just as an axial and tangential force. Tibia is then subjected to only transverse and normal forces at the knee and ankle joints, and the bending moment comes mostly from the activity of soleus (545.72N), vastus lateralis (67.6N), vastus medialis (2.39N) and rectus femoris (0.08N) muscles. All of those muscles while contracting create bending moment which creates compressive strains at the posterior site of tibia. Tensile strains from bending moment actually reduce the compressive loading from axial force at the anterior site of tibia. This example demonstrates 1) that internal forces should not be neglected in estimating skeletal loads and 2) that muscle contraction can provide  Table 4. Total axial stresses, σ, and total strains, ǫ, at the vertices of the polygon defining the cross-section of the bone obtained from forward dynamics model stress shielding for the bone Sverdlova & Witzel (2010). When it comes to the external forces they were predicted correctly by all three models.

Conclusion
Estimating lower limb skeletal loading may be useful e.g. in association to planning osteogenic exercise interventions. Although estimating skeletal loading from ground reaction forces is tempting due to the simplicity of the measurement, more information may be gained from more elaborate approaches. As was demonstrated by the numerical example, accounting for mass distribution in the model (as shown in model 2) can lead to 143% difference in strains. Further difference of up to 19 times can be noted when internal forces are considered in addition to external forces (i.e. joint torque actuators are replaced with muscle models). Even though all of the methods predicted almost identical stabilization force and produced the same ground reaction forces, internal loading was very different. For example in implant design only relatively realistic strain estimates may be expected to be useful. Implant loosening is a well known phenomenon (van Loon et al., 1999) and designing the implant in such a way that appropriate loading levels are attained can be therefore argued to be quite relevant. Furthermore, if a specific (for example clinically relevant) bone site is to be targeted with an exercise intervention, it may be necessary to the estimate the loading in a site specific manner. During last couple of years there has been an increasing recognition that problems arising in biology or related to medicine really need a multidisciplinary approach. For this reason some special branches of both applied theoretical physics and mathematics have recently emerged such as biomechanics, mechanobiology, mathematical biology, biothermodynamics. This first section of the book, General notes on biomechanics and mechanobiology, comprises from theoretical contributions to Biomechanics often providing hypothesis or rationale for a given phenomenon that experiment or clinical study cannot provide. It deals with mechanical properties of living cells and tissues, mechanobiology of fracture healing or evolution of locomotor trends in extinct terrestrial giants. The second section, Biomechanical modelling, is devoted to the rapidly growing field of biomechanical models and modelling approaches to improve our understanding about processes in human body. The last section called Locomotion and joint biomechanics is a collection of works on description and analysis of human locomotion, joint stability and acting forces.