Development and implementation of a coupled computational muscle force optimization bone shape adaptation modeling method.

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN BIOMEDICAL ENGINEERING Int. J. Numer. Meth. Biomed. Engng. 2015; e02699 Published online 18 February 2015 in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/cnm.2699

Development and implementation of a coupled computational muscle force optimization bone shape adaptation modeling method C. S. Florio*,† Department of Mechanical and Industrial Engineering, Newark College of Engineering, New Jersey Institute of Technology, University Heights, Newark, NJ 07102, U.S.A.

SUMMARY Improved methods to analyze and compare the muscle-based influences that drive bone strength adaptation can aid in the understanding of the wide array of experimental observations about the effectiveness of various mechanical countermeasures to losses in bone strength that result from age, disuse, and reduced gravity environments. The coupling of gradient-based and gradientless numerical optimization routines with finite element methods in this work results in a modeling technique that determines the individual magnitudes of the muscle forces acting in a multisegment musculoskeletal system and predicts the improvement in the stress state uniformity and, therefore, strength, of a targeted bone through simulated local cortical material accretion and resorption. With a performance-based stopping criteria, no experimentally based or systembased parameters, and designed to include the direct and indirect effects of muscles attached to the targeted bone as well as to its neighbors, shape and strength alterations resulting from a wide range of boundary conditions can be consistently quantified. As demonstrated in a representative parametric study, the developed technique effectively provides a clearer foundation for the study of the relationships between muscle forces and the induced changes in bone strength. Its use can lead to the better control of such adaptive phenomena. Copyright © 2015 John Wiley & Sons, Ltd. Received 20 June 2014; Revised 14 November 2014; Accepted 17 November 2014 KEY WORDS: cortical bone adaptation; muscle force magnitude; gradientless shape optimization; gradient

projection optimization; parametric studies

1. INTRODUCTION The investigation of the strength adaptations of bone and the application of the resulting insight learned about these phenomena to the prevention of the detrimental effects of decreased bone strength has been the focus of much research. Numerous clinical investigations, experimental studies, and computational models have been developed to examine the interrelationships between alterations to the mechanical environment of a bone and the ensuing changes to its shape and density and, therefore, strength. Such research has resulted in a wide array of proposed techniques to counteract the losses of bone strength that arise from age or disuse with an equally wide range of reported effectiveness. Despite the significant knowledge gained about the effects of mechanical factors on bone strength from these studies, the direct comparison of various mechanical-based bone loss mitigation methods has been limited. The modeling method developed in this work enhances and expands the set of currently available tools with which mechanical means of targeted bone strengthening can be better analyzed, understood, and designed. *Correspondence to: C. S. Florio, Department of Mechanical and Industrial Engineering, Newark College of Engineering, New Jersey Institute of Technology, University Heights, Newark, NJ 07102, U.S.A. † E-mail: [email protected] Copyright © 2015 John Wiley & Sons, Ltd.

(1 of 29) e02699

e02699 (2 of 29)

C. S. FLORIO

Experimental investigations of mechanically induced bone strength changes using animal subjects have employed both the direct control of loading conditions and the direct measurement of the induced strains on the bone surfaces. The mechanical environment of the bone of interest in these studies is often altered through the removal of all loads [1], the removal of specific muscle forces [2], or the application of external forces through surgically implanted mechanisms [3–6]. The effects of more realistic loading patterns, modes, and activities have also been studied by subjecting animals to controlled ambulation through the use of treadmills, running wheels, or swimming tanks and comparing the bones of the exercised subjects to those of immobilized or non-exercised ones [7–10]. In many of these animal experiments, strain gages have been surgically attached to the bones for the direct measurement of the local stress state during load application [5, 11–14]. The sacrifice of the subjects after a particular amount of time and the study of cadaveric animal bones have allowed for the direct measurement of changes in bone structure and shape [1, 3, 9, 15, 16]. Despite these invasive controlled-load applications and measurements, the identification of absolute trends and correlations from these animal experiments has been difficult because of the various ages, species, sizes, bones, and study durations used [1, 4, 9]. Because the loss of strength in the long bones of adult humans is typically the clinical issue intended to be addressed, the examination of the targeted population can reduce much of the variation that has limited the application of observations from animal experiments. Disuse and specified loading patterns have been studied both on Earth and in Space through the examination of those subjected to bed rest, casting, or microgravity, enrollees in prescribed exercise programs over a particular study period, and long-term sports participants [17–25]. While these studies have been able to identify resistance exercises as beneficial to reducing bone loss or even improving bone strength [18, 26, 27], the specific activities that are most effective at strengthening fracture-prone regions of bone have yet to be identified through these clinical investigations. While strain gages have been surgically attached to human bones [28], the loading intensity in these clinical investigations is most often calculated using non-invasive global measurements such as whole body ground reaction forces [25] and accelerations [22]. Mixed loading mode regimes are often employed [21–23], preventing the isolation of the effects of singular loading modes on bone strength adaptations. Additionally, difficulties with subject compliance to the prescribed exercise programs [21, 22, 29] and the inability to control external influences and variations in normal activity levels between subjects [29] have prevented these clinical studies from providing quantitative comparisons of the strengthening effectiveness of the exercise routines examined. Additionally, methods to non-invasively and simultaneously measure local and global changes in bone strength, composition, and shape are scarce. Computational models have often been used in the analysis of inert mechanical systems as a means for controlled comparisons of the effects of system parameter values because they eliminate the variations inherent in physical experiments and allow phenomena to be revealed that are often difficult to measure experimentally. Therefore, in the study of the effects of mechanical factors on bone strength adaptation, numerical modeling methods can provide the repeatability lacking in animal experiments and human clinical studies and can allow for more detailed measurements of changes on bone strength, shape, or composition throughout the entire system studied, which may be difficult to achieve through imaging, histological or morphological evaluations of test subjects. Additionally, computational investigations allow for more efficient and thorough parametric evaluations, as system properties and boundary conditions can easily be altered, and study times are not dependent upon rates of biological phenomena. Despite these benefits, the use of computational modeling methods in the comparison and evaluation of mechanically based bone strengthening methods has been limited. Because of the complexity of the geometry and boundary conditions in the musculoskeletal systems of interest, significant simplifications are often applied. The bone geometry is reduced in detail [30, 31] and extent of region studied [32, 33], often to a single bone within the system [34, 35] or even a circular or elliptical two-dimensional transverse cross-section [36–39]. Boundary conditions are simplified by combining muscle forces [33, 34, 40] or reducing the muscle force contact region to a single point [19, 35]. However, because bone strength adaptation is driven by the stress distribution over the bone studied, the physical representation of the examined system directly impacts the prediction of the adaptive response. Copyright © 2015 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Biomed. Engng. 2015; e02699 DOI: 10.1002/cnm

COUPLED MUSCULOSKELETAL OPTIMIZATION

e02699 (3 of 29)

More complete representations of the boundary conditions of the bone have been shown to be more important than more complete representations of the bone geometry in improving the accuracy of these computational simulations [33, 34, 40]. Often, however, computational models of bone strength adaptation include either detailed representations of the geometry or boundary conditions, but not both. Recommended modeling features from various studies can be combined for improved model accuracy. Individual, rather than grouped, muscle forces with non-singular point regions of application are preferable, as more gradual spatial stress gradients, generally lower peak stresses, and, therefore, more conservative adaptation predictions result [33, 41]. Because muscles often act on multiple bones and neighboring bones transfer loads and stresses throughout the system, multibone models remove the need for the artificial constraints at cut bone sections or joint ends that are required in more simplified models [34, 42, 43]. The representation of axial variations along the bone length is important either through three-dimensional models [44] or twodimensional models using longitudinal, rather than transverse, cross-sections. Finally, both cortical and cancellous bone tissues should be included for more realistic simulations of the mechanical response of the bone structure to applied forces [35, 45]. Many computational bone strength adaptation models are often heavily reliant upon experimentally based parameters. In the development of these models, parameters are commonly altered until the model predictions match experimental observations [35, 39, 46, 47]. However, a model developed in this fashion is typically only applicable to the system and conditions used in deriving its parameters, often preventing its use in comparing varying system conditions, loads, and configurations, such as in the wide range of currently available mechanically based bone strengthening techniques. A computational modeling method is developed in this work that addresses the current limitations in the study of bone strength adaptation through a more complete representation of the system studied and a more consistent means of predicting the strength adaptations, independent of experimentally-based or condition-based parameters. Through the coupled use of numerical optimization and structural analysis methods, a more universally applicable tool is created with which to quantitatively analyze, compare, and develop techniques to counteract the losses of bone strength that occur with disuse or age. Bone shape adaptation is included in this work, rather than bone density adaptation, because such shape changes are thought to have a more significant and stable impact on bone strength [48, 49]. While the modeling methods developed may be applied to any system, to understand their capabilities, the techniques are implemented to examine the strength changes in the tibia bone within a human leg system under isometric loading conditions. The leg is the focus of study in this work because its muscles can generate a wide range of force magnitudes, thereby evaluating the potential span of modeling ability. In addition to such parametric studies, the methods created in this work can be used to help improve the understanding of the relationships between boundary conditions and bone strength adaptations. 2. MODEL DEVELOPMENT The computational modeling method developed in this work couples two models previously developed by the author which predict the individual muscle force magnitudes in a multisegment system that create a given resultant isometric force [45] and simulate mechanically based bone shape adaptation [50]. The model functions completely within a commercial finite element program with significant customization through user-defined subroutines. The basic model function is shown in the flowchart in Figure 1. The system definition includes identification of the included bone segments and their geometries, material properties, and their contact behavior at the joints, pertinent individual muscles and their physiological cross-sectional areas and locations and areas of attachment on the included bone, as well as all included joint angles. Boundary conditions include displacement constraints based on the conditions simulated and the desired net isometric resultant force. The magnitudes of the included individual muscle forces required to generate this resultant force are then calculated. These muscle forces and displacement constraints are directly applied to a discretized model Copyright © 2015 John Wiley & Sons, Ltd.


e02699 (4 of 29)

C. S. FLORIO

Figure 1. Basic function of developed modeling method.

of the bone system, and a static structural finite element analysis is performed. Based on the resulting stress distributions over the periosteal (outer) and endosteal (inner) surfaces of the targeted bone, the amounts of local material accretion or removal on these surfaces are calculated, and their profiles are changed by moving individual surface nodes. Because this alters the shape of the elements containing the surface nodes, which can adversely affect the accuracy of the finite element structural analysis, node smoothing methods are implemented both for the surface nodes that are moved by the shape strength adaptation model and for the interior nodes that connect the changing outer and inner surfaces. The process is repeated until predefined convergence or stopping criteria is reached. The model can be used to study the effects of the repeated application of a single net resultant force or a sequential series of such forces in any multisegment musculoskeletal system. Copyright © 2015 John Wiley & Sons, Ltd.



e02699 (5 of 29)

2.1. Two main modeling components The function of the modeling techniques used to determine the individual muscle force magnitudes and the mechanically based bone strength/shape adaptation will be summarized. Details can be found elsewhere [45, 50, 51]. Each of these developed methods can be implemented independently of each other based on the desired modeling goals. However, their coupled usage, as in the current work, allows for a more automated and complete study of the mechanical factors that influence bone strength adaptation. 2.1.1. Determination of individual muscle force magnitudes. The resolution of the individual muscle force magnitudes in a typical musculoskeletal system cannot be solved with standard conservation of momentum equations because the number of muscles usually exceeds the number of degrees of freedom in the studied system. Therefore, mathematical optimization methods are often employed. In the current work, a gradient-based computational optimization method was implemented based on a previously developed method [52]. A code was written to calculate the set of muscle force magnitudes in a given system that minimize the sum of the squares of the muscle stresses induced in the generation of a net isometric resultant force [45]. The method was selected based on a review of previously published studies of similar conditions. Comparisons of muscle force magnitudes predicted using different optimization methods and goals have found that the particular optimization method used in this work is relatively insensitive to the inclusion of higher order terms [52, 53] and that the minimization of squares, rather than higher orders, of muscle stresses produces superior correlations to experimental electromyograph measurements when estimating muscle forces created in isometric exercises [54]. Muscle stresses were defined as the ratio of the magnitude of the force generated by the muscle along its axial length to its physiological cross-sectional area, the average cross-sectional area of the muscle normal to its axis. The optimization problem is written as follows: 2 n Fl (1) Minimize Z ðFÞ ¼ ∑ l¼1 PCSAl Constrained by gðFÞ≤ 0 where gðFÞ is in the form gðFÞ ¼ AT F b

(2a)

(2b)

The Z(F) is the optimization goal, and F is the array of unknown muscle force magnitudes. g(F) is an array of length n + m that contains the system constraints. The first n (inequality) constraints restrict the muscle force magnitudes to nonnegative values and are of the form: g q ðFÞ≤ 0; q ¼ 1…n

(3a)

where g q ðFÞ ¼ Fl

(3b)

The last m (equality) constraints ensure that the momentum balance equations are satisfied and are of the following form: g q ðFÞ ¼ 0;

q ¼ n þ 1…n þ m;

(4a)

where gq(F) represents the following moment balance: n

Mm ¼ ∑ rlm Fl

(4b)

l¼1

Copyright © 2015 John Wiley & Sons, Ltd.


e02699 (6 of 29)

C. S. FLORIO

When examining static or quasi-static activities as in the current study, this optimization formulation has been shown to eliminate the need for additional constraints to limit the maximum values of predicted muscle forces [52], thus removing the undesirable artificial effects on predicted muscle force distributions that has sometimes been known to occur when using such restricting constraints [55, 56]. A gradient projection optimization method is employed in this work because it has the properties of maintaining feasibility in each iterative solution with a feasible initial guess and search direction and of finding the global optimum when the optimization goal Z(F) is convex, as in this case (Equation (1)). It is a form of the steepest descent method, which defines the search direction as normal to the optimization goal Z(F) and a function of its gradient ∇Z(F). By projecting this gradient onto the surface of the space of possible solutions, only the extreme locations of the feasible solution space are searched in this method, providing the sufficient conditions for any point along this search direction to be an optimum. This also reduces the inequality problem to an equality one, so that only the constraints that are exactly satisfied (equal to zero within tolerance) at the possible optimal point are included in each iterative solution. Because the constraints are linear, Rosen’s method [57–59] is used to determine the step size along the search direction, which was defined through adjusting the set of inequality constraints that are exactly satisfied to find the next nearest potential optimal solution point. The solution proceeds iteratively until a global optimum is found when the projection of the gradient of the optimization function at the solution point is equal to zero (within tolerance). This is a representation of the Kuhn–Tucker conditions, and it provides the necessary conditions that the solution point found is an optimal solution. The code was written [45, 51] to numerically implement the gradient projection method [59–61] using MATLAB (Mathworks, Natick, MA, USA) [62] software to take advantage of its extensive built-in mathematical functions, especially its linear algebra libraries. To improve numerical stability, tolerance bands were used in the selection of “exactly satisfied” constraints and QR factorization was employed to ease susceptibility to ill-conditioning and roundoff error. Input parameters include the number of individual muscles in the system studied, their physiological cross-sectional areas, and their moment arms to the center of each system joint they cross (such as hip, knee, and ankle) and to the location of the net resultant force generated. An initial feasible solution is also required. The angular momentum balance equations are integral to the code and are the only part which must be manually altered based on the particular system studied. The individual muscle force magnitudes and the net joint moments are returned upon completion of the solution. The developed code could be used independently to study the muscle activity under various conditions, or it can be combined with subsequent analyses of the effects of these individual muscle forces. 2.1.2. Simulation of bone shape strength adaptation. Similar to many structural shape optimization models and previously developed biological adaptation models, a gradientless optimization method is employed to simulate the changes to the bone shape in response to the stress distributions over the bone surfaces that result from the applied muscle forces within the multisegment system studied. In such methods, the shape iteratively changes at discrete points over the bone’s surface as a measure of the local stress state at each of these points progresses towards a pre-defined goal, which is often a reduction in the variation of the local stress state over the bone’s surface. While gradient-based methods similar to the one used to determine the individual muscle force magnitudes have been employed to simulate these phenomena, their implementation is difficult in finite element and boundary element-based numerical shape optimization models and prohibitive for large design regions, such as the three-dimensional geometries in this work [63–65]. Therefore, while gradientless optimization methods sacrifice mathematical rigor, solution uniqueness, and repeatability for ease of implementation [58, 59, 66], because many bone shape adaptation and structural shape optimization models have initial geometries similar to their optimized designs, an object with enhanced strength often results [67]. The basic form of these computational bone adaptation models iteratively changes the location x of each discrete point l on the optimizing surface by an amount U in the following manner:

xkþ1 ðlÞ ¼ xk ðlÞ þ U k ðlÞ U k ðlÞ ¼ α σ k ðlÞ σ ref Copyright © 2015 John Wiley & Sons, Ltd.

(5a) (5b)



e02699 (7 of 29)

In this work, the nodal location xk(l) is reduced to a single component of the position vector, the local radial coordinate value, Uk(l) is referred to as the growth driver expression, and the index k is used to track the iteration count. While many different measures of the local stress state σ k(l) have been used [68–70] and much research has been undertaken to determine which measure is most likely to be used by the bone cells themselves to invoke the adaptive phenomena [13, 71, 72], trends in shape changes with varying mechanical loading modes are consistently predicted regardless of the particular measure employed [73]. In the current adaptation simulation, the strain energy density is used because, as a scalar, it eliminates the need to account for the relative orientations of the bone geometry, growth direction, and stress/strain tensor components, and, as an energy measure, it represents the energy imparted to the system that drives the mechanically based cellular bone strength adaptation phenomena [74, 75]. Because of the mathematically lax nature of this gradientless type of optimization method, the model parameters, σ ref, α, and the stopping criteria are chosen by the model developer. Often, these parameters are based on experimental observations [35, 39, 46, 76], the specific conditions studied [64, 77, 78], or to control model stability and mesh distortion [79–81]. However, this makes the model suited only for the conditions under which the parameters were selected, requiring different parameter values not only for the study of different conditions but sometimes also for different portions of the design region, as multiple orders of magnitude variations in the stress measure are typical for large design regions, such as the surfaces of long bones [50, 51]. Because these gradientless optimization methods do not converge to a mathematically optimal solution but merely one closer to the optimization goal than the initial configuration, the use of system-specific or experimentally based parameters limits their universal applicability, preventing quantitative comparisons from being derived from model results. With the modeling intent in the current work of a systematic comparison of the effects of a wide variety of mechanically based bone strengthening methods, more global applicability was required. A condition-independent, system-independent method to consistently define the reference value, σ ref, and scaling factor, α, in Equation (5) was developed [50, 51]. Rather than the “universal” value often desired by many bone adaptation models [82], a surface-averaged value of the current stress state (Equation (7)), similar to many structural shape optimization models [78, 83, 84], is used, ensuring that the reference value is always suited to the current conditions modeled. While normalizing the growth driver expression in Equation (5b) by the reference value σ ref [79, 85] aids in limiting the variation in calculated nodal growth over the design region (bone surface), node-specific scaling factors, α(l), must often be implemented to prevent mesh distortion. To eliminate this arbitrary scaling of nodal growth over the bone surface in the current work, the growth driver expression Uk(l) is normalized in the developed model by the standard deviation of the local nodal stress state over the design region (bone surface) (Equation (8)). This bone adaptation model, based on statistical standardization of an arbitrary set of data, reduces the nodal growth over a large design surface to a more consistent range that no longer requires individual nodal modifications based on current, local conditions. Thus, the bone shape adaptation model used in this work is of the following form: xkþ1 ðlÞ ¼ xk ðlÞ þ αU k ðlÞ U k ðlÞ ¼ α

SEDk ðlÞ SEDEleAWSurfAvg k SEDkEleAWStDev

(6a) (6b)

where the element area weighted surface average and standard deviation of the local strain energy density are defined as in Equations (7) and (8). M

∑ SEDk ðeÞAk ðeÞ ¼ e¼1 SEDEleAWSurfAvg k

M

∑ Ak ðeÞ

(7)

e¼1



e02699 (8 of 29)

C. S. FLORIO

SEDEleAWStDev k

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uX 2 u M EleAWSurfAvg u A ð e Þ SED ð e Þ SED k k k u u ¼ u e¼1 M u X u ð M1 Þ Ak ðeÞ t

(8)

i¼1 M

With N X

SEDk ðeÞ ¼

SEDk ðnÞ

n¼1

N

(9)

The surface average used in the reference value and the standard deviation used to normalize the growth driver expression in Equation (6b) are weighted by the area A(e) of the exposed face of each of the elements e on the optimizing surface in the design region. Each of the M “elemental” strain energy densities, SED(e), in the design region are calculated as the average of the strain energy densities, SED(n), at each of the N nodes that form this element face. A mapped mesh using eight-node hexahedral elements is applied to the initially cylindrical design region of a typical long bone in this work. The use of a cylindrical coordinate system simplifies many calculations in the developed model. For example, A(e) is calculated as follows: AðeÞ ¼ r avg ΔθΔz

r avg ¼

4 1X ri 4 i¼1

(10a)

(10b)

Δθ ¼ θj θjþ1

(10c)

Δz ¼ ðθk θkþ1 Þ

(10d)

where i indicates the nodes on the four corners of the face of the element on the optimizing, bounding surface., j refers to neighboring nodes along the theta-direction, and k denotes neighboring nodes along the z-direction. Because the mapped mesh is regular in theta-coordinate and z-coordinate, the element surface areas remain rectangular throughout the optimization, and this method of area calculation remains valid. Because the variation in nodal growth over the bone surface is controlled through the definition of the reference value and normalization, the scaling factor, α, in Equation (6) is no longer dependent upon the local nodal conditions as it is in Equation (5). Instead, to maintain mesh integrity during the optimization, this value was defined as 5% of the initial radial element size in the design region, similar to suggestions in the literature [79, 86]. The code was written in FORTRAN to perform this shape optimization/bone strength adaptation simulation from directly within a customized version of a commercial finite element program, ANSYS Version 12.1 (ANSYS Inc., Canonsburg, PA, USA) [87]. The local nodal and elemental strain energy densities and the ensuing amounts of nodal growth Uk(l) are calculated in the userdeveloped code from nodal coordinate and stress/strain tensor data computed by the built-in finite Copyright © 2015 John Wiley & Sons, Ltd.



e02699 (9 of 29)

element structural analysis routines. Nodal coordinates are then adjusted such that positive Uk(l) values result in material accretion and negative Uk(l) values result in material resorption. Output data including nodal coordinates, the amounts of nodal growth per iteration, and local nodal and global surface-averaged strain energy densities, are recorded for tracking of the optimization progress and for subsequent calculations. Because the nodal coordinates are recorded by the commercial code with respect to the global system coordinate system and the nodal growth is calculated in the user subroutine normal to the bone surface in a local coordinate system, coordinate transformations [88] are performed as needed. The adaptation process is simulated in the developed model for both the endosteal (inner) and periosteal (outer) bone surfaces and proceeds iteratively as the boundary conditions are repeatedly applied to the changing bone geometry and new stress state distributions are revealed until a predefined stopping criteria is achieved. Because the positions of the surface nodes in the design region are moved in this model to represent the shape/strength adaptations, extreme mesh distortion can occur as the simulation progresses. This can affect the accuracy of the structural finite element analysis that is used to drive the adaptation predictions. Therefore, like many similar models [39, 79, 81, 89–92], mesh adjustment methods are employed on the exterior surfaces and in the interior volume of the bone geometry design region. 2.2. Nodal smoothing methods Mesh adjustment methods used in previous bone adaptation and structural shape optimization models usually involve remeshing [39, 93] and/or nodal smoothing [30, 94]. Because the developed routine is based on the changes in the positions of the surface nodes, it is desirable that the number of surface nodes remained constant in order to easily quantify the changes to the surface profile. Therefore, nodal smoothing methods, which adjust the locations of existing nodes rather than adding and removing nodes, employed in the current work. Two distinct methods are used, one for adjusting the interior nodes to maintain a more uniform distribution throughout the bone volume as the surfaces change shape and one for adjusting the nodes in the design region to accommodate for any large differences among the calculated growths of neighboring nodes. Each of the mesh adjustment methods discussed is computed and executed through author-developed user-defined subroutines within the customized version of the commercial finite element code used. 2.2.1. Interior node smoothing. Although the outer surface nodes are the focus of the modeling technique development presented in this work, the importance of modifying the locations of the interior nodes is directly related to the accuracy of the structural finite element analysis from which the predictions of the surface profile changes are made [81, 89, 95]. Without such smoothing, extreme distortion of the elements near the adapting surfaces can occur (Figure 2). A spring-based smoothing method [96] is employed in the developed model to maintain the relative spacing of interior nodes despite changes to the surface profiles during the iterative optimization/adaptation process. In such a method, the nodes of the meshed volume are imagined

Figure 2. Mesh distortions due to growth of outer surface without interior node smoothing. Copyright © 2015 John Wiley & Sons, Ltd.


e02699 (10 of 29)

C. S. FLORIO

to be connected via a system of springs that run along the element edges. The initial mesh configuration, before any shape changes, is assumed to be in an equilibrium state. This state is perturbed as the surface nodes are moved, simulating material removal or accretion, and the individual springs connected to these surface nodes are “compressed” or “stretched”, respectively. Dimensionless “forces” are generated on the surface nodes by these perturbed springs and then propagated through this imaginary spring system as follows: Fi ¼

ni X

k ij Δxj Δxi

(11)

j

The index i represents the node with the applied force Fi, generated by the movement of a surface node i, while j represents the index of a neighboring node i. Each node has ni connected neighboring nodes. The effective spring constant, kij, of the spring connecting node i to node j is inversely proportional to the pre-smoothed, pre-perturbed distances between neighboring nodes: 1 k ij ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi xi xj • xi xj

(12)

With the amount of displacement of the exterior nodes known from the bone shape adaptation model, the required changes to the positions of the interior nodes resulting from this smoothing of the meshed geometry are calculated through a system of equations patterned after force equilibrium equations for a system of connected linear springs. This ensures that the initial spacing between interior nodes is maintained throughout the shape optimization (adaptation) process. In the current work, the mesh in the adapting region is composed of eight-node hexahedral elements. Shape changes are limited to the radial direction on the inner and outer surfaces of the hollowed region of the cylindrical long bone geometry. Thus, the adjustments of the interior nodal positions are limited to radial lines connecting the inner and outer surface nodes that have common local axial and angular cylindrical coordinates. The imaginary spring system in the interior smoothing method employed is, therefore, a one-dimensional set where each interior node is connected to two springs, and each surface node is connected to one spring. Therefore, the coefficient matrix in the system of force equilibrium equations (Equation (11)) is reduced to a simple band matrix. An example of the implementation of the node smoothing routine is shown in Figure 3. Either net growth or decay can occur in a three-dimensional initially uniform circular cylinder based on the value of an applied uniform surface pressure when the top and bottom surfaces are fixed. A cross-section through the mid-plane of the cylinder demonstrates the maintenance of the uniform nodal spacing.

Figure 3. Interior node smoothing with element size (a) reduction and (c) expansion. The initial mesh is shown in (b). Copyright © 2015 John Wiley & Sons, Ltd.



e02699 (11 of 29)

2.2.2. Exterior node smoothing. The developed bone shape adaptation model in Equation (6) allows each node to move independently of the conditions at any other node, even its direct neighbors. While the majority of the calculated nodal growth values per optimization iteration are within the same order of magnitude in local regions, stress concentrations or changes in boundary conditions can lead to significant mesh distortions. A Laplacian smoothing method is implemented in this work to reduce the extreme growth variations that may arise between neighboring nodes. This kind of method moves nodes based on a weighted effect of nearby conditions [97, 98] following the basic form: N X

Xk

Smooth

ðl Þ ¼

Xk ðiÞwk ðiÞ

i¼1 N X

(13) wk ðiÞ

i¼1

where XSmooth ðlÞ is the smoothed position of the node l. The position of node l is adjusted by a weighted k effect wk(i) of the positions of the Xk(i) of the neighboring N elements or nodes. The weighting factor can be any parameter that varies smoothly, such as positions, displacements, strains, or stresses, or it can be an arbitrarily selected value 0 < wk(i) ≤ 1. The method is repeated iteratively until a sufficient amount of smoothing is achieved. In this work, the weighting factor in the Laplacian smoothing of the outer surface nodes is selected through a previously developed method based on the gradients of the nodal strain energy densities [94]. While the previously developed method adjusts the nodes in each of the three Cartesian coordinates, because the shape adaptation model in this work only alters the local radial coordinates of the surface nodes, the smoothing method is applied to the calculated change in radial position Uk(l) in Equation (6b) rather than the nodal coordinates themselves Xk(l) in Equation (6a). However, to capture the complete variation of the strain energy density that could lead to the discontinuities in the predicted growth at neighboring nodes, the magnitude of the full strain energy density gradient at each node is used in calculating the weighting factor. Thus, the surface smoothing routine implemented is of the following form, with l as the index of the node being smoothed and i, the index of the N neighboring nodes plus the node itself: Nþ1 X

ðlÞ ¼ U Smooth k

wk ðiÞU orig k ðiÞ

i¼1 Nþ1 X

(14a) wk ðiÞ

i¼1

Where:

wk ðiÞ ¼

8 1 > > >
> > : ∑ jð∇SEDk ðiÞÞj

if maxi jð∇SEDk ðiÞÞj > 0

(14b)

i¼1

Because of the cylindrical bone geometry in this work, the "magnitude" of the strain energy density gradient is calculated in local cylindrical (θ, r z) coordinates as follows: ∇SEDðiÞ ¼

∂SEDðiÞ ∂SEDðiÞ ∂SEDðiÞ θþ rþ z ∂θ ∂r ∂z

s ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ∂SEDðiÞ 2 ∂SEDðiÞ 2 ∂SEDðiÞ 2 þ þ j∇SEDðiÞj ¼ ∂θ ∂r ∂z Copyright © 2015 John Wiley & Sons, Ltd.

(15a)

(15b)


e02699 (12 of 29)

C. S. FLORIO

For most nodes, a fourth-order accurate finite difference calculation of the directional gradient is used [99]. However, in some cases, such as for the radial gradient as well as the gradient in the z-direction for nodes at the ends of the design region, a second-order accurate finite difference or central difference calculation is employed [99]. As in a typical Laplacian smoothing method, the developed surface smoothing routine is iteratively implemented. The smoothing is stopped when the maximum change in the value of Uk(l) from the previous smoothing iteration over all of the nodes being smoothed is less than 10% of a limit imposed on nodal growth per iteration. In some regions, this Laplacian smoothing method is not sufficient to reduce the large discontinuities in growth calculated for neighboring nodes using Equation (6b), especially near non-moving nodes at the boundaries of the design region. Therefore, the shape adaptation optimization model is not applied in these regions, and the locations of these nodes are simply adjusted to create a smooth transition from the fixed, non-moving nodes to the nodes in optimization design region. The lengths of the transition regions on the outer and inner surfaces are equal to the initial outer or inner bone diameters, respectively, following Saint-Venant’s principle of the effect of boundary conditions on local stress distributions [100]. Because these discontinuities tend to be more significant on the outer surface because of directly applied boundary conditions, a quadratic fit is used on the outer surface while a linear fit is used on the inner surface. These fits are created only along axial lines of the mapped mesh using nodes with common angular coordinates. The implemented smoothing methods successfully removed the mesh distortions that resulted from the independent calculation of the local adaptation and movement of individual surface nodes based on the gradientless shape optimization method developed in this work. The effectiveness of these methods is demonstrated in the test case shown in Figure 4 where four surface forces result in the regions of large growth on an initially uniform circular cylinder with a fixed base. 2.3. Stopping the simulation The described process of calculating the stress distributions in the design region of the bone due to the applied boundary conditions, predicting the local amounts of material accretion or resorption, and adjusting the mesh to reduce distortions imposed during the alteration of the surface profiles is repeated iteratively. The bone shape strength adaptation simulation is completed when either a

Figure 4. Exterior node smoothing (a) before and (b) after implementing the developed smoothing routines. Copyright © 2015 John Wiley & Sons, Ltd.



e02699 (13 of 29)

predefined size limit is reached or a sufficient improvement of the uniformity of the stress state in the design region is achieved. Size limits are imposed to restrict the amount of bone growth or decay that could reasonably be predicted by the developed model. Four unique limits are defined: maximum outer diameter, minimum inner diameter, and maximum and minimum cortical (wall) thickness. Each of these limits, except the minimum inner diameter, is expected to be selected from published extreme measurements of the bone of interest for the particular system studied. Because of the large amount of growth that is likely to be calculated by the developed model (Equation (6)) in regions of directly applied forces, the maximum outer diameter and maximum cortical wall thickness reported in the literature is scaled by a 50% safety factor so that the simulation is not stopped prematurely because of these isolated regions of relatively large growth. The minimum inner diameter is used to prevent the interference of the elements on the inner surface should so much growth occur that the hollow cylinder become completely filled. A minimum inner diameter limit equal to the initial element size in the radial direction is selected in this model. Gradientless optimization methods, like the one used in this model, do not produce a unique solution, often preventing their use in comparative design studies. Instead, these methods are often employed to study specific design problems with a predefined required amount of strength improvement [85]. When used to understand the adaptive changes under given boundary conditions, these models typically are run to a predefined number of iterations [101]. Additionally, when applied to systems with a large number of design variables, the extreme stress states in localized regions may obscure the identification of the improved uniformity of the surface stress distribution. A method was developed to overcome these challenges through the creation of a quantifiable measure of the change in stress state uniformity of the design region in any system without the artificial influences from nodal extremes and stress concentrations [50]. Effective, rather than absolute, maximum and minimum measures of the local stress state in the design region are calculated as the average of the upper 25th percentile of nodal values and the average of the lower 25th percentile, respectively. The range of the stress variation is then defined as the difference between these averages. Calculated at each optimization iteration, this convergence measure is used to determine the amount of improvement in the uniformity of the design region stress state resulting from the shape optimization. In the current work, a 50% reduction from the conditions in the initially uniform cylindrical geometry indicates achievement of this optimization goal, providing a consistent means for measuring “convergence” and identifying “sufficiently” improved uniformity. The methods for stopping the simulation were implemented through the code written within an ANSYS script file that controlled the execution of the bone shape/strength adaptation simulation using the nodal strain energy densities calculated in the growth model FORTRAN subroutine and the coordinates of the inner and outer surface nodes taken from the ANSYS model database. In this model, the inner (endosteal) and outer (periosteal) surface optimizations function independently. The entire simulation is stopped only when both surfaces meet the convergence criteria or if any of the size limits are reached. 2.4. Model validation Two validation studies were performed to ensure that the function of the developed cortical bone shape adaptation model is consistent with the remodeling behavior observed in real bones. These studies followed previously published numerical simulations [46] of well-known animal studies of load-induced cortical bone adaptation in the radius bone that result from the osteotomy of the ulna bone of mature sheep [5] and juvenile swine [102]. In the current model simulating these experimental conditions, a three-dimensional initially circular hollow cylinder composed completely of cortical bone tissue, which represents diaphyseal region of the animal radius bone, was used. The bottom surface was completely constrained, while loads were applied to the top surface based on those described in the previously published simulation [46]. Trends in cross-sectional shape changes calculated by the current model at locations identified in the animal studies were compared with both the reported experimental observations [5, 102] and adaptations predicted by the previously published computational model [46]. Copyright © 2015 John Wiley & Sons, Ltd.


e02699 (14 of 29)

C. S. FLORIO

The offset loads on the radius bone induced under walking conditions because of the removal of the ulna bone in both validation studies showed regional periosteal cortical growth on the side of the bone facing the site of the removed ulna. This growth was centered on the line connecting the centroids of the radius and ulna bone prior to ulnar removal (Figure 5). Because of the differences in the relative orientations of the radial and ulnar centroids, the growth in the sheep model resulted in a diagonal shift in the bone shape, while the swine model caused a more symmetric growth pattern. This remodeling behavior predicted by the current model follows the trends in periosteal growth predicted by the previously published computational model [46] and the two experimental studies [5, 102]. Additionally, the current computational model captured the periosteal decay on the side of the bone opposite to that where growth occurred, a phenomenon which is observed in published images of both the experimentally adapted sheep and swine bones. Finally, the endosteal remodeling trends predicted by the current computational model corresponded to those in published images of the experimental studies. Specifically, in the sheep study, endosteal decay was noted on the side of the bone that had significant periosteal growth, and in the swine study, cortical wall thinning occurred on the side of the bone opposite that saw significant periosteal growth. The similarities in the cortical bone shape remodeling features predicted by the current model to those observed in experimental animal studies as well as calculated through previously published

Figure 5. Validation of developed shape adaptation model. Current model results at top. Initial geometry indicated by solid black line, and final geometry indicated by meshed area. (a) Simulation of offset loading on radius bone in experimental sheep ulnar osteotomy study of [5] and numerical model of [46] (b) Simulation of offset loading on radius bone in experimental swine ulnar osteotomy study of [102] and numerical model of [46]. Straight lines indicate those connecting the centroids of radius and ulna in the experimental studies. The lower images enclosed in boxes are directly reprinted with permission from [46], [5], and [102]. “Initial” and “Final” labels were added for clarity. Copyright © 2015 John Wiley & Sons, Ltd.



e02699 (15 of 29)

numerical models under common loading conditions, as shown through Figure 5, provide confidence in the ability of the developed model to predict realistic adaptive behavior. Such consistency of predicted trends in shape adaptation despite significant differences in numerical model formulation follows similar observations from previously published comparative studies [71, 73, 103]. The developed modeling components are capable of either functioning independently or as part of another bone shape adaptation simulation, structural shape optimization model, or examination of muscle activity. However, their coupled application can present a more inclusive representation of the mechanical influences of applied muscle force on bone strength adaptation. 2.5. Coupling of modeling components Because one of the most influential modeling features in the simulation of bone strength adaptation is an accurate depiction of the boundary condition, the complete analysis method developed in this work couples the determination of the individual muscle force magnitudes that, together, generate a net resultant force by a multibone, multimuscle system with the prediction of the ensuing changes to the shape and, consequently, strength of a bone of interest in the studied system. Only static isometric forces are considered, as they have been experimentally and clinically shown to induce the greatest changes in muscle and bone strength. The interactions between the bone segments are represented through the use of a built-in contact model, and the muscle forces are applied over areas representing the typical tendon cross-sectional areas for the particular musculoskeletal system of interest. The simulation is executed through an ANSYS script file that controls the input parameter values needed by each modeling component and directs the flow of the execution of each modeling component from the determination of the muscle force boundary conditions to the prediction of the amounts of local material accretion or resorption on the cortical bone surfaces to the nodal smoothing for mesh integrity, and finally, to the calculation of convergence or stopping measures. The script file also directs the retention of output data for subsequent analyses. The meshed bone geometry must be created externally to the developed program. Basic system parameters, such as stopping criteria, initial bone diameters, lengths of the design region, and muscle sizes and locations of attachment to the bone system (from which the force vector directions and joint moment arms are calculated) are input into the model through this script. The data obtained and stored throughout the optimization process include local and global measures of the stress state (stress/strain tensors, strain energy density, and principal stresses/strain), nodal positions, and measures of stress state uniformity. Modeling features were included to improve run efficiency and prevent extreme mesh distortion during the simulation. Because sorting of the large set of nodal data is often required for various operations in the developed subroutines, an efficient heap sorting method is employed [104]. Because the interior node smoothing routine, which is performed for each radial group of nodes in the meshed design region, is computationally intensive, it is carried out only once in every five optimization iterations. Surface node smoothing, however, is performed every optimization iteration. Additionally, limits are imposed on the change in nodal position per iteration, as a 20% change in the radial dimension of an element was found to create unacceptable distortion and potential element inversion. With the interior node smoothing performed every five iterations, the growth at each surface node per optimization iteration is uniquely limited to 4% of the pre-interior-nodesmoothed distance between the surface node and the next interior node with common local angular and axial coordinates. The model is not intended to predict exact amounts of material accretion and resorption that may be observed on the cortical surface of a bone because of mechanical environmental changes. The iterative application of forces in the model does not correspond to a particular amount of time. Instead, the model represents what would happen to the bone if a sufficient amount of time passed so that the adaptive shape changes resulted in a 50% improvement in the targeted bone’s surface stress state uniformity under the applied boundary conditions. Using the developed model, data that are difficult to obtain through physical experiments can be revealed. Trends over the entire optimization process or only comparing initial and final configurations can be found at any point, line or crosssection in the design region or throughout each adapted bone surface. Additionally, the effects of Copyright © 2015 John Wiley & Sons, Ltd.


e02699 (16 of 29)

C. S. FLORIO

one bone’s adaptations on the stresses in other bones and regions within the system studied can be examined. Used in parametric studies, the effects of different loading patterns can be directly compared. Such data can help identify conditions that invoke desired changes in targeted regions of the bone of interest. 3. MODEL IMPLEMENTATION – TEST CASE The effectiveness of the developed modeling method was examined through a representative study. The changes in the shape and stress distribution in a tibia bone were compared for the generation of various prescribed net resultant isometric forces by a leg system. 3.1. System representation The flexion/extension function of the leg was the focus of this study. Four bone segments: the pelvis, femur, tibia, and foot; three joints: the hip, knee, and ankle; and ten major muscles: the iliacus (IL), sartorius (SART), rectus femoris (RF), vastii (VAST), tibialis anterior (TA), soleus (SOL), gastrocnemius (GAST), short and long heads of the biceps femoris (SHBF and LHBF), and the gluteus maximus (GM) were included (Figure 6). The bone geometry was three-dimensional, symmetric about a sagittal mid-plane, and the exclusion of “out of plane” muscles simplified the determination of the muscle force magnitudes [105, 106]. The system geometry was based on a healthy, average-sized adult male. The boundary conditions consisted of the ten individual muscle forces and two displacement constraints. The muscle forces acted directly on the bone surfaces at their locations of attachment. The top of the pelvis was constrained in all degrees of freedom to represent the effect of the upper half of

Figure 6. System examined in representative parametric study. Copyright © 2015 John Wiley & Sons, Ltd.



e02699 (17 of 29)

the body. The end of the toe region was also fully constrained to represent its fixed contact with an immovable surface, where the net resultant force is generated. The interaction between the bone segments at the joints was simulated through an augmented Lagrange multiplier contact model [87] that allowed compliance but not relative motion, consistent with the isometric conditions studied [45]. The leg was arranged with the hip, knee, and ankle joints in their neutral positions. The shape strength adaptations in the tibia were examined. This bone was selected because disuse based bone strength losses have been found to be greater at distal locations [19, 107–109], and, therefore, targeted bone strengthening techniques for the tibia may be especially beneficial. The model of this bone system was generated using commercially-available three-dimensional computer aided design (CAD) software (Pro/Engineer Wildfire 4.0, PTC, Needham, MA, USA) [110]. 3.1.1. Muscles. The muscles were not physically modeled, but their effects were included through individual muscle forces acting on the bone system studied. Physical attributes of the included muscles, their physiological cross-sectional areas, and their locations of attachment to the bones, were obtained from a survey of the published literature. The physiological cross-sectional area of each muscle is employed as a weighting factor in the numerical optimization model used to determine the magnitudes of the forces generated by the included muscles that create the desired net resultant force. Because the value of this parameter depends on the size and strength of the subject from which the measurement was taken, significant variation has been reported. The values used in this work were an average of those reported in nine published studies, excluding minimum and maximum reported measurements, [51] and are listed in Table I. The attachment locations were taken directly from published physical measurements [111] identified based on physically significant and uniquely defined local coordinate systems [111]. These local coordinate systems were recreated in the bone system models developed in the current work, and the locations of muscle force applications were subsequently defined through a series of coordinate transformations. Because of the simplified bone geometries used in the current model, when necessary, the published locations of muscle attachment were projected onto the sagittal plane of symmetry. The muscles were assumed to act along straight lines connecting their attachment points. The moment arms necessary to determine the individual muscle force magnitudes and the directions in which each muscle force acted were found using these locations and the developed CAD model of the bone system. Details about the calculated values of the muscle attachment locations, force direction vectors, and moment arms can be found elsewhere [51]. The muscle forces acted over regions representing the typical cross-sectional area of a human leg muscle tendon. The tendon area used in this work had a width of 9 mm, was the same for each muscle, and was calculated as the average of four physical measurements reported in the literature [45]. The muscle force was distributed in this area so that it was greatest at the central nodes and transitioned linearly to the unloaded region surrounding the attachment area, reducing the resulting stress gradients.

Table I. Muscle physiological cross-sectional areas used (cm2). Muscle

PCSA

IL SART RF VAST TA SOL GAST SHBF LHBF

23.0 3.1 31.8 78.1 18.1 97.8 39.2 9.9 27.9



e02699 (18 of 29)

C. S. FLORIO

3.1.2. Bones. Similar to the locations of the muscle attachment, the bone geometry in this model was derived from reported physical measurements [111, 112], projected, when necessary, onto the sagittal plane of symmetry. The femur and tibial were simplified as uniform circular cylinders. The foot bones were simplified into a unified geometry. Because only isometric conditions with fixed joint angles were considered, the joint geometries were reduced to basic representations of the articular surfaces. Cortical bone tissue was modeled as transversely isotropic in the femur and tibia, with the properties different in the axial direction [113] (Table II), and fully isotropic in the pelvis and foot, with an elastic modulus of 16.35 GPa [114] and a Poisson’s ratio of 0.34 [115]. Fully isotropic cancellous bone tissue was included with a Poisson’s ratio of 0.3 and an elastic modulus of 11 GPa, the midpoint of the typically reported range of 9–15 GPa [116], as the specific modulus value used has been shown to have little effect on the resulting stress distributions [45]. Cancellous bone tissue completely filled the cortical shell in the pelvis. The forefoot was solid cortical bone tissue, while the region from the calcaneus to the navicular tubercle was filled with cancellous bone tissue [111]. The femur and the tibia were hollow in the central, diaphyseal region and had cancellous bone tissue filling end regions each comprising 12% of the total bone length [1, 117]. Shape adaptation in the tibia was simulated in the entire hollow region on the inner (endosteal) surface and on the initially uniform cylindrical region (up to the joints) on the outer (periosteal) surface. The initial outer and inner diameters of both the femur and tibia were 30 mm and 15 mm, respectively [118]. The size limits in the shape adaptation simulation were derived from a maximum reported outer diameter of 35.5 mm [119] and cortical thickness of 12.5 mm [120, 121]. The minimum allowable cortical thickness was 2 mm [122]. The CAD-developed three-dimensional bone geometries were imported into the finite element software (ANSYS Version 12.1, ANSYS Inc., Cannonsburg, PA, USA) [87] so that each of the bone segments remained independent with unique articular surfaces. In contrast, cortical and cancellous tissue regions were bound through shared surfaces. 3.1.3. Mesh. The geometry was meshed using a combination of 8-node hexahedral elements and 20-node tetrahedral elements, with pyramidal elements transitioning between regions of each type of element. The element size was limited by the assumption that the bone material could be represented as a continuum. For such a continuum assumption, the minimum element size must be at least an order of magnitude greater than the characteristic size of the bone material [123]. With the focus of the study on cortical bone adaptation, the characteristic material size was chosen as that of an osteon, 100–300 μm [116], and thus, the minimum element size in this model was 1 mm. The foot, pelvis, and joint regions of the long bones were discretized with tetrahedral elements using the built-in automated mesh generation function and refined as necessary near the applied boundary conditions, in the muscle attachment areas, and on the joint surfaces. The hexahedral elements were manually created as a mapped mesh in the cylindrical regions of the femur and tibia bone, except for a small central core of tetrahedral elements required to fill the solid cancellous region. The hexahedral mapped mesh was necessary for the proper function of the developed shape adaptation model. Nodes were evenly spaced in the radial direction in the cortical regions, decreasing Table II. Transversely isotropic material properties for the cortical bone in femur and tibia (all in GPa except ν). Material property

Value

Ex Ey Ez Gxy Gyz Gxz ν

9.55 9.55 16.61 3.28 4.74 4.74 0.37




e02699 (19 of 29)

in size moving to the center of the cancellous bone tissue filled the interior volumes. In the angular coordinate direction, nodes were spaced closer together anteriorly and posteriorly, near muscle force application regions, transitioning to greater medial and lateral spacing. The axial direction featured similar nodal spacing as the angular direction in the regions of muscle force application, creating nearly cubic elements, with gradual transitions to larger spacing away from the boundary conditions. Eightnode elements were used because midside nodes are not accessible from the ANSYS database and, therefore, cannot be updated during the shape optimization to prevent mesh distortion. The final mesh consisted of 184,067 elements and 242,066 nodes. It was disassociated from the underlining solid model so that the surface nodes could be freely moved during the shape adaptation simulations. 3.2. Representative parametric study of effect of resultant force direction A parametric study was performed on the effect of varying the direction of the prescribed isometric resultant force generated on the predicted tibial shape adaptations. The alteration of this loading parameter was found to produce wide variations in the bone stresses [45] and could be easily implemented to target the strengthening of a specific fracture prone bone region. In this study, a 50 N resultant force, directed either anteriorly, inferiorly, posteriorly, or superiorly, was generated at the toe by a straight leg against a fixed surface. A 64-bit LINUX RedHat machine with dual quad core 2.40 GHz processors and 12 GB RAM was used. Run times were on the order of one week. 3.3. Typical resulting data Because the developed computational model determines absolute muscle force magnitudes and bone stresses throughout the three-dimensional multisegment system studied and this data, as well as the nodal coordinates, can be recorded at every optimization iteration, a more complete analysis of the mechanically based bone strength shape optimization process can result than may be achieved at that discrete locations and times that are typically measured in experimental studies. Local and global data can be compared for the conditions studied at a particular node, along a curve, at a cross-section, or on a surface both in the design region and in neighboring joints and bone segments. Additionally, this data can be used in subsequent calculations of mechanical, statistical, or geometric measures. The data presented is representative, but not inclusive, of comparisons that may be gained from such a parametric study. The applied muscle force magnitudes can be directly compared either in absolute values or in relative intensities. Because the strength shape adaptations of cortical bone are driven by the variation of the stress state over the bone’s surface, normalizing the muscle force magnitudes for a given set of conditions by the greatest muscle force magnitude calculated for that set can reveal the dominant active muscles and the relative contributions of the others (Table III). For example, while the same set of muscles is active to create the anteriorly directed and superiorly directed resultant forces (and, similarly, the interiorly directed and posteriorly directed resultants), the relative intensities of the forces generated by these muscles vary. Anteriorly attached muscles drive the anteriorly directed Table III. Relative muscle force intensities. Muscle

Attach

Ant.

Inf.

Post.

Sup.

IL SART RF VAST TA SOL GAST SHBF LHBF GM

Anterior Anterior Anterior Anterior Anterior Posterior Posterior Posterior Posterior Posterior

0.003 0.010 1.000 0.133 0.037 0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.240 0.000 1.000 0.092 0.866 0.113

0.000 0.000 0.000 0.000 0.368 0.000 0.605 0.095 1.000 0.556

0.003 0.009 0.754 1.000 0.380 0.000 0.000 0.000 0.000 0.000



e02699 (20 of 29)

C. S. FLORIO

and superiorly directed resultant force cases, with the knee flexor rectus femoris providing the majority of the anteriorly directed force, while a more distributed effort of the rectus femoris and vastii and, to a lesser extent, the tibialis anterior muscles generates the superiorly directed force. Similarly, the inferiorly directed and posteriorly directed resultant forces are both created by a common set of muscles. Yet, the gastrocnemius dominates the interiorly directed resultant, while the posteriorly directed resultant requires more significant activity of the hip extensors, namely the long head of the biceps femoris and the gluteus maximus. When analyzed in combination with the relative orientations of the muscle force vectors, constrained surfaces, included joints, and resultant force vectors in the straight leg arrangement of the system studied, the feasibility of the model-predicted muscle activity is realized. This type of analysis of the muscle force data generated by the current model can help guide the development of exercise regimens to target particular muscles or bone regions for strengthening and can have implications on the effects of the induced muscle activity on strength changes in other bone regions near less dominant muscle forces. Comparisons of the stress distributions induced in the bone system by these muscle forces and the ensuing bone shape adaptations can complement the study of the muscle force intensities to reveal local and global system effects. Stress distributions throughout the multi-bone system studied present differences in its response to the mechanical loads applied by the compared conditions and the changes that occur because of cortical remodeling of the targeted bone (Figure 7). Local shape changes resulting from the simulated accretion and removal of material on the cortical bone surfaces

Figure 7. von Mises stress distributions (MPa) for each of the load cases studied before and after tibial bone shape adaptation simulation. Red arrows indicate active muscle forces. Dashed line indicates location of further study in Figure 9. Copyright © 2015 John Wiley & Sons, Ltd.



e02699 (21 of 29)

can be qualitatively compared (Figure 8) or quantified by the change in surface node coordinates or the change in cross-sectional geometric measures (Figure 9). While standard measures of crosssectional geometry can provide insight into changes in the ability of the analyzed bone region to resist particular basic loading modes, such as bending and torsional stresses, the effectiveness of the shape adaptations to improve the bone’s response to the particular mechanical conditions to which it has been optimized (each load case studied) is better revealed by comparisons of local von Mises stresses (Figure 10). Additionally, the progress of each simulation towards the defined optimization goal, remodeling-induced improved strain energy density uniformity, can be indicated by changes in local nodal values, such as along axial lines (Figure 11), or by global measures, such as averages over each surface being optimized (Figure 12). Comparisons of the amounts of these changes can, therefore, be made between the various loading cases studied to indicate their relative bone strengthening effectiveness. Each of these analysis methods may be used alone or in combination to gain insight into the conditions studied. Common changes in the overall tibial bone profiles in the representative study of this work (Figure 8) follow the similarities in trends of muscle force activity for the posteriorly directed and inferiorly directed and for the anteriorly directed and superiorly directed resultant forces (Table III). However, more subtle features in the resulting shape adaptations are seen which may not be directly predicted from study of the muscle force magnitudes alone. For example, there is a common region of increased periosteal growth on the posterior proximal side of the bone at a location about two-thirds of the tibial length from the ankle for the anteriorly directed and posteriorly directed resultant forces despite significant differences in muscle activities. The developed model does predict large amounts of growth at the locations of concentrated stress because of the direct application of the muscle forces. This may, in part, be the result of the simplified initially uniform circular cylindrical bone geometry. Modeling a more realistic initial geometry, with widened cross-sections near the ends of the bones where many muscles attach, may alleviate the need for this large growth in order to achieve the desired improvement in stress state uniformity. Nonetheless, in regions away from the effects of the directly applied

Figure 8. Comparison of shape changes on periosteal surface for the conditions studied.

Figure 9. Comparison of geometric measures at the particular cross-section indicated by the dashed line in Figure 7, normal to the longitudinal axis, and far from the applied boundary conditions. Copyright © 2015 John Wiley & Sons, Ltd.


e02699 (22 of 29)

C. S. FLORIO

Figure 10. Comparison of change in von Mises stress (MPa) contours at the same cross-section as Figure 9 due to shape adaptation. Note different scales for each load case.

Figure 11. Inferiorly directed resultant force: change in nodal strain energy density along a line on the posterior side of the bone (x coordinate) from the ankle to the knee (a) outer (periosteal) surface and (b) inner (endosteal) surface demonstrating improved uniformity.

muscle forces, the more subtle changes predicted by the current model indicate likely representations of the growth variations that may be induced by the exercise conditions studied. Relative comparisons of the von Mises stress distributions help to understand the effects that the direction of the resultant force generated by the leg system have on the system studied even before tibial bone adaptations occur (Figure 7) and demonstrate the wide range of conditions that may be Copyright © 2015 John Wiley & Sons, Ltd.



e02699 (23 of 29)

Figure 12. Comparison of change in the average and standard deviation of the elemental strain energy density on (a) outer (periosteal) and (b) inner (endosteal) surfaces. Note the wide range of conditions analyzed with the developed technique and significant decrease in variation due to the shape adaptation modeled even as the average increases under some conditions.

evaluated with the developed model. Local surface bone stresses vary several orders of magnitude because of the alteration of the direction, not the magnitude, of the resultant force. This model data can help identify the bone’s initial susceptibility to fracture under each of the load conditions examined. The effects of the adaptive shape changes on the local tibial stresses as well as on the bone stresses at other system locations are also determined by the model (Figure 8). For example, the alterations in the tibial bone shape when adapting to the superiorly directed resultant force successfully decrease the proximal anterior tibial stresses from those in the initial geometry. However, the adaptive changes in the tibia also result in increased stresses in the anterior proximal femur near the hip joint. The identification of such secondary effects can prompt further analysis of the local femoral stress values to indicate if fracture at this location is of increased concern under these conditions induced by the tibial remodeling. Comparisons of local adaptive changes can be quantified. For example, cross-sectional geometric measures at a location one-third of the tibial length from the ankle (Figure 9) show model-predicted increases in area and moments of inertia only under the inferiorly directed resultant force. The bones adapted to the anteriorly directed and superiorly directed resultants have smaller moments of inertia and area at this axial location than the initially circular cylinder despite little change in the maximum radial measurement (Figure 9). The effects of these geometric modifications under each of the studied conditions are revealed through the nodal von Mises stress (Figure 10), a measure typically used in mechanical component failure analyses. When evaluated in combination with trends in changes in nodal coordinates on both the periosteal and endosteal surfaces, correlations can be found between local shape adaptations and local stresses alterations. For example, at the studied cross-section, significant medial and lateral thinning result in increased anterior and posterior stresses under the anteriorly directed and superiorly directed resultant forces. However, the posteriorly directed and inferiorly directed resultant forces create adaptive changes that improve stress state uniformity at this cross-section under the respective loading conditions. To create these nearly Copyright © 2015 John Wiley & Sons, Ltd.


e02699 (24 of 29)

C. S. FLORIO

uniform-stress regions, the inferiorly directed resultant force causes varying amounts of stress reductions throughout the entire cross-section while the posteriorly directed resultant force increases the local anterior stresses and decreases the local posterior stresses. In both of these loading cases, the local stress changes are a result of endosteal and periosteal growth, especially on the posterior side of the tibia, helping the studied cross-section of the tibia bone to better resist the applied mechanical loads than could the initially circular geometry. Improved strain energy density uniformity is found locally (Figure 11) and globally (Figure 12), indicating successful achievement of the optimization goal. Typical changes in nodal strain energy density along an axially directed line on the periosteal and endosteal bone surfaces as a result of the shape adaptations show how the predicted growth patterns significantly reduce the large stresses induced by the directly applied muscle forces (large spike in Figure 11(a)) and how the strain energy density in the optimized geometry becomes nearly uniform, aside from effects of local constraints and material property changes. All of the load cases studied revealed reduced surface stress state variation (Figure 12), though only the posteriorly directed and inferiorly directed resultants also caused overall reduced average surface stresses. Interestingly, the significantly larger magnitudes of the average strain energy density and surface variation under the posteriorly directed resultant force induced similar shape adaptations to that of the less intense stresses caused by the inferiorly directed resultant force. Similar observations are found for the superiorly directed and anteriorly directed resultants. Such data derived from the developed computational model may invoke further computational, experimental, or clinical studies to find exercises that produce desired adaptive strength changes with reduced risk of exercise-induced fracture.

4. CONCLUSIONS A complete numerical model to simulate the mechanical factors that influence bone shape/strength adaptation is developed in this work, and the procedure for its implementation is presented. The model uses a multibone/multimuscle representation of the system investigated. It calculates the individual muscle forces developed in the system to generate a prescribed interaction with its external environment (such as a force between the system and a fixed surface) and the ensuing local amounts of material accretion and resorption that reduce the variation of the stress state within the bone of interest. The model has been shown to predict bone adaptation trends that are consistent with experimentally observed cortical bone surface remodeling behavior. The representative implementation of the developed techniques presented in this work demonstrates the model’s effectiveness in quantitatively comparing the local and global effects of significantly different sets of boundary conditions on the bone strength adaptation in the musculoskeletal system studied. More detailed, complex, or even different systems can be easily studied using the developed technique. Similarly, more detailed bone geometry, joint contact models, material properties, or representations of muscle-bone attachment can be applied, as these model features are created independently of the developed methods. Because of its strong link to improving bone strength, the focus of the current model development and function was on isometric bone strengthening regimens. With its implementation through a series of user-developed subroutines, other conditions can readily be evaluated by the presented modeling technique through the development of appropriate muscle force optimization or structural analysis methods. Further model refinements could be added in the future to more directly address the effects of gravitational forces, the general reduction, rather than the alteration, of muscle-induced forces, and the study of dynamic loading conditions. It should be emphasized that the model is intended as a tool for comparative studies of the relative effects of modified loading conditions on local bone strength, not for replication of quantitative experimental measurements. Because it is based on a gradientless optimization method, an absolute global optimal shape that minimizes stress variations over the bone surface is not found. Instead, it is assumed that the “improved” conditions found by the model are not radically different than the initial starting point, both geometrically and mechanically. The model provides a means for the consistent definition of such an improvement so that direct comparisons can be made among varying exercise-based conditions. The benefit of the developed model is the ability to study a range Copyright © 2015 John Wiley & Sons, Ltd.



e02699 (25 of 29)

of alterations from a given mechanical state without concern of the effects of selected model parameters so that quantitative relative assessments of the local effects of mechanical environmental changes on cortical bone shape and, therefore, strength can be realized. Functioning within this scope, the tools developed in this work can be useful in the analysis and design of mechanically based bone strengthening techniques. They can be used to target particular fracture-prone regions of bone for strengthening or to understand the effects of particular exercises or activities on both increases and decreases in local bone strength, as areas that are most significantly affected by changes to system loading may be readily identified. Based on methods often employed to study inert mechanical components, this computational modeling method provides capabilities that are often lacking in many current experimentally based investigations of the relationships between alterations in a bone’s mechanical environment and its strength adaptations. ACKNOWLEDGEMENTS

This work was partially supported by Amelia Earhart Fellowship from the Zonta International Foundation. REFERENCES 1. Uhthoff HK, Jaworski ZFG. Bone loss in response to long-term immobilisation. Journal of Bone and Joint Surgery - Series B 1978; 60B:420–429. 2. Be’ery-Lipperman M, Gefen A. Contribution of muscular weakness to osteoporosis: computational and animal models. Clinical Biomechanics 2005; 20(9):984–997. 3. Liskova M, Hert J. Reaction of bone to mechanical stimuli. Part 2. Periosteal and endosteal reaction of tibial diaphysis in rabbit to intermittent loading. Folia Morphologica 1971; 29(3):301–317. 4. Hert J, Liskova M, Landa J. Reaction of bone to mechanical stimuli. Part 1. Continuous and intermittent loading of tibia in rabbit. Folia Morphologica 1971; 19(3):290–371. 5. Lanyon LE, Goodship AE, Pye CJ et al. Mechanically adaptive bone remodelling. Journal of Biomechanics 1982; 15(3):141–154. 6. Lanyon LE, Rubin CT. Static vs dynamic loads as an influence on bone remodelling. Journal of Biomechanics 1984; 17(12):897–905. 7. Borer KT, Kuhns LR. Radiographic evidence for acceleration of skeletal growth in adult hamsters by exercise. Growth 1977; 41(1):1–13. 8. Steinberg ME, Trueta J. Effects of activity on bone growth and development in the rat. Clinical Orthopaedics and Related Research 1981; 156:52–60. 9. Woo SL-Y, Kuei SC, Amiel D et al. The effect of prolonged physical training on the properties of long bone: a study of Wolff’s law. Journal of Bone and Joint Surgery - Series A 1981; 63A(5):780–787. 10. Warner SE, Shea JE, Miller SC et al. Adaptations in cortical and trabecular bone in response to mechanical loading with and without weight bearing. Calcified Tissue International 2006; 79(6):395–403. 11. Lanyon LE. The measurement of bone strain ‘in vivo’. Acta Orthopaedica Belgica 1976; 42(sup1):98–108. 12. Lanyon LE, Goodship AE, Baggott DG. The significance of bone strain ‘in vivo’. Acta Orthopaedica Belgica 1976; 42(sup1):109–122. 13. Judex S, Gross TS, Bray RC et al. Adaptation of bone to physiological stimuli. Journal of Biomechanics 1997; 30(5):421–429. 14. Adams DJ, Spirt AA, Brown TD et al. Testing the daily stress stimulus theory of bone adaptation with natural and experimentally controlled strain histories. Journal of Biomechanics 1997; 30(7):671–678. 15. Gross TS, Edwards JL, Mcleod KJ et al. Strain gradients correlate with sites of periosteal bone formation. Journal of Bone and Mineral Research 1997; 12(6):982–988. 16. Skedros JG, Mason MW, Nelson MC et al. Evidence of structural and material adaptation to specific strain features in cortical bone. Anatomical Record 1996; 246(1):47–63. 17. Deitrick JE, Whedon GD, Shorr E. Effects of immobilization upon various metabolic and physiologic functions of normal men. The American Journal of Medicine 1948; 4(1):3–36. 18. Greenleaf JE, Bulbulian R, Bernauer EM et al. Exercise-training protocols for astronauts in microgravity. Journal of Applied Physiology 1989; 67(6):2191–2204. 19. Whalen R. Musculoskeletal adaptation to mechanical forces on earth and in space. Physiologist 1993; 36(1 Suppl): S127–S130. 20. Pearson OM, Lieberman DE. The aging of Wolff’s "law": Ontogeny and responses to mechanical loading in cortical bone. Yearbook of Physical Anthropology 2004; 47:68–99. 21. Engelke K, Kemmler W, Lauber D et al. Exercise maintains bone density at spine and hip. EFOPS: a 3-year longitudinal study in early postmenopausal women. Osteoporosis International 2006; 17(1):133–142. Copyright © 2015 John Wiley & Sons, Ltd.


e02699 (26 of 29)

C. S. FLORIO

22. Vainionpaa A, Korpelainen R, Sievanen H et al. Effect of impact exercise and its intensity on bone geometry at weight-bearing tibia and femur. Bone 2007; 40:604–611. 23. Guadalupe-Grau A, Fuentes T, Guerra B et al. Exercise and bone mass in adults. Sports Medicine 2009; 39(6): 439–468. 24. Nikander R, Kannus P, Dastidar P et al. Targeted exercises against hip fragility. Osteoporosis International 2009; 20(8):1321–1328. 25. Crossley K, Bennell KL, Wrigley T et al. Ground reaction forces, bone characteristics, and tibial stress fracture in male runners. Medicine and Science in Sports and Exercise 1999; 31(8):1088–1093. 26. Shackelford LC, Leblanc AD, Driscoll TB et al. Resistance exercise as a countermeasure to disuse-induced bone loss. Journal of Applied Physiology 2004; 97(1):119–129. 27. Smith SM, Heer MA, Shackelford LC et al. Benefits for bone from resistance exercise and nutrition in longduration spaceflight: evidence from biochemistry and densitometry. Journal of Bone and Mineral Research 2012; 27(9):1896–1906. 28. Burr DB, Milgrom C, Fyhrie D et al. In vivo measurement of human tibial strains during vigorous activity. Bone 1996; 18(5):405–410. 29. Ruff C, Holt B, Trinkaus E. Who’s afraid of the big bad Wolff?: “Wolff’s law” and bone functional adaptation. American Journal of Physical Anthropology 2006; 129:484–498. 30. Mittlmeier T, Mattheck C, Dietrich F. Effects of mechanical loading on the profile of human femoral diaphyseal geometry. Medical Engineering and Physics 1994; 16(1):75–81. 31. Fischer KJ, Jacobs CR, Carter DR. Computational method for determination of bone and joint loads using bone density distributions. Journal of Biomechanics 1995; 28(9):1127–1135. 32. Jang IG, Kim IY. Computational simulation of simultaneous cortical and trabecular bone change in human proximal femur during bone remodeling. Journal of Biomechanics 2010; 43(2):294–301. 33. Bitsakos C, Kerner J, Fisher I et al. The effect of muscle loading on the simulation of bone remodelling in the proximal femur. Journal of Biomechanics 2005; 38(1):133–139. 34. Speirs AD, Heller MO, Duda GN et al. Physiologically based boundary conditions in finite element modelling. Journal of Biomechanics 2007; 40(10):2318–2323. 35. Kumar NC, Dantzig JA, Jaisuk IM et al. Numerical modeling of long bone adaptation due to mechanical loading: correlation with experiments. Annals of Biomedical Engineering 2010; 38(3):594–604. 36. Cowin SC, Firoozbakhsh K. Bone remodeling of diaphysial surfaces under constant load: theoretical predictions. Journal of Biomechanics 1981; 7:471–484. 37. Cowin SC. Bone remodeling of diaphyseal surfaces by torsional loads: theoretical predictions. Journal of Biomechanics 1987; 20(11–12):1111–1120. 38. Van Der Meulen MCH, Beaupre GS, Carter DR. Mechanobiologic influences in long bone cross-sectional growth. Bone 1993; 14:635–342. 39. Levenston ME, Beaupre GS, Carter DR. Loading mode interactions in simulations of long bone cross-sectional adaptation. Computer Methods in Biomechanics and Biomedical Engineering 1998; 1(4):303–319. 40. Fischer KJ, Jacobs CR, Levenston ME et al. Different loads can produce similar bone density distributions. Bone 1996; 19(2):127–135. 41. Fernandez J, Sartori M, Lloyd D et al. Bone remodelling in the natural acetabulum is influenced by muscle forceinduced bone stress. International Journal for Numerical Methods in Biomedical Engineering 2014; 30(1): 28–41. 42. Viceconti M, Testi D, Taddei F et al. Biomechanics modeling of the musculoskeletal apparatus: status and key issues. Proceedings of the IEEE 2006; 94(4):725–738. 43. Ezquerro F, Simon A, Prado M et al. Combination of finite element modeling and optimization for the study of lumbar spine biomechanics considering the 3D thorax-pelvis orientation. Medical Engineering and Physics 2004; 26(1):11–22. 44. Marom SA, Linden MJ. Computer aided stress analysis of long bones utilizing computed tomography. Journal of Biomechanics 1990; 23(5):399–404. 45. Florio CS, Narh KA. Development of a modeling technique for the investigation of muscle activity and its effect on bone stresses in the human leg during an isometric exercise. Simulation 2011; 87(4):313–333. 46. Cowin SC, Hart RT, Balser JR et al. Functional adaptation in long bones: establishing in vivo values for surface remodeling rate coefficients. Journal of Biomechanics 1985; 18(9):665–671. 47. Fridez P, Rakotomanana L, Terrier A et al. Three dimensional model of bone external adaptation. Computer Methods in Biomechanics and Biomedical Engineering 1998; 2:189–196. 48. Frost HM. A 2003 update of bone physiology and Wolff’s law for clinicians. Angle Orthodontist 2004; 74(1):3–15. 49. Ferretti JL, Frost HM, Gasser JA et al. Perspectives on osteoporosis research: its focus and some insights from a new paradigm. Calcified Tissue International 1995; 57(6):399–404. 50. Florio CS. Development of a widely applicable gradientless shape optimization based bone adaptation model for comparative parametric studies. Structural and Multidisciplinary Optimization 2015. doi:10.1007/s00158-015-1227-y 51. Florio CS. Computational optimization methods for modeling the effect of muscle forces on bone strength adaptation. Doctoral Dissertation, New Jersey Institute of Technology, Newark, NJ, 2014. 52. Crowninshield RD, Brand RA. A physiologically based criterion of muscle force prediction in locomotion. Journal of Biomechanics 1981; 14(11):793–801.




e02699 (27 of 29)

53. Arjmand N, Shirazi-Adl A. Sensitivity of kinematics-based model predictions to optimization criteria in static lifting tasks. Medical Engineering & Physics 2006; 28(6):504–514. 54. Van Bolhuis BM, Gielen CCM. A comparison of models explaining muscle activation patterns for isometric contractions. Biological Cybernetics 1999; 81(3):249–261. 55. Dul J, Townsend MA, Shiavi R et al. Muscular synergism—I. On criteria for load sharing between synergistic muscles. Journal of Biomechanics 1984; 17(9):663–673. 56. Tsirakos D, Baltzopoulos V, Bartlett R. Inverse optimization: functional and physiological considerations related to the force-sharing problem. Critical Reviews in Biomedical Engineering 1997; 25(4–5):371–407. 57. Rosen JB. The gradient projection method for nonlinear programming Part 1. Linear constraints. Journal of the Society for Industrial and Applied Mathematics 1960; 8(1):181–217. 58. Fox RL. Optimization Methods for Engineering Design. Addison-Wesley Publishing Company: Reading, MA, 1971. 59. Haftka RT, Gurdal Z, Kamat MP. Elements of Structural Optimization. Kluwer Academic Publishers: Boston, MA, 1990. 60. Haug EJ, Arora JS. Applied Optimal Design: Mechanical Structures and Systems. John Wiley and Sons: New York, 1979. 61. Gill PE, Murray W, Wright MH. Practical Optimization. Elsevier Academic Press: New York, NY, 1986. 62. Mathworks Inc. MATLAB (2009a edn). Mathworks Inc: Natick, MA, USA, 2009. 63. Zienkiewicz OC, Campbell JS. Shape optimization and sequential linear programming. In Optimum Structural Design: Theory and Applications, Gallagher RH, Zienkiewicz OC (eds.). Wiley: New York, NY, 1973. 64. Ding Y. Shape optimization of structures: a literature survey. Computers and Structures 1986; 24(6): 985–1004. 65. Van Keulen F, Haftka RT, Kim NH. Review of options for structural design sensitivity analysis. Part 1: linear systems. Computer Methods in Applied Mechanics and Engineering 2005; 194:3213–3243. 66. Kolda TG, Lewis RM, Torczon V. Optimization by direct search: new perspectives on some classical and modern methods. SIAM Review 2003; 45(3):385–482. 67. Das R, Jones R, Peng D. Optimization of damage tolerant structures using a 3D biological algorithm. Engineering Failure Analysis 2006; 13:362–379. 68. Kummer BKF. Biomechanics of bone: mechanical properties, functional structure, and functional adaptation. In Biomechanics: Its Foundations and Objectives, Fung YC, Perrone N, Anliker M (eds.). Prentice Hall: Englewood Cliffs, NJ, 1972. 69. Cowin SC, Van Buskirk WC. Surface bone remodeling induced by a medullary pin. Journal of Biomechanics 1979; 12(4):269–276. 70. Huiskes R, Weinans H, Grootenboer HJ et al. Adaptive bone-remodeling theory applied to prosthetic-design analysis. Journal of Biomechanics 1987; 20(11–12):1135–1150. 71. Brown TD, Pedersen DR, Gray ML et al. Toward an identification of mechanical parameters initiating periosteal remodeling: a combined experimental and analytic approach. Journal of Biomechanics 1990; 23(9):893–897. 72. Xu W, Robinson K. X-ray image review of the bone remodeling around an osseointegrated trans-femoral implant and a finite element simulation case study. Annals of Biomedical Engineering 2008; 36(3):435–443. 73. Florio CS, Narh KA. Effect of modeling method on prediction of cortical bone strength adaptation under various loading conditions. Meccanica 2013; 48(2):393–413. 74. Morris CE. Mechanosensitive ion channels. Journal of Membrane Biology 1990; 113(2):93–107. 75. Ko KS, Mcculloch CG. Partners in protection: interdependence of cytoskeleton and plasma membrane in adaptations to applied forces. Journal of Membrane Biology 2000; 174(2):85–95. 76. Terrier A, Rakotomanana L, Ramaniraka AN et al. Adaptation models of anisotropic bone. Computer Methods in Biomechanics and Biomedical Engineering 1997; 1(1):47–59. 77. Wu Z. An efficient approach for shape optimization of components. International Journal of Mechanical Sciences 2005; 47(10):1595–1610. 78. Le C, Bruns T, Tortorelli D. A gradient-based parameter-free approach to shape optimization. Computational Methods in Applied Mechanics and Engineering 2011; 200:985–996. 79. Schnack E. Optimization procedure for stress concentrations by the finite element technique. International Journal for Numerical Methods in Engineering 1979; 14(1):115–124. 80. Tekkaya AE, Guneri A. Shape optimization with the biological growth method: a parameter study. Engineering Computations (Swansea, Wales) 1996; 13(8):4–18. 81. Meske R, Sauter J, Schnack E. Nonparametric gradient-less shape optimization for real-world applications. Structural and Multidisciplinary Optimization 2005; 30(3):201–218. 82. Frost HM. Bone’s mechanostat: a 2003 update. Anatomical Record - Part A Discoveries in Molecular, Cellular, and Evolutionary Biology 2003; 275(2):1081–1101. 83. Heller M, Kaye R, Rose LRF. Gradientless finite element procedure for shape optimization. Journal of Strain Analysis for Engineering Design 1999; 34(5):323–336. 84. Yang RJ, Chen CJ. Stress-based topology optimization. Structural Optimization 1996; 12:98–105. 85. Umetani Y, Hirai S. An adaptive shape optimization method for structural material using the growing-reforming procedure. In The 1975 Joint JSME-ASME Applied Mechanics Western Conference. The Japan Society of Mechanical Engineers: Ilikai Hotel, Honolulu, Hawaii, 1975; 359–365.



e02699 (28 of 29)

C. S. FLORIO

86. Bedair O. The application of indirect boundary element method to optimum shape design. Computers and Structures 1997; 65(5):651–668. 87. ANSYS Inc. ANSYS. In (12.1 edn). ANSYS, Inc.: Canonsburg, PA, USA. 88. Florio PJ. Eg-202: Computer Graphics Package. NJIT Bookstore: Newark, NJ, 1983. 89. Imam MH. Three-dimensional shape optimization. International Journal for Numerical Methods in Engineering 1982; 18(5):661–673. 90. Zhang S, Belegundu AD. Mesh distortion control in shape optimization. AIAA Journal 1993; 31(7):1360–1362. 91. Banichuk NV, Barthold FJ, Falk A et al. Mesh refinement for shape optimization. Structural Optimization 1995; 9(1):46–51. 92. Wilke DN, Kok S, Groenwold AA. A quadratically convergent unstructured remeshing strategy for shape optimization. International Journal for Numerical Methods in Engineering 2006; 65:1–17. 93. Qian Y, Liu Z, Fan Y. Numerical simulation of canine bodily movement. International Journal for Numerical Methods in Biomedical Engineering 2010; 26(2):157–163. 94. Luo Y. R-adaptation algorithm guided by gradients of strain energy density. International Journal for Numerical Methods in Biomedical Engineering 2010; 26(8):1077–1086. 95. Chang K-H, Choi KK. An error analysis and mesh adapation method for shape design of structural components. Computers and Structures 1992; 44(6):1275–1289. 96. Batina JT. Unsteady Euler airfoil solutions using unstructured dynamic meshes. AIAA Journal 1990; 28(8): 1381–1388. 97. George PL. Automatic Mesh Generation: Application to Finite Element Methods. John Wiley and Sons: New York, NY, 1991. 98. Salagame RR, Belegundu AD. Distortion, degeneracy and rezoning in finite elements – a survey. Sadhana 1994; 19(2):311–335. 99. Peyret R, Taylor TD. Computational methods for fluid flow (vol). New York. Springer-Verlang: New York, NY, 1983. 100. Ugural AC, Fenster SK. Advanced Strength and Applied Elasticity (3rd edn). Prentice Hall PTR: Upper Saddle River, NJ, 1995. 101. Beaupre GS, Orr TE, Carter DR. An approach for time-dependent bone modeling and remodeling – theoretical development. Journal of Orthopaedic Research 1990; 8:651–661. 102. Goodship AE, Lanyon LE, Mcfie H. Functional adaptation of bone to increased stress. An experimental study. Journal of Bone and Joint Surgery - Series A 1979; 61(4):539–546. 103. Koontz JT, Charras GT, Guldberg RE. A microstructural finite element simulation of mechanically induced bone formation. Journal of Biomechanical Engineering 2001; 123:607–612. 104. Press WH, Flannery BP, Teukolsky SA et al. Numerical Recipes: The art of Scientific Computing. Cambridge University Press: New York, NY, 1986. 105. Pandy MG. Computer modeling and simulation of human movement. Annual Review of Biomedical Engineering 2001; 3:245–273. 106. Erdemir A, Mclean S, Herzog W et al. Model-based estimation of muscle forces exerted during movements. Clinical Biomechanics 2007; 22(2):131–154. 107. Jaworski ZFG, Liskova-Kiar M, Uhthoff HK. Effect of long term immobilization on the pattern of bone loss in older dogs. Journal of Bone and Joint Surgery - Series B 1980; 62B(1):104–110. 108. Covertino VA. Considerations for an exercise prescription. NASA, 1989. Report no.: 19910001272. 109. Vico L, Collet P, Guignandon A et al. Effects of long-term microgravity exposure on cancellous and cortical weight-bearing bones of cosmonauts. Lancet 2000; 355(9215):1607–1611. 110. Ptc. Pro/engineer. In (Wildfile 4.0 edn). PTC, Inc.: Needham, MA, USA. 111. White SC, Yack J, Winter DA. A three-dimensional musculoskeletal model for gait analysis, anatomical variability estimates. Journal of Biomechanics 1989; 22(8/9):885–893. 112. Reynolds H, Snow C, Young J. Spatial geometry of the human pelvis. Protection and Survival Laboratory, Civil Aeromedical Institute, Mike Monroney Aeronautical Center, Federal Aviation Administration: Oklahoma City, OK, 1981. Report no.: Memorandum Report AAC-119-81-5. 113. Dong XN, Guo XE. The dependence of transversely isotropic elasticity of human femoral cortical bone on porosity. Journal of Biomechanics 2004; 37(8):1281–1287. 114. Hoffmeister BK, Smith SR, Handley SM et al. Anisotropy of Young’s modulus of human tibial cortical bone. Medical and Biological Engineering and Computing 2000; 38:333–338. 115. Ebacher V, Tang C, Mc Kay H et al. Strain redistribution and cracking behavior of human bone during bending. Bone 2007; 40:1265–1275. 116. Cowin SC (ed.). Bone Mechanics Handbook. Informa Healthcare: New York, NY, 2008. 117. Varghese BA, Miller ME, Hangartner TN. Estimation of bone strength from pediatric radiographs of the forearm. Journal of Musculoskeletal Neuronal Interactions 2008; 8(4):379–390. 118. Stephenson P, Seedhom BB. Cross-sectional geometry of the human femur in the mid-third region. Proceedings of the Institution of Mechanical Engineers. Part H, Journal of Engineering in Medicine 1999; 213(2): 159–166. 119. Noble PC, Box GG, Kamaric E et al. The effect of aging on the shape of the proximal femur. Clinical Orthopaedics and Related Research 1995; 316:31–44.




e02699 (29 of 29)

120. Croker SL, Clement JG, Donlon D. A comparison of cortical bone thickness in the femoral midshaft of humans and two non-human mammals. HOMO - Journal of Comparative Human Biology 2009; 60(6):551–565. 121. Uzel A-P, Deloumeaux J, Rouvillain J-L et al. Comparative study of femoral diaphyseal morphometry in two male populations, in France and a French West Indies island: an example of clinical relevance of comparative anatomy for orthopedic practice. Surgical and Radiologic Anatomy 2011; 33(3):235–240. 122. Lee SC, Coan BS, Bouxsein ML. Tibial ultrasound velocity measured in situ predicts the material properties of tibial cortical bone. Bone 1997; 21(1):119–125. 123. Knets I. Peculiarities of the structure and mechanical properties of biological tissues. Meccanica 2002; 37(4–5): 375–384.



Copyright of International Journal for Numerical Methods in Biomedical Engineering is the property of John Wiley & Sons, Inc. and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.

A coupled mechano-biochemical model for bone adaptation.

Implementation of design of experiments for optimization of forced degradation conditions and development of a stability-indicating method for furosemide.

Time course of muscle adaptation after high force eccentric exercise.

Skeletal development and bone functional adaptation.

Non-orthogonal spin-adaptation of coupled cluster methods: A new implementation of methods including quadruple excitations.

Effect of Phase Response Curve Shape and Synaptic Driving Force on Synchronization of Coupled Neuronal Oscillators.

Improving the Q2MM method for transition state force field modeling.

Optimization of scaffold design for bone tissue engineering: A computational and experimental study.

Form follows function: a computational simulation exercise on bone shape forming and conservation.

Concurrent subspace width optimization method for RBF neural network modeling.

A feasibility investigation for modeling and optimization of temperature in bone drilling using fuzzy logic and Taguchi optimization methodology.

Implementation of intrinsic lumped parameter modeling into computational fluid dynamics studies of cardiopulmonary bypass.

Computational modeling of seizure dynamics using coupled neuronal networks: factors shaping epileptiform activity.

Tetanic force development of adductor pollicis muscle in anesthetized man.

Computational modeling of peripheral pain: a commentary.

Curve-based shape modeling a tutorial.

Taking dietary habits into account: A computational method for modeling food choices that goes beyond price.

Computational Advances for the Development of Allosteric Modulators and Bitopic Ligands in G Protein-Coupled Receptors.

A simple network thermodynamic method for modeling series-parallel coupled flows. I. The linear case.

Computational modeling of sphingolipid metabolism.

Computational modeling of membrane proteins.

Development of a 10-year-old full body geometric dataset for computational modeling.

bone interface using simulated bone.

Development and evaluation of a method of calibrating medical displays based on fixed adaptation.