- Email: [email protected]
Bone remodelling adjacent to intramedullary stems: an optimal structures approach
Biomatedals 17 (1996) 223-232 0 1996 Elsevier Science Limited Printed in Great Britain. All rights reserved 014%9612/9FJ/$15.00
ELSEVIER
Bone remodelling adjacent to intramedullary stems: an optimal structures approach Timothy P. Harrigan”, James J. Hamilton+, Jefhey D. Reuben*, Aldo Toni’ and Marco Vicecont? ‘Department of Orthopedic Surgery, University of Texas Health Center at Houston, 6431 Fannin, Suite 6.154, Houston, TX 77030, USA; +Department of Orthopedic Surgery, Medical School, University of Missouri-Kansas 2301 Holmes Street, Kansas City, MO 64108, USA; t lnstituto Rizzoli, Via Barbiano 1110, Bologna 40736, ltaly
The internal structure bone zation
theories,
analytical
finite
element
around
cal results.
we have
This leads
results
different,
within
formula bone
due to Huiskes
practical
application
are not clearly
overall
remodelling
simulation.
The analytical
The predictions
interpretation
for internal remodelling
of the optimizing
between
rule.
In this study,
can predict formula
of the formula
predicts
which
are clinically
around based
remodelling
theory
we assess
whether
when
intramedullary
bone will
out in the numeri-
This study implants
The study
the authors
of a large-
in the remodelling
relevant.
have
of
optimi-
are borne
on optimization.
class
structural
the behaviour
of one of the two parameters
remodelling algorithm
to the bony
a connection
optimization
characteristics
related
For a restricted
on structural
implant.
some
often
they are implemented.
found
a continuum-level
based
to a physical
stable
theories in which
previously
remodelling
also show
numerically
remodelling
the simulations
an intramedullary
The results
earlier direct
from
and the parameters
remain used.
in bone
results
remodelling
a simplified scale
parameters
which
City,
rule
extends
by using
a
also shows
a
developed
previously. Keywords:
Finite
element
Received 23 November
models,
bone
analytical
predictions
1994; accepted 6 June 1995
structure is unable to assess the mechanisms of remodelling per se. However, stability studies can establish sensible bounds on the overall continuumlevel behaviour of a remodelling scheme proposed to act within the microstructure of bone. That is, a celllevel remodelling scheme is usually generalized to a continuum-level remodelling scheme so that a larger scale structural study can be accomplished. The stability results from continuum-level remodelling studies can provide sensible bounds on the cell-level remodelling hypothesis or on the generalization to a continuum-level remodelling scheme. For example, if strain energy density/density were the actual stimulus in vivo, and the bone obeyed the remodelling rule, then the bone density would be unstable in vivo and it would diverge from its resting state given a very small stimulus. Since this type of behaviour does not happen in vivo, then that remodelling rule does not even qualitatively characterize bone remodelling in vivo. If the remodelling simulation is stable, then at least one qualitative aspect of the behaviour is correct. Also, a remodelling scheme which is stable and well-verified will be a useful tool for the design of long-term orthopaedic implants. One problem with continuum-level remodelling
Human and animal bones have become a paradigm for demonstrating the ability of a biological structure to match form to function. The modelling and remodelling processes which create and maintain bone have been explored in many experimental and theoretical studies, but substantial questions regarding these processes remain. Bone remodelling simulations have been attempted by many investigators to simulate the observation that overloaded bone generally becomes thicker and denser, and underloaded bone generally becomes more porous and thinner. Bone remodelling has been approached on two primary levels. In one approach, the cellular processes that remodel bone have been studied recently in some detail by biologists and engineers’-3. A different approach has also been taken in which the interaction between the changes in bone tissue at a continuumaveraged point due to mechanical stress and the characteristics of a whole bone as a structure is studied4-a. This paper is concentrated on the latter approach. Study of the interaction between bone density changes and the stability of the resulting remodelled Correspondence
adaptation,
to Dr T.P. Harrigan. 223
Biomaterials
1996.
Vol.
17 No. 2
Optimal
224
algorithms is the difficulty of assessing how significant the internal parameters in the remodelling rule are. This is made more difficult by the large numbers of parameters in some models, and by remodelling rules which are very complicated. Although the actual remodelling process could indeed be as complicated as the rate equations postulated, the interpretation of several parameters in such a rule is very difficult. Thus, we have used a simple remodelling rule thus far, and we have concentrated on a detailed assessment of the internal parameters in the model. We have recently developed a specific mathematical connection between bone remodelling as a continuumlevel material adjustment to mechanical strain energy density and bone remodelling as an overall structural optimization processlO. This has led to simple conditions for a unique and stable solution for longterm bone density, for a class of hypothetical bone remodelling equations”. However, the physical significance of this connection is not obvious, and the usefulness of such a connection is not clear. In this paper we develop a global interpretation of the parameters in the local remodelling rule. This is very important for comparison with experimental results, because if we understand how the local parameters relate to global behaviour, then we can assess whether adjustments in local parameter values could account for differences between experimental results and model predictions. Understanding how local parameters influence results will also allow an informed assessment of how well the numerical model is behaving. Remodelling around intramedullary stems has been explored by Huiskes et al.‘” but the ‘internal remodelling’ algorithm displayed numerical difficulties. This study revisited the internal remodelling simulations of Huiskes et al., using a newly developed remodelling algorithm. This study also was intended to show a practical application of the remodelling algorithm which was developed in mathematical detail in our prior publications. The remodelling rule used in this study is similar to that of Huiskes et aI.” and Weinans et al:‘” but several important distinctions exist. The difference between this model and that of Huiskes is simply that this model arrives at unique, stable solutions which optimize a known cost function for the entire structure. The approach is new, because the exponent m was left as a parameter and chosen to make the global behaviour stable, whereas Huiskes et al. chose m = 1 so that the local parameters were well defined. The results are new, because prior to this work, remodelling approaches relied on global or local optimization, but there was no direct connection between the two approaches. The results are also new because this approach allows the solution of very refined finite element meshes in clinically relevant geometries, whereas if the model of Huiskes is refined, a continuous solution for the density distribution exists, but it is unstable and cannot be found in practice. The approach we have taken in our remodelling rule has been to studv the global behaviour of the remodelling simulaiion, and to derive limits on the Biomaterials
1996, Vol. 17 No. 2
remodelling
around
intramedullary
stems:
J.P. Harrigan
et al.
stimulus within the microstructure so that the solutions we arrive at are continuous and stable. There is no justification based on a remodelling signal within the microstructure. That is, we are letting the global behaviour of the model establish limits on the acceptable stimulus within the microstructure, and we are not being specific about what the stimulus means microscopically. The purpose of this study was to test how well a simple characterization of bone remodelling as an optimizing process can predict the behaviour of a large-scale finite element-based bone remodelling simulation. A secondary purpose was to gain a better understanding of the parameters in this and other finite element bone remodelling simulations.
METHODS Our bone remodelling theory is based on two descriptive relationships: a simple characterization of the bone remodelling stimulus and a structureproperty relationship for bone. Both are idealizations of known data, but they are reasonable approximations, and they are mathematically simple. In the future, more complicated remodelling rules and prosthesis geometries can be explored in detail based on the development here. This optimization study was used to derive a prediction of whether the cortical shell of bone will remain in the long term around a loaded intramedullary implant. Many studies have shown a drastic loss of bone around uncemented total hip replacements and porous-coated intramedullary rods, and fibrereinforced composite materials are proposed as a way to decrease the stiffness of a prosthesis enough to maintain bone. lising an optimization function developed in previous publications, we can make a simple prediction of the needed prosthesis size and material stiffness to maintain bone around the implant for axial loads, and transverse bending moments. The bone remodelling algorithm used in this study is described in our related publications1”.1’ and we will only summarize the basics of our methodology here. The method rests on a simple relationship between bone density and elastic modulus, and another simple relationship between strain energy density, volumetric density of hard tissue, and the rate of change of density. We assume isotropy, a material Poisson’s ratio that is constant, and a power law relationship between volumetric hard tissue density, +, and elastic modulus’“: E = E,, (I,”
(1)
with n the material exponent, usually taken as 2 (Rice et al.‘“) or 3 (Carter and Hayes’“). The bone remodelling rule we use is detailed in the literature, and it uses a remodelling stimulus which is strain energy density divided by 9”‘. Here m is an exponent which is chosen to make the simulations stable. This is a fundamental difference in approach from that of Huiskes et nl.‘“, Weinans et al.“’ and Beaupr6 et ~21.~.The difference between this stimulus and a set point is taken as proportional to the rate of
Optimal
remodelling
change of remodelling
around
intramedullary
4. This relationship rate equation which is
stems:
leads
to
et al.
a
a geometrical
bone
(21 with Q0 the stimulus set point. In a previous papeP, we showed that changing the state variable used in the remodelling equation allows that equation to be part of an optimizing process. By using density to the power m we can define a variable Y as li’ ZZ4”
(3)
and by changing
the rate equation
to be
we can satisfy the requirements for an optimizing function. This will have the same steady-state solutions as the previous equation. We have shown that this remodelling rule gives the same steady-state density distributions as Equation 2, and we have also shown that finding a distribution of y which makes the rate of change in Equation 4 zero is equivalent to finding stationary points in the function H =
(
:
>
225
T.P. Harrigan
HI + H,
with
optimization, rather than a density-based optimization, but the density-based optimization is actually general enough to include this type of optimization as a subset. In terms of the remodelling rule we have developed above, this analytical procedure corresponds to choosing a restricted class of density distributions in the optimizing process. That is, we use density distributions which are fully dense (i.e. >J= 1) over a bone area Ah and zero outside that area. This is a subset which is included in the set of all possible density distributions. Thus, minimizing H with respect to bone area fits within the framework developed for density-based remodelling if we consider only density distributions which are a restricted subset of all possible bone density distributions. Substituting for the strain in terms of the total force on the section yields H=
(m/n)F2L
(91
+\16AbL
2@,4, +&Ad
Since this is the strain energy only in a section of the structure, we cannot expect that minimizing H by manipulating Ab will give an exact answer, but since most femoral prostheses contain long stems, we expect this representative section to allow us to make an approximate generalization. Minimizing H with respect to Ah yields the following equation for Ah: Ab = (s)‘”
(&)
-ZAP
(10)
and for Ab to be greater than zero, we must have
(6) E,A,
e=1
H, is the total strain energy in the structure (bone and prosthesis included). In this study, we consider axial bending loads on a straight cylindrical prosthesis within a straight bone. We treat the loading cases separately for clarity, but the combined effect of these loads is easily incorporated.
(11)
This requirement can allow us to give a limit on the axial stiffness of an implant to maintain cortical bone around the shaft. In the simulations below, we used a fixed geometry, and we varied implant elastic modulus and the value of Q0 to assess cortical remodelling. The limit on the value of ‘16 needed for maintenance of a cortex is
(12) These requirements are reflected in the results of the simulations in which the cortex is allowed to remodel.
Axial loading For axial loads, if we consider a section of bone containing a prosthesis and we assume that the strain within the section is uniform, we can write Has H = ((m/n) (EPA, +Et,AtJ
mEb II2 __ 2n$ ( >
k2/2) + %AdL
(8)
where Ep and Eb are the elastic moduli of prosthesis and bone, respectively, A, and Ah are the crosssectional areas of prosthesis and bone, respectively, and L is the length of the section. We can now minimize the function H with respect to bone area, and arrive at a prediction of the stem stiffness required to maintain bone adjacent to the stem. On first glance, this analytical procedure seems to be
Transverse moments We developed two simplified models for bone adaptation to an intramedullary implant under a transverse bending moment in this study. Both models were aimed at assessing conditions for bone to remain around the implant. The models assumed that the bending radius of curvature for the implant and the surrounding bone were the same, and that bone exists in a concentric ring around an implant, with an interior radius ri and an exterior radius r, concentric with the implant section which is characterized by the flexural rigidity EpZ,,. The inner bone radius is not Biomaterials
1996.
Vol. 17 No. 2
226
Optimal
constrained to be the same as the outer prosthesis radius, since in many clinical cases the bone does not touch the implant in a given cross-section. Given this geometry and an equal radius of curvature for bone and implant, we can write the strain energy of the implant and prosthesis as TSE =
2M”L rcE,,(r; ~
$1 + 4EpIp
(13)
with M the transverse bending moment and EI, the bone elastic modulus. Both approximate approaches used this expression for total strain energy in the bone and implant. In the first bending optimization model, the density of the bone was taken as 1.0 and thus the restricted set of density distributions used in the analytical model consisted of circularly symmetric regions of fully dense bone, with the outer radius taken as fixed, and the inner radius variable. Thus, the index which selects the density distribution in the optimization model is ri, the inner bone radius. Minimizing the indicator functional over this restricted set of bone density distributions leads to the following equation: (14) and to assess when cortical bone will be adjacent to the implant, we need to assess the limit as ri approaches r,,. In that case, the bone flexural stiffness becomes small, and we can replace the first term on the right with f: and solve for c. For cortical bone to remain around the implant, we require r; < r, and this results in
(15) and we can take the limiting value E = 0. This results in a limiting value of flexural stiffness for an intramedullary implant for maintenance of bone, using the singleparameter restricted set of density distributions. In the second analytical optimization model we used the same geometric analysis but we took volumetric density as a fixed value less than unity. This resulted in the limit
This limit has some physical implications, since the exponent on y is negative for m > n as 7 goes to zero, the value of $‘I ~17):2~7~ will become infinite, but since the magnitude of the exponent is small, this limit can be accounted for in a practical sense. For n = 3.0 and m = 3.1, the exponent is equal to -0.016129 and if the actual density 4 = lo-” (below the lower bound for the numerical simulations) then yen “li!Zr” will equal only 1.58. This factor was used to predict the behaviour of the simulations.
Finite element
simulations
To implement the remodelling rule (Equation 4) we used a modification of the Euler backward timestepping method described in a prior publication”. Biomaterials
1996, Vol. 17 No. 2
remodelliw
around
intramedullary
stems:
T.P. Harrigan et al.
This time-stepping simulation methodology is based entirely on the remodelling rate equation, and does not depend on the optimization approach described in prior publications. The input data for these simulations were adapted from an input file for ADINA that was generated by ADINA-INIH.l”, the input processor for ADINA. A Euler backward time-stepping method, similar to a previously described method17, was implemented using Equation 4. Comparison of steady-state results between this method and the one used previously”” shows the results to be the same, within the numerical tolerances used for convergence. For axial loading simulations a two-dimensional axisymmetric model was used. The prosthesis diameter was 15 mm, with 5 mm of cancellous bone modelled between the outer surface of the implant and the cortex, which was 5 mm thick. Axial loads were applied to the top of the prosthesis, and the bottom of the bone was fixed. All interfaces were assumed to be perfectly bonded. We assumed that the distal tip of the intramedullary stem was spherical, and the distance from the distal tip of the prosthesis to the upper level of the femur was 107.5 mm. The prosthesis extended 10 mm beyond the end of the cortex. We chose to use a fine finite element discretization so we could assess the distribution of bone as predicted by the model in detail. Eight-noded elements, which assume that displacements vary quadratically with position, were used. The finite element mesh used 10 elements through the thickness of the cancellous bone and there were 100 rows of elements along the length of the prosthesis. There were 125 elements representing cancellous bone in the region of the prosthesis tip. The cortex was 5 elements thick and 120 elements long. To simulate the prosthesis, we used a T-element wide by l?O-element long shaft portion, and a 45-element tip. A total of 8287 nodes were used in the model. Convergence checks on the elasticity results of the remodelling simulations were not done, although similar studies in our prior publication’7 indicate that the mesh discretization is adequate for approximating the density. That is, we have not assessed the relationship between mesh size and the accuracy of the results. We have not completed a rigorous assessment of the various possible convergence criteria for the elasticity problem; this is under way. The results plotted in Figures I and z are not interpolated, however: they are made up of even blocks of colour which indicate element density. The rather smooth results seem to indicate that the mesh used in this study is sufficiently fine to answer the research question posed in this study. Also, since Huiskes et al.” obtained useful results for remodelling around an intramedullary stem with less than one-tenth the number of elements we LW: here, we feel our mesh was sufficient to characterize the elasticity problem. Ten axial load simulations were run in this study, with both the cortex and the cancellous bone allowed to remodel. Within the simulations, the prosthesis elastic moduli were chosen using a scaling procedure. We took the elastic modulus of fully dense bone as 17 GPa, and the prosthesis moduli to be 17, 34, 85, 136 and 204 GPa. We scaled the total load applied to
Optimal
remodelling
around
intramedullary
stems:
T.P. Harrigan
Figure 2 Remodelling results with $0 below the limits in Equation 72. In these cases, bone is maintained. The remodelling cases from Table 7 are (from left to right) A6, A7, A8, A9, and AlO. The central elastic moduli in the prostheses are, from left to right, 17, 34, 85, 136, and 206 GPa. As in figure 7, the shading indicates bone density, and each element is plotted as an even shade. Bone persists in all cases, with the initial cortical demarcation faintly apparent in cases A6, A9, and AlO.
Figure 1 Remodelling results with s above the limit in Equation 72. These results indicate that for @,, above the limit in Equation 6, wasting of the bone adjacent to the implant is to be expected. The remodelling cases from Table 7 are (from left to right) Al, A2, A3, A4, and A5 The elastic moduli for the prosthesis material were (from left to right) 17, 34, 85, 136, and 204 GPa, respectively. The shading shows bone density from 0.0 (darkest at right) to 1.0 (lightest at right). Notice that no interpolation is used; a constant density is displayed for each remodelling finite element. Resorption is clear in cases A2, A3, and A4, and partial in cases Al and A5
The finite element model was generated using ADINA-IN (ADINA Engineering, Watertown, MA, USA) and the finite element data were re-formatted for our simulation program using a translation program and a text editor. The axial and transverse load cases we used in this study were chosen to test the limits on prosthesis elastic modulus points and remodelling set simultaneously. For each prosthesis elastic modulus, a limiting value of q0 was calculated for axial loading from Equation 12. Six values of q0 were established so that we could test the limiting value of prosthesis stiffness, as predicted from Equation 12. Table 1 shows how the ten finite element runs in this study were used to test Equation 11, and Table 2 shows how eight of these ten cases were used to test Equation 12. Similarly, Equations 15 and 16 were used to generate limits on !QOand prosthesis stiffness which were tested using the seven large-scale three-dimensional cases shown in Table 3. Sets of three runs were arranged to test the two limits posed by Equations 15 and 16, with one case below the lower limit (on Ep or @J, one case between the two limits, and one case above both limits. Thus, by comparing cases B3, B2 and B6, or cases B5, B4 and B7, the results of progressive increases in Q0 for a given Ep were assessed, and by
1000 N, and we scaled the elastic moduli so that the ‘set point’ in our simulations was 0.5. Thus, the load which is relevant to the data in Table 1 is 1000 N, applied vertically as a point load to the top of the intramedullary stem. The cortical bone was fully dense (y = 1.0) and cancellous bone was simulated by taking y = 0.5. Seven three-dimensional finite element simulations were run in this study to assess the remodelling response of bone adjacent to an intramedullary stem which was loaded by a transverse bending moment. The finite element model consisted of a cylindrical stem with a spherical end-cap which was surrounded by bone with a cylindrical outer surface. The cylindrical portion of the stem was 100 mm long and the stem was 15 mm in diameter. The bone outer diameter was 35 mm. In the model, bone extended past the stem tip a distance of 10 mm and extended proximally to the top of the stem. The finite element model incorporated the entire stem and bone (i.e. a full 360” model was used) so that we could check for numerical problems by comparing left and right halves of the model.
Table 1 Elastic
17 34 85 136 204
Remodelling modulus
(GPa)
cases used to test predicted
limits on E,
Upper limit on Q0 from Equation (J mm3) 973.1 243.3 38.9 15.2 6.8
227
et al.
73
Simulation (case no.) 1055.8 422.3 84.5 33.8 8.4
(Al) (A2) (A3) (A4) (A5)
value above
Simulation (case no.) 422.3 84.5 33.8 8.4 4.2
value below
(A6) (A7) (A8) (A9) (AlO)
Biomaterials
1996, Vol. 17 No. 2
228
Optimal
Table 2
Remodelling
cases used to test predicted
Value of Q0 (J mm3)
Table 3
Q0 q0 q0 Q0
= = = =
Remodelling
500 J me3 1000 J mm3 2500 J mm3 5000 J me3
Equation 72
(GPa)
Simulation (case no)
25.8 57.7 91.2 182.5
17 34 85 136
cases used to test limits in pure bending
E, = 85 GPa
E, = 115 GPa
Bl
83 B2 B6
E, = 200 GPa B5 B4 87
comparing cases Bl, B2 and B7, the effect of increasing Ep for a given value of bone parameters was assessed. The elastic modulus of cancellous bone was taken as proportional to density cubed (n = 3), and the exponent on volumetric density that was used to divide the strain energy density (m in Equation 4) was taken as 3.1~. The density was limited to be no greater than 1.0 or less than 0.01 in the simulations. The norm of an increment in density within a given time step was computed by summing up the square of each density increment, dividing by the number of remodelling elements, and taking the square root of that quantity. The same was done to compute the norm of the remodelling rate and the norm of the overall density distribution. These norms included only the elements in which the density was greater than 0.01 and less than 1.0. Convergence to a solution at each time step was assumed when the norm of the density increment was 0.01 times the norm of the density. Convergence to a steady-state solution was assumed when the size of the remodelling rate (as scaled for simulation) was 0.005 times the size of the density. The finite element runs used in this study were coded in C and run on the Cray C-90 at the Pittsburgh Supercomputer Center. An adjustable time step value which increased in four intervals from 0.1 to 3.0 (in normalized time) was used to save on execution time. As a part of the numerical procedure, the stability of the simulations was checked at the end of the simulations by forming the matrix T + P, factorizing it numerically into LDLT form, and assessing whether negative diagonal entries resulted”. In all cases, the system of density evolution equations was found to be stable, as predicted analytically.
RESULTS The finite element remodelling results agreed in general with the predictions of the simplified analytical optimization results, and show that X& has the effect anticipated by the analytical predictions. The five axial load cases in which q’,, was chosen Biomaterials 1996,
Vol. 17 No. 2
around
intramedullary
stems:
T.P.
Harriman et al.
limits on q,,
Upper limit on E from
422.3 84.5 33.8 8.4
remodelling
(8) (7) (8) (9)
value below
Simulation (case no) 34 85 136 204
value above
(2) (3) (4) (5)
above the limit for cortical maintenance showed clear cortical resorption (Figure z), and the five axial load cases where Q0 was below the limit showed cortical maintenance (Figure 2). The implant stiffness increases from left to right in Figures 1 and 2. These plots show half of a section through the mid-plane of an intramedullary rod, and the colours are used to depict bone density. The colour scales on the right of Figures 1 and 2 show grey scales corresponding to densities that vary linearly from zero (at the bottom] to 1.0 (at the top). We have shaded the implant evenly. Cases Al and A5 (the farthest left and right in Figure 1) were only slightly over the threshold for resorption, based on Equations 11 and 12, so the persistence of some bone in the simulations is not surprising. Thus, the results of the large-scale simulations corresponded to the analytical results, and the set point Q(, has a substantial effect on whether bone will persist around an intramedullary stem. The results here also agree with the predictions of Equations 11 and 12. For a given value of $:, the cases noted in the latter part of Table 1 show cases which are above and below the limiting prosthesis stiffness in Equation 12. As in the cases which tested Equation 6, the remodelling results (Figures z and 2) support the derived limit on implant stiffness. The results of the bending simulations show that the predictions of Equations lli and 16 are borne out. As would be expected, the assumption that the density distributions are cylindrically symmetric is not accurate. Nevertheless, the density distributions followed the predictions of the limits over much of the prosthesis. Figures 3, 4 and 5 show remodelling cases used to test changes in density distributions due to changes in the remodelling set point and prosthesis elastic modulus. These figures are quarter sections of the density distribution around an intramedullary implant, and the direction of the applied bending moment is clear. In the leftmost cases, the prosthesis modulus and the remodelling set point are chosen so that both Equations 15 and 16 predict full resorption. The middle case in these figures has resorption predicted by Equation 15 but not Equation 16, and the rightmost case has bone maintenance predicted by both Equations 15 and 16. The distributions show that Equation 15 can predict when fully dense bone will occur around an implant, and Equation 26 can predict when bone will effectively disappear. Some of the features of the final density distributions have implications to follow-up studies of bone remodelling around intramedullary devices, and to finite element studies. The simulation cases in which the optimization approach appeared to be only marginally correct can lead to some mechanical
Optimal
remodelling
around
intramedullary
stems:
T.P. Harrigan
229
et al.
Figure 3 Three three-dimensional cases illustrating the changes in bone density for E, = 200 GPa and & = 2500 J mm3 (left, case B7), 1000 Jm3 (middle, case B4), and 500 Jmd3 (right, case 85). For this implant stiffness, the limits on Q0 from Equations 75 and 76, respectively, are 685.4 J mm3 and 1082.9 J mm3. Both analytical models predict resorption in case 87 and maintenance in case B5. The bone density in case 84 at the maximal radius was less than 1.0, indicating that Equation 75 predicts the occurrence of fully dense bone correctly. The correspondence between shading and bone density is as in Figures 7 and 2.
Figure 4 Three three-dimensional cases illustrating the changes in bone density for E, = 115 GPa and q0 = 5000 J mm3 (left, case B6), 2500 J m3 (middle, case B2), and 1000 Jmm3 (right, case 83). For this implant stiffness, the limits on Q0 from Equations 75 and 76, respectively, are 2174 J mm3 and 3435 J mm3. Both analytical models predict resorption in case B6 and maintenance in case B3. The bone density in case 82 at the maximal radius, along the shaft of the implant, was less than 1.0, indicating that Equation 75 predicts the occurrence of fully dense bone correctly. The correspondence between shading and bone density is as in Figures 7 and 2.
interpretations. Cases Al and A5 show density variations at steady state which extend from the proximal end of the intramedullary rod a substantial distance along the stem. These cases indicate that Saint Venant’s principle should be applied very carefully in bone remodelling simulations. These simulations indicate that in these cases, there is no ‘central region’ in which the load sharing between the bone and implant is uniform. In three of the cases in which bone persists at steady state (cases A8, A9 and AlO), the initial corticalcancellous interface can be seen, faintly. Previous simulations in which we used a more lax tolerance for convergence of the model to a steady state showed more pronounced residual cortical-cancellous demarcations, so we anticipate that still tighter convergence tolerances would further eliminate these demarcations. Also, in cases A6, A7 and A8, as well as in all cases where a pure bending load is applied, there are regions of low bone density adjacent to the prosthesis. These areas could appear as radiolucencies on a follow-up radiograph. In bending, the low strain in regions nearer to the neutral axis of the boneprosthesis composite beam is probably the cause for this effect. The low-density areas could be mistaken for regions where either bone ingrowth has not occurred, or where interface failure is occurring. That is, the low-density areas around a prosthesis seen on
X-ray could be a remodelling genuine interface failures.
effect
or they could
be
DISCUSSION The remodelling simulations carried out here were aimed at a restricted set of design parameters for intramedullary stems, and used a restricted set of physical effects. In this study, we concentrated on the limited goals set out so we could maintain control of the parameters in the simulation, and so we could assess the remodelling rule used in isolation. This allowed us to build upon detailed knowledge concerning the stability and the optimizing characteristics of the remodelling simulation used”. The results of this study show that a very coarse set of approximations can predict the behaviour of a largescale remodelling simulation around a loaded intramedullary implant. The use of very coarse approximations such as these is similar to the use of the Rayleigh Quotient in vibration analysis”‘. We have verified these predictions under torsional loads as well (unpublished data) and we are currently assessing the accuracy of simple optimization predictions in multiple loading situations. The accuracy of the approximate relationships derived from these coarse models seems good, and the relationships should give Biomaterials
1996, Vol. 17 No. 2
230
Optimal
Figure 5 Three three-dimensional cases illustrating the changes in bone density for q0 = 2500 J mm3 and with the implant elastic modulus taken as 200 GPa (left, case B7), 115 GPa (middle, case B2) and 85 GPa (right, case Bl). Both analytical models predict resorption in case 67 and maintenance in case Bl. The bone density in case B2 at the maximal Fadius was less than 1.0, indicating that Equation 75 predicts the occurrence of fully dense bone correctly. The correspondence between shading and bone density is as in Figures 7 and 2.
simple and useful data for the design of implants and as a guide for more accurate bone remodelling simulations. In its current form, our remodelling stimulus has only two parameters, and we can readily test whether these two parameters can characterize natural bone density distributions. Here we have developed some mechanical intuition for q’,, as used in Equation 4. This parameter is the set point of the remodelling stimulus, but it can also be interpreted as the weighting factor on bone mass (or a measure of bone mass) in a global optimization approach to bone remodelling. As is apparent from the remodelling rate equation, $, represents the amount of stimulation needed to maintain bone, and this study provides quantitative criteria for the intuitive idea that if the prosthesis is too stiff (or QO is too high) the bone will waste away around an intramedullary implant. We expect that in a individuals will have different biological setting, combinations of the variables in Equations 11, 12, 15 or ~6, and that there will be variable results around intramedullary stems, but there should be a value of 4, which most individuals do not exceed. This will give an upper limit on the prosthesis stiffness for bone maintenance. Although convergence checks were not made in this study (since they are not well developed for finite element problems such as these), we are confident that the models used are refined enough to assess whether the optimization approach is worthwhile. Previous Biomaterials
19%. Vol. 17 No. 2
remodelling
around
intramedullary
stems:
T.P.
Harriaan et al.
studies have shown that progressive refinement in these models results in closer and closer approach to a continuous solution, so that further refinement of the finite element meshes in this study will not change the salient results. The internal corner in the axially symmetric finite element mesh below the implant generates an artefact which will not be generally seen. However, since this study was concerned with a prediction of bone maintenance (or loss) adjacent to the cortex along the shaft of the stem, we feel that the presence of this artefact is not a serious problem in this study. Note that this artefact is due to the shape of the remodelling region, rather than the mesh discretization. The parameters we know of which are not in this model are as important for interpreting our results as those which are in the model. We assume that there are no non-mechanical stimuli which depend on position. That is, we assume that the non-mechanical effects are small, or that they are applied evenly throughout the modelled area. We do not assume that there is a non-responsive range around the @04353H.We have assessed the response to two isolated loading situations and we have not included the effect of strain rate in the model. Based on these results, we can assess the relative influence of individually added remodelling features on the overall response to implants in the future. The remodelling prediction we arrive at is similar to one derived by Huiskes et a1.12 for bone around an implant. Huiskes et aI.” used a remodelling theory based on strain energy density which modified outer cortical geometry in ‘external remodelling’ and which elastic modulus in ‘internal modified tissue remodelling’. This study, by comparison, is concerned only with internal remodelling. Here we confirm and extend the results of Huiskes et a1.l’ for internal remodelling. Huiskes et aI.‘” noted that “A homeostatic end configuration was reached only in the case of the fixed implant, the most flexible stem, and assuming the maximal threshold level s = 30(/o, in which only the proximal l/3 of the bone was fully demineralized”. This indicates that the method used by Huiskes et al. results in a very difficult numerical convergence problem. Since Huiskes et al. only reported one case in which the internal remodelling converged, the limit derived in that study for bone maintenance in internal bone remodelling could not be fully tested. That is, one case in which bone resorption was predicted and bone maintenance was another case in which predicted were not reported. Our method for internal remodelling results in a numerically stable formulation which allows a very fine finite element discretization without a lazy zone. The method is directly related to an optimizaand this study shows that the tion procedure, internal remodelling parameters can be interpreted physically in a practically useful model. Thus, we can confirm and extend the results of Huiskes et remodelling with a remodelling al.” for internal which has several theoretical and procedure practical advantages. The utility of a numerically robust finite element
Optimal
remodelling
around
intramedullary
stems:
T.P.
Harrigan et a/.
simulation methodology which allows continuous solutions of very refined finite element meshes is important, and the distinction between this approach and that of Weinans et al.‘” is noteworthy. As detailed in prior publications, if we choose m large enough, then a continuous, stable, unique solution is found to the posed mathematical problem for bone remodelling. As shown by Weinans et ~l.l”, taking m = 1 leads to a solution behaviour in which a continuous solution is approached, but then a discontinuous solution results. In fact, the continuous solutions approached by Weinans et al. are valid solutions for the remodelling response, but these solutions are unstable. If a numerical method was developed to arrive at those solutions, then the discussion of a discontinuous solution would be moot. The resulting solutions which diverge from the continuous solution are not unique, and the checkerboard patterns which result are entirely dependent on the element pattern in the finite element mesh. If the approach of Mullander et aI.” is taken, the patterns of divergence from the continuous solution are interesting, but the continuous solution still exists. Thus, the trabeculating patterns are interesting as possible trabecular patterns, but the background density distributions (i.e. the density as averaged over several trabeculae) are probably close to the continuous (but unstable) solution. As detailed in a prior publication16 the optimization function is derived from the remodelling stimulus, so that the predictions should be well represented. Note, however, that the finite element remodelling simulations do not depend on the optimization approach, they are based solely on Equation 4. Thus, the general agreement between the optimization predictions and the simulations seemed to confirm the usefulness of the optimization approach. This optimization approach is related to, but distinct from, other optimization approaches in which (a) the arrangement of a fixed mass of bone tissue is optimized to minimize strain energy in the structure or (b) weight is minimized for a given structural stiffness”. In the optimization scheme used in this study, a weighted sum of overall potential energy and a measure of the mass in the structure is used as the optimization function. This is similar to the method used by Subbarayan and Bartel’“. The distinction is primarily in the parameters used to define the problem. In approach (a), mass is fixed, and a Lagrange multiplier is used which corresponds roughly to ‘16 in the related continuum-level remodelling rule. In approach (b), the compliance is fixed, and a Lagrange multiplier is also used which is related to Q(, as well. However, in both (a) and (b) the Lagrange multipliers are unknowns in the optimization process. That is, the value of the Lagrange multiplier in either (a) or (b) is a result of the optimization, and is not specified at the outset. This corresponds to leaving QO as an unknown in the remodelling rate equation, and finding a value which satisfies the optimization criterion. Although, for a fortuitous choice of parameters, similar structures can be found in these approaches, the related bone remodelling algorithms are not the same.
231
ACKNOWLEDGEMENTS This research was supported in part by Pittsburgh Supercomputing Center grant Number 1 P41 RR06009 from the NIH National Center for Research Resources, and by Grant Ricerca Finalizzata of the Italian Ministry of Health, del. CIPE 13/l@% “Studio dei fattori the governano il fenomeno di adattamento funzionale de1 tessuto osseo ai carichi meccanicci”.
REFERENCES 1
Currey
2
Harrigan TP, Hamilton JJ. Bone strain sensation via transmembrane potential changes in surface osteoblasts: loading rate and microstructural implications. 1 Biomech 1993; 26: 183-200. Martin RB, Burr DB. Structure, Function, and Adaptation of Compact Bone. New York: Raven Press, 1989. Beaupre GS, Orr TE, Carter DR. An approach for timedependent bone modeling and remodeling - theoretical development. 1 Orthop Res 1990: 8: 651-661. Carter DR, Orr TE, Fyhrie DP. Relationships between loading history and femoral cancellous bone architecture. J Biomech 1989; 22: 231-244. Harrigan TP, Hamilton JJ. An analytical and numerical study of the stability of bone remodelling theories: dependance on microstructural stimulus. J Biomech 1992; 25(5): 477-488 (Erratum 26: 365-366). Hart RT, Davy DT. Theories of bone modelling and remodelling. In: Cowin SC, ed. Bone Mechanics. Boca Raton, FL: CRC Press, 1989; 253-277. Huiskes R, Weinans H, Reitbergen BV. The relationship between stress shielding and bone resorption around total hip stems and the effects of flexible materials.
3
4
5
6
7
a
JO. The Mechanical Adaptations of Bones. Princeton, NJ: Princeton University Press, 1984.
Clin Orthop Rel Res 1992; 274: 124. 9
10
11
Orr TE, Beaupre GS, Carter DR, Schurman DJ. Computer predictions of bone remodeling around porous-coated implants. J Arthroplasty 1990; 5(3): 191-200. Harrigan TP, Hamilton JJ. Bone remodeling and structural optimization. J Biomech 1994; 27(3): 323328. Harrigan TP, Hamilton JJ. Necessary and sufficient conditions for global stability and uniqueness in finite element simulations of adaptive bone remodeling. Int J
Solids Struct 1994; 31(l): 97-107. 12
Huiskes R, Weinans H, Grootenboer HJ, Dalstra M, Fudala B, Sloof TJ. Adaptive bone remodeling theory applied to prosthetic design analysis. J Biomech 1987;
13
Weinans adaptive
14
Carter DR, Hayes WC. The compressive behavior of bone as a two-phase porous structure. J Bone Joint Surg 1977; 59A: 954-962. Rice JC, Cowin SC, Bowman JA. On the dependence of the elasticity and strength of cancellous bone on apparent density. J Biomech 1988; 21: 155-168. Harrigan TP, Hamilton JJ. Optimality conditions for finite element simulation of adaptive bone remodelling. Int J Solids Struct 1992; 29(23): 2897-2906. Harrigan TP, Hamilton JJ. Finite element simulation of adaptive bone remodelling: a stability criterion and a time stepping method. Int J Numer Meth Engng 1993: 36: 837-854.
20: 1135-1150. H, Huiskes R, Grootenboer HJ. The behavior of bone remodeling simulation models. JBiomech
1992; 25: 1425-1441.
15
16
17
Biomaterials
1996, Vol. 17 No. 2
232
18
19 20
Optimal Adina R&D Inc. ADINA-IN for ADINA User’s Manual (Report ARD 87-4). ADINA R&D Inc., Watertown, MA, USA, 1987. Bathe KJ. Finite Element Procedures in Engineering Analysis. Englewood Cliffs, NJ: Prentice-Hall, 1982. Harrigan TP, Hamilton JJ. Bone remodeling adjacent to total hip replacements: a naturally occurring material design problem. I Comput Aided Mater Des 1993; 1: 27-40.
Biomaterials
1996, Vol. 17 No. 2
remodelling
21
22 23
around
intramedullary
stems:
T.P. Harrigan
et al.
Crandall SH, Karnopp DC, Kurtz EF, PridmoreBrown DC. Dynamics of Mechanical and Electromechanical Systems. New York: McGraw-Hill, 1968. Strang G, Kohn RV. Optimal design in elasticity and plasticity. Int J Num Mefh Engng 1986; 22: 183-188. Subbarayan G, Bartel DL. Bone remodeling around a hole: a comparison of theoretical and experimental results. Trans Orthop Res Sot 1991; 16:425.
ELSEVIER
Bone remodelling adjacent to intramedullary stems: an optimal structures approach Timothy P. Harrigan”, James J. Hamilton+, Jefhey D. Reuben*, Aldo Toni’ and Marco Vicecont? ‘Department of Orthopedic Surgery, University of Texas Health Center at Houston, 6431 Fannin, Suite 6.154, Houston, TX 77030, USA; +Department of Orthopedic Surgery, Medical School, University of Missouri-Kansas 2301 Holmes Street, Kansas City, MO 64108, USA; t lnstituto Rizzoli, Via Barbiano 1110, Bologna 40736, ltaly
The internal structure bone zation
theories,
analytical
finite
element
around
cal results.
we have
This leads
results
different,
within
formula bone
due to Huiskes
practical
application
are not clearly
overall
remodelling
simulation.
The analytical
The predictions
interpretation
for internal remodelling
of the optimizing
between
rule.
In this study,
can predict formula
of the formula
predicts
which
are clinically
around based
remodelling
theory
we assess
whether
when
intramedullary
bone will
out in the numeri-
This study implants
The study
the authors
of a large-
in the remodelling
relevant.
have
of
optimi-
are borne
on optimization.
class
structural
the behaviour
of one of the two parameters
remodelling algorithm
to the bony
a connection
optimization
characteristics
related
For a restricted
on structural
implant.
some
often
they are implemented.
found
a continuum-level
based
to a physical
stable
theories in which
previously
remodelling
also show
numerically
remodelling
the simulations
an intramedullary
The results
earlier direct
from
and the parameters
remain used.
in bone
results
remodelling
a simplified scale
parameters
which
City,
rule
extends
by using
a
also shows
a
developed
previously. Keywords:
Finite
element
Received 23 November
models,
bone
analytical
predictions
1994; accepted 6 June 1995
structure is unable to assess the mechanisms of remodelling per se. However, stability studies can establish sensible bounds on the overall continuumlevel behaviour of a remodelling scheme proposed to act within the microstructure of bone. That is, a celllevel remodelling scheme is usually generalized to a continuum-level remodelling scheme so that a larger scale structural study can be accomplished. The stability results from continuum-level remodelling studies can provide sensible bounds on the cell-level remodelling hypothesis or on the generalization to a continuum-level remodelling scheme. For example, if strain energy density/density were the actual stimulus in vivo, and the bone obeyed the remodelling rule, then the bone density would be unstable in vivo and it would diverge from its resting state given a very small stimulus. Since this type of behaviour does not happen in vivo, then that remodelling rule does not even qualitatively characterize bone remodelling in vivo. If the remodelling simulation is stable, then at least one qualitative aspect of the behaviour is correct. Also, a remodelling scheme which is stable and well-verified will be a useful tool for the design of long-term orthopaedic implants. One problem with continuum-level remodelling
Human and animal bones have become a paradigm for demonstrating the ability of a biological structure to match form to function. The modelling and remodelling processes which create and maintain bone have been explored in many experimental and theoretical studies, but substantial questions regarding these processes remain. Bone remodelling simulations have been attempted by many investigators to simulate the observation that overloaded bone generally becomes thicker and denser, and underloaded bone generally becomes more porous and thinner. Bone remodelling has been approached on two primary levels. In one approach, the cellular processes that remodel bone have been studied recently in some detail by biologists and engineers’-3. A different approach has also been taken in which the interaction between the changes in bone tissue at a continuumaveraged point due to mechanical stress and the characteristics of a whole bone as a structure is studied4-a. This paper is concentrated on the latter approach. Study of the interaction between bone density changes and the stability of the resulting remodelled Correspondence
adaptation,
to Dr T.P. Harrigan. 223
Biomaterials
1996.
Vol.
17 No. 2
Optimal
224
algorithms is the difficulty of assessing how significant the internal parameters in the remodelling rule are. This is made more difficult by the large numbers of parameters in some models, and by remodelling rules which are very complicated. Although the actual remodelling process could indeed be as complicated as the rate equations postulated, the interpretation of several parameters in such a rule is very difficult. Thus, we have used a simple remodelling rule thus far, and we have concentrated on a detailed assessment of the internal parameters in the model. We have recently developed a specific mathematical connection between bone remodelling as a continuumlevel material adjustment to mechanical strain energy density and bone remodelling as an overall structural optimization processlO. This has led to simple conditions for a unique and stable solution for longterm bone density, for a class of hypothetical bone remodelling equations”. However, the physical significance of this connection is not obvious, and the usefulness of such a connection is not clear. In this paper we develop a global interpretation of the parameters in the local remodelling rule. This is very important for comparison with experimental results, because if we understand how the local parameters relate to global behaviour, then we can assess whether adjustments in local parameter values could account for differences between experimental results and model predictions. Understanding how local parameters influence results will also allow an informed assessment of how well the numerical model is behaving. Remodelling around intramedullary stems has been explored by Huiskes et al.‘” but the ‘internal remodelling’ algorithm displayed numerical difficulties. This study revisited the internal remodelling simulations of Huiskes et al., using a newly developed remodelling algorithm. This study also was intended to show a practical application of the remodelling algorithm which was developed in mathematical detail in our prior publications. The remodelling rule used in this study is similar to that of Huiskes et aI.” and Weinans et al:‘” but several important distinctions exist. The difference between this model and that of Huiskes is simply that this model arrives at unique, stable solutions which optimize a known cost function for the entire structure. The approach is new, because the exponent m was left as a parameter and chosen to make the global behaviour stable, whereas Huiskes et al. chose m = 1 so that the local parameters were well defined. The results are new, because prior to this work, remodelling approaches relied on global or local optimization, but there was no direct connection between the two approaches. The results are also new because this approach allows the solution of very refined finite element meshes in clinically relevant geometries, whereas if the model of Huiskes is refined, a continuous solution for the density distribution exists, but it is unstable and cannot be found in practice. The approach we have taken in our remodelling rule has been to studv the global behaviour of the remodelling simulaiion, and to derive limits on the Biomaterials
1996, Vol. 17 No. 2
remodelling
around
intramedullary
stems:
J.P. Harrigan
et al.
stimulus within the microstructure so that the solutions we arrive at are continuous and stable. There is no justification based on a remodelling signal within the microstructure. That is, we are letting the global behaviour of the model establish limits on the acceptable stimulus within the microstructure, and we are not being specific about what the stimulus means microscopically. The purpose of this study was to test how well a simple characterization of bone remodelling as an optimizing process can predict the behaviour of a large-scale finite element-based bone remodelling simulation. A secondary purpose was to gain a better understanding of the parameters in this and other finite element bone remodelling simulations.
METHODS Our bone remodelling theory is based on two descriptive relationships: a simple characterization of the bone remodelling stimulus and a structureproperty relationship for bone. Both are idealizations of known data, but they are reasonable approximations, and they are mathematically simple. In the future, more complicated remodelling rules and prosthesis geometries can be explored in detail based on the development here. This optimization study was used to derive a prediction of whether the cortical shell of bone will remain in the long term around a loaded intramedullary implant. Many studies have shown a drastic loss of bone around uncemented total hip replacements and porous-coated intramedullary rods, and fibrereinforced composite materials are proposed as a way to decrease the stiffness of a prosthesis enough to maintain bone. lising an optimization function developed in previous publications, we can make a simple prediction of the needed prosthesis size and material stiffness to maintain bone around the implant for axial loads, and transverse bending moments. The bone remodelling algorithm used in this study is described in our related publications1”.1’ and we will only summarize the basics of our methodology here. The method rests on a simple relationship between bone density and elastic modulus, and another simple relationship between strain energy density, volumetric density of hard tissue, and the rate of change of density. We assume isotropy, a material Poisson’s ratio that is constant, and a power law relationship between volumetric hard tissue density, +, and elastic modulus’“: E = E,, (I,”
(1)
with n the material exponent, usually taken as 2 (Rice et al.‘“) or 3 (Carter and Hayes’“). The bone remodelling rule we use is detailed in the literature, and it uses a remodelling stimulus which is strain energy density divided by 9”‘. Here m is an exponent which is chosen to make the simulations stable. This is a fundamental difference in approach from that of Huiskes et nl.‘“, Weinans et al.“’ and Beaupr6 et ~21.~.The difference between this stimulus and a set point is taken as proportional to the rate of
Optimal
remodelling
change of remodelling
around
intramedullary
4. This relationship rate equation which is
stems:
leads
to
et al.
a
a geometrical
bone
(21 with Q0 the stimulus set point. In a previous papeP, we showed that changing the state variable used in the remodelling equation allows that equation to be part of an optimizing process. By using density to the power m we can define a variable Y as li’ ZZ4”
(3)
and by changing
the rate equation
to be
we can satisfy the requirements for an optimizing function. This will have the same steady-state solutions as the previous equation. We have shown that this remodelling rule gives the same steady-state density distributions as Equation 2, and we have also shown that finding a distribution of y which makes the rate of change in Equation 4 zero is equivalent to finding stationary points in the function H =
(
:
>
225
T.P. Harrigan
HI + H,
with
optimization, rather than a density-based optimization, but the density-based optimization is actually general enough to include this type of optimization as a subset. In terms of the remodelling rule we have developed above, this analytical procedure corresponds to choosing a restricted class of density distributions in the optimizing process. That is, we use density distributions which are fully dense (i.e. >J= 1) over a bone area Ah and zero outside that area. This is a subset which is included in the set of all possible density distributions. Thus, minimizing H with respect to bone area fits within the framework developed for density-based remodelling if we consider only density distributions which are a restricted subset of all possible bone density distributions. Substituting for the strain in terms of the total force on the section yields H=
(m/n)F2L
(91
+\16AbL
2@,4, +&Ad
Since this is the strain energy only in a section of the structure, we cannot expect that minimizing H by manipulating Ab will give an exact answer, but since most femoral prostheses contain long stems, we expect this representative section to allow us to make an approximate generalization. Minimizing H with respect to Ah yields the following equation for Ah: Ab = (s)‘”
(&)
-ZAP
(10)
and for Ab to be greater than zero, we must have
(6) E,A,
e=1
H, is the total strain energy in the structure (bone and prosthesis included). In this study, we consider axial bending loads on a straight cylindrical prosthesis within a straight bone. We treat the loading cases separately for clarity, but the combined effect of these loads is easily incorporated.
(11)
This requirement can allow us to give a limit on the axial stiffness of an implant to maintain cortical bone around the shaft. In the simulations below, we used a fixed geometry, and we varied implant elastic modulus and the value of Q0 to assess cortical remodelling. The limit on the value of ‘16 needed for maintenance of a cortex is
(12) These requirements are reflected in the results of the simulations in which the cortex is allowed to remodel.
Axial loading For axial loads, if we consider a section of bone containing a prosthesis and we assume that the strain within the section is uniform, we can write Has H = ((m/n) (EPA, +Et,AtJ
mEb II2 __ 2n$ ( >
k2/2) + %AdL
(8)
where Ep and Eb are the elastic moduli of prosthesis and bone, respectively, A, and Ah are the crosssectional areas of prosthesis and bone, respectively, and L is the length of the section. We can now minimize the function H with respect to bone area, and arrive at a prediction of the stem stiffness required to maintain bone adjacent to the stem. On first glance, this analytical procedure seems to be
Transverse moments We developed two simplified models for bone adaptation to an intramedullary implant under a transverse bending moment in this study. Both models were aimed at assessing conditions for bone to remain around the implant. The models assumed that the bending radius of curvature for the implant and the surrounding bone were the same, and that bone exists in a concentric ring around an implant, with an interior radius ri and an exterior radius r, concentric with the implant section which is characterized by the flexural rigidity EpZ,,. The inner bone radius is not Biomaterials
1996.
Vol. 17 No. 2
226
Optimal
constrained to be the same as the outer prosthesis radius, since in many clinical cases the bone does not touch the implant in a given cross-section. Given this geometry and an equal radius of curvature for bone and implant, we can write the strain energy of the implant and prosthesis as TSE =
2M”L rcE,,(r; ~
$1 + 4EpIp
(13)
with M the transverse bending moment and EI, the bone elastic modulus. Both approximate approaches used this expression for total strain energy in the bone and implant. In the first bending optimization model, the density of the bone was taken as 1.0 and thus the restricted set of density distributions used in the analytical model consisted of circularly symmetric regions of fully dense bone, with the outer radius taken as fixed, and the inner radius variable. Thus, the index which selects the density distribution in the optimization model is ri, the inner bone radius. Minimizing the indicator functional over this restricted set of bone density distributions leads to the following equation: (14) and to assess when cortical bone will be adjacent to the implant, we need to assess the limit as ri approaches r,,. In that case, the bone flexural stiffness becomes small, and we can replace the first term on the right with f: and solve for c. For cortical bone to remain around the implant, we require r; < r, and this results in
(15) and we can take the limiting value E = 0. This results in a limiting value of flexural stiffness for an intramedullary implant for maintenance of bone, using the singleparameter restricted set of density distributions. In the second analytical optimization model we used the same geometric analysis but we took volumetric density as a fixed value less than unity. This resulted in the limit
This limit has some physical implications, since the exponent on y is negative for m > n as 7 goes to zero, the value of $‘I ~17):2~7~ will become infinite, but since the magnitude of the exponent is small, this limit can be accounted for in a practical sense. For n = 3.0 and m = 3.1, the exponent is equal to -0.016129 and if the actual density 4 = lo-” (below the lower bound for the numerical simulations) then yen “li!Zr” will equal only 1.58. This factor was used to predict the behaviour of the simulations.
Finite element
simulations
To implement the remodelling rule (Equation 4) we used a modification of the Euler backward timestepping method described in a prior publication”. Biomaterials
1996, Vol. 17 No. 2
remodelliw
around
intramedullary
stems:
T.P. Harrigan et al.
This time-stepping simulation methodology is based entirely on the remodelling rate equation, and does not depend on the optimization approach described in prior publications. The input data for these simulations were adapted from an input file for ADINA that was generated by ADINA-INIH.l”, the input processor for ADINA. A Euler backward time-stepping method, similar to a previously described method17, was implemented using Equation 4. Comparison of steady-state results between this method and the one used previously”” shows the results to be the same, within the numerical tolerances used for convergence. For axial loading simulations a two-dimensional axisymmetric model was used. The prosthesis diameter was 15 mm, with 5 mm of cancellous bone modelled between the outer surface of the implant and the cortex, which was 5 mm thick. Axial loads were applied to the top of the prosthesis, and the bottom of the bone was fixed. All interfaces were assumed to be perfectly bonded. We assumed that the distal tip of the intramedullary stem was spherical, and the distance from the distal tip of the prosthesis to the upper level of the femur was 107.5 mm. The prosthesis extended 10 mm beyond the end of the cortex. We chose to use a fine finite element discretization so we could assess the distribution of bone as predicted by the model in detail. Eight-noded elements, which assume that displacements vary quadratically with position, were used. The finite element mesh used 10 elements through the thickness of the cancellous bone and there were 100 rows of elements along the length of the prosthesis. There were 125 elements representing cancellous bone in the region of the prosthesis tip. The cortex was 5 elements thick and 120 elements long. To simulate the prosthesis, we used a T-element wide by l?O-element long shaft portion, and a 45-element tip. A total of 8287 nodes were used in the model. Convergence checks on the elasticity results of the remodelling simulations were not done, although similar studies in our prior publication’7 indicate that the mesh discretization is adequate for approximating the density. That is, we have not assessed the relationship between mesh size and the accuracy of the results. We have not completed a rigorous assessment of the various possible convergence criteria for the elasticity problem; this is under way. The results plotted in Figures I and z are not interpolated, however: they are made up of even blocks of colour which indicate element density. The rather smooth results seem to indicate that the mesh used in this study is sufficiently fine to answer the research question posed in this study. Also, since Huiskes et al.” obtained useful results for remodelling around an intramedullary stem with less than one-tenth the number of elements we LW: here, we feel our mesh was sufficient to characterize the elasticity problem. Ten axial load simulations were run in this study, with both the cortex and the cancellous bone allowed to remodel. Within the simulations, the prosthesis elastic moduli were chosen using a scaling procedure. We took the elastic modulus of fully dense bone as 17 GPa, and the prosthesis moduli to be 17, 34, 85, 136 and 204 GPa. We scaled the total load applied to
Optimal
remodelling
around
intramedullary
stems:
T.P. Harrigan
Figure 2 Remodelling results with $0 below the limits in Equation 72. In these cases, bone is maintained. The remodelling cases from Table 7 are (from left to right) A6, A7, A8, A9, and AlO. The central elastic moduli in the prostheses are, from left to right, 17, 34, 85, 136, and 206 GPa. As in figure 7, the shading indicates bone density, and each element is plotted as an even shade. Bone persists in all cases, with the initial cortical demarcation faintly apparent in cases A6, A9, and AlO.
Figure 1 Remodelling results with s above the limit in Equation 72. These results indicate that for @,, above the limit in Equation 6, wasting of the bone adjacent to the implant is to be expected. The remodelling cases from Table 7 are (from left to right) Al, A2, A3, A4, and A5 The elastic moduli for the prosthesis material were (from left to right) 17, 34, 85, 136, and 204 GPa, respectively. The shading shows bone density from 0.0 (darkest at right) to 1.0 (lightest at right). Notice that no interpolation is used; a constant density is displayed for each remodelling finite element. Resorption is clear in cases A2, A3, and A4, and partial in cases Al and A5
The finite element model was generated using ADINA-IN (ADINA Engineering, Watertown, MA, USA) and the finite element data were re-formatted for our simulation program using a translation program and a text editor. The axial and transverse load cases we used in this study were chosen to test the limits on prosthesis elastic modulus points and remodelling set simultaneously. For each prosthesis elastic modulus, a limiting value of q0 was calculated for axial loading from Equation 12. Six values of q0 were established so that we could test the limiting value of prosthesis stiffness, as predicted from Equation 12. Table 1 shows how the ten finite element runs in this study were used to test Equation 11, and Table 2 shows how eight of these ten cases were used to test Equation 12. Similarly, Equations 15 and 16 were used to generate limits on !QOand prosthesis stiffness which were tested using the seven large-scale three-dimensional cases shown in Table 3. Sets of three runs were arranged to test the two limits posed by Equations 15 and 16, with one case below the lower limit (on Ep or @J, one case between the two limits, and one case above both limits. Thus, by comparing cases B3, B2 and B6, or cases B5, B4 and B7, the results of progressive increases in Q0 for a given Ep were assessed, and by
1000 N, and we scaled the elastic moduli so that the ‘set point’ in our simulations was 0.5. Thus, the load which is relevant to the data in Table 1 is 1000 N, applied vertically as a point load to the top of the intramedullary stem. The cortical bone was fully dense (y = 1.0) and cancellous bone was simulated by taking y = 0.5. Seven three-dimensional finite element simulations were run in this study to assess the remodelling response of bone adjacent to an intramedullary stem which was loaded by a transverse bending moment. The finite element model consisted of a cylindrical stem with a spherical end-cap which was surrounded by bone with a cylindrical outer surface. The cylindrical portion of the stem was 100 mm long and the stem was 15 mm in diameter. The bone outer diameter was 35 mm. In the model, bone extended past the stem tip a distance of 10 mm and extended proximally to the top of the stem. The finite element model incorporated the entire stem and bone (i.e. a full 360” model was used) so that we could check for numerical problems by comparing left and right halves of the model.
Table 1 Elastic
17 34 85 136 204
Remodelling modulus
(GPa)
cases used to test predicted
limits on E,
Upper limit on Q0 from Equation (J mm3) 973.1 243.3 38.9 15.2 6.8
227
et al.
73
Simulation (case no.) 1055.8 422.3 84.5 33.8 8.4
(Al) (A2) (A3) (A4) (A5)
value above
Simulation (case no.) 422.3 84.5 33.8 8.4 4.2
value below
(A6) (A7) (A8) (A9) (AlO)
Biomaterials
1996, Vol. 17 No. 2
228
Optimal
Table 2
Remodelling
cases used to test predicted
Value of Q0 (J mm3)
Table 3
Q0 q0 q0 Q0
= = = =
Remodelling
500 J me3 1000 J mm3 2500 J mm3 5000 J me3
Equation 72
(GPa)
Simulation (case no)
25.8 57.7 91.2 182.5
17 34 85 136
cases used to test limits in pure bending
E, = 85 GPa
E, = 115 GPa
Bl
83 B2 B6
E, = 200 GPa B5 B4 87
comparing cases Bl, B2 and B7, the effect of increasing Ep for a given value of bone parameters was assessed. The elastic modulus of cancellous bone was taken as proportional to density cubed (n = 3), and the exponent on volumetric density that was used to divide the strain energy density (m in Equation 4) was taken as 3.1~. The density was limited to be no greater than 1.0 or less than 0.01 in the simulations. The norm of an increment in density within a given time step was computed by summing up the square of each density increment, dividing by the number of remodelling elements, and taking the square root of that quantity. The same was done to compute the norm of the remodelling rate and the norm of the overall density distribution. These norms included only the elements in which the density was greater than 0.01 and less than 1.0. Convergence to a solution at each time step was assumed when the norm of the density increment was 0.01 times the norm of the density. Convergence to a steady-state solution was assumed when the size of the remodelling rate (as scaled for simulation) was 0.005 times the size of the density. The finite element runs used in this study were coded in C and run on the Cray C-90 at the Pittsburgh Supercomputer Center. An adjustable time step value which increased in four intervals from 0.1 to 3.0 (in normalized time) was used to save on execution time. As a part of the numerical procedure, the stability of the simulations was checked at the end of the simulations by forming the matrix T + P, factorizing it numerically into LDLT form, and assessing whether negative diagonal entries resulted”. In all cases, the system of density evolution equations was found to be stable, as predicted analytically.
RESULTS The finite element remodelling results agreed in general with the predictions of the simplified analytical optimization results, and show that X& has the effect anticipated by the analytical predictions. The five axial load cases in which q’,, was chosen Biomaterials 1996,
Vol. 17 No. 2
around
intramedullary
stems:
T.P.
Harriman et al.
limits on q,,
Upper limit on E from
422.3 84.5 33.8 8.4
remodelling
(8) (7) (8) (9)
value below
Simulation (case no) 34 85 136 204
value above
(2) (3) (4) (5)
above the limit for cortical maintenance showed clear cortical resorption (Figure z), and the five axial load cases where Q0 was below the limit showed cortical maintenance (Figure 2). The implant stiffness increases from left to right in Figures 1 and 2. These plots show half of a section through the mid-plane of an intramedullary rod, and the colours are used to depict bone density. The colour scales on the right of Figures 1 and 2 show grey scales corresponding to densities that vary linearly from zero (at the bottom] to 1.0 (at the top). We have shaded the implant evenly. Cases Al and A5 (the farthest left and right in Figure 1) were only slightly over the threshold for resorption, based on Equations 11 and 12, so the persistence of some bone in the simulations is not surprising. Thus, the results of the large-scale simulations corresponded to the analytical results, and the set point Q(, has a substantial effect on whether bone will persist around an intramedullary stem. The results here also agree with the predictions of Equations 11 and 12. For a given value of $:, the cases noted in the latter part of Table 1 show cases which are above and below the limiting prosthesis stiffness in Equation 12. As in the cases which tested Equation 6, the remodelling results (Figures z and 2) support the derived limit on implant stiffness. The results of the bending simulations show that the predictions of Equations lli and 16 are borne out. As would be expected, the assumption that the density distributions are cylindrically symmetric is not accurate. Nevertheless, the density distributions followed the predictions of the limits over much of the prosthesis. Figures 3, 4 and 5 show remodelling cases used to test changes in density distributions due to changes in the remodelling set point and prosthesis elastic modulus. These figures are quarter sections of the density distribution around an intramedullary implant, and the direction of the applied bending moment is clear. In the leftmost cases, the prosthesis modulus and the remodelling set point are chosen so that both Equations 15 and 16 predict full resorption. The middle case in these figures has resorption predicted by Equation 15 but not Equation 16, and the rightmost case has bone maintenance predicted by both Equations 15 and 16. The distributions show that Equation 15 can predict when fully dense bone will occur around an implant, and Equation 26 can predict when bone will effectively disappear. Some of the features of the final density distributions have implications to follow-up studies of bone remodelling around intramedullary devices, and to finite element studies. The simulation cases in which the optimization approach appeared to be only marginally correct can lead to some mechanical
Optimal
remodelling
around
intramedullary
stems:
T.P. Harrigan
229
et al.
Figure 3 Three three-dimensional cases illustrating the changes in bone density for E, = 200 GPa and & = 2500 J mm3 (left, case B7), 1000 Jm3 (middle, case B4), and 500 Jmd3 (right, case 85). For this implant stiffness, the limits on Q0 from Equations 75 and 76, respectively, are 685.4 J mm3 and 1082.9 J mm3. Both analytical models predict resorption in case 87 and maintenance in case B5. The bone density in case 84 at the maximal radius was less than 1.0, indicating that Equation 75 predicts the occurrence of fully dense bone correctly. The correspondence between shading and bone density is as in Figures 7 and 2.
Figure 4 Three three-dimensional cases illustrating the changes in bone density for E, = 115 GPa and q0 = 5000 J mm3 (left, case B6), 2500 J m3 (middle, case B2), and 1000 Jmm3 (right, case 83). For this implant stiffness, the limits on Q0 from Equations 75 and 76, respectively, are 2174 J mm3 and 3435 J mm3. Both analytical models predict resorption in case B6 and maintenance in case B3. The bone density in case 82 at the maximal radius, along the shaft of the implant, was less than 1.0, indicating that Equation 75 predicts the occurrence of fully dense bone correctly. The correspondence between shading and bone density is as in Figures 7 and 2.
interpretations. Cases Al and A5 show density variations at steady state which extend from the proximal end of the intramedullary rod a substantial distance along the stem. These cases indicate that Saint Venant’s principle should be applied very carefully in bone remodelling simulations. These simulations indicate that in these cases, there is no ‘central region’ in which the load sharing between the bone and implant is uniform. In three of the cases in which bone persists at steady state (cases A8, A9 and AlO), the initial corticalcancellous interface can be seen, faintly. Previous simulations in which we used a more lax tolerance for convergence of the model to a steady state showed more pronounced residual cortical-cancellous demarcations, so we anticipate that still tighter convergence tolerances would further eliminate these demarcations. Also, in cases A6, A7 and A8, as well as in all cases where a pure bending load is applied, there are regions of low bone density adjacent to the prosthesis. These areas could appear as radiolucencies on a follow-up radiograph. In bending, the low strain in regions nearer to the neutral axis of the boneprosthesis composite beam is probably the cause for this effect. The low-density areas could be mistaken for regions where either bone ingrowth has not occurred, or where interface failure is occurring. That is, the low-density areas around a prosthesis seen on
X-ray could be a remodelling genuine interface failures.
effect
or they could
be
DISCUSSION The remodelling simulations carried out here were aimed at a restricted set of design parameters for intramedullary stems, and used a restricted set of physical effects. In this study, we concentrated on the limited goals set out so we could maintain control of the parameters in the simulation, and so we could assess the remodelling rule used in isolation. This allowed us to build upon detailed knowledge concerning the stability and the optimizing characteristics of the remodelling simulation used”. The results of this study show that a very coarse set of approximations can predict the behaviour of a largescale remodelling simulation around a loaded intramedullary implant. The use of very coarse approximations such as these is similar to the use of the Rayleigh Quotient in vibration analysis”‘. We have verified these predictions under torsional loads as well (unpublished data) and we are currently assessing the accuracy of simple optimization predictions in multiple loading situations. The accuracy of the approximate relationships derived from these coarse models seems good, and the relationships should give Biomaterials
1996, Vol. 17 No. 2
230
Optimal
Figure 5 Three three-dimensional cases illustrating the changes in bone density for q0 = 2500 J mm3 and with the implant elastic modulus taken as 200 GPa (left, case B7), 115 GPa (middle, case B2) and 85 GPa (right, case Bl). Both analytical models predict resorption in case 67 and maintenance in case Bl. The bone density in case B2 at the maximal Fadius was less than 1.0, indicating that Equation 75 predicts the occurrence of fully dense bone correctly. The correspondence between shading and bone density is as in Figures 7 and 2.
simple and useful data for the design of implants and as a guide for more accurate bone remodelling simulations. In its current form, our remodelling stimulus has only two parameters, and we can readily test whether these two parameters can characterize natural bone density distributions. Here we have developed some mechanical intuition for q’,, as used in Equation 4. This parameter is the set point of the remodelling stimulus, but it can also be interpreted as the weighting factor on bone mass (or a measure of bone mass) in a global optimization approach to bone remodelling. As is apparent from the remodelling rate equation, $, represents the amount of stimulation needed to maintain bone, and this study provides quantitative criteria for the intuitive idea that if the prosthesis is too stiff (or QO is too high) the bone will waste away around an intramedullary implant. We expect that in a individuals will have different biological setting, combinations of the variables in Equations 11, 12, 15 or ~6, and that there will be variable results around intramedullary stems, but there should be a value of 4, which most individuals do not exceed. This will give an upper limit on the prosthesis stiffness for bone maintenance. Although convergence checks were not made in this study (since they are not well developed for finite element problems such as these), we are confident that the models used are refined enough to assess whether the optimization approach is worthwhile. Previous Biomaterials
19%. Vol. 17 No. 2
remodelling
around
intramedullary
stems:
T.P.
Harriaan et al.
studies have shown that progressive refinement in these models results in closer and closer approach to a continuous solution, so that further refinement of the finite element meshes in this study will not change the salient results. The internal corner in the axially symmetric finite element mesh below the implant generates an artefact which will not be generally seen. However, since this study was concerned with a prediction of bone maintenance (or loss) adjacent to the cortex along the shaft of the stem, we feel that the presence of this artefact is not a serious problem in this study. Note that this artefact is due to the shape of the remodelling region, rather than the mesh discretization. The parameters we know of which are not in this model are as important for interpreting our results as those which are in the model. We assume that there are no non-mechanical stimuli which depend on position. That is, we assume that the non-mechanical effects are small, or that they are applied evenly throughout the modelled area. We do not assume that there is a non-responsive range around the @04353H.We have assessed the response to two isolated loading situations and we have not included the effect of strain rate in the model. Based on these results, we can assess the relative influence of individually added remodelling features on the overall response to implants in the future. The remodelling prediction we arrive at is similar to one derived by Huiskes et a1.12 for bone around an implant. Huiskes et aI.” used a remodelling theory based on strain energy density which modified outer cortical geometry in ‘external remodelling’ and which elastic modulus in ‘internal modified tissue remodelling’. This study, by comparison, is concerned only with internal remodelling. Here we confirm and extend the results of Huiskes et a1.l’ for internal remodelling. Huiskes et aI.‘” noted that “A homeostatic end configuration was reached only in the case of the fixed implant, the most flexible stem, and assuming the maximal threshold level s = 30(/o, in which only the proximal l/3 of the bone was fully demineralized”. This indicates that the method used by Huiskes et al. results in a very difficult numerical convergence problem. Since Huiskes et al. only reported one case in which the internal remodelling converged, the limit derived in that study for bone maintenance in internal bone remodelling could not be fully tested. That is, one case in which bone resorption was predicted and bone maintenance was another case in which predicted were not reported. Our method for internal remodelling results in a numerically stable formulation which allows a very fine finite element discretization without a lazy zone. The method is directly related to an optimizaand this study shows that the tion procedure, internal remodelling parameters can be interpreted physically in a practically useful model. Thus, we can confirm and extend the results of Huiskes et remodelling with a remodelling al.” for internal which has several theoretical and procedure practical advantages. The utility of a numerically robust finite element
Optimal
remodelling
around
intramedullary
stems:
T.P.
Harrigan et a/.
simulation methodology which allows continuous solutions of very refined finite element meshes is important, and the distinction between this approach and that of Weinans et al.‘” is noteworthy. As detailed in prior publications, if we choose m large enough, then a continuous, stable, unique solution is found to the posed mathematical problem for bone remodelling. As shown by Weinans et ~l.l”, taking m = 1 leads to a solution behaviour in which a continuous solution is approached, but then a discontinuous solution results. In fact, the continuous solutions approached by Weinans et al. are valid solutions for the remodelling response, but these solutions are unstable. If a numerical method was developed to arrive at those solutions, then the discussion of a discontinuous solution would be moot. The resulting solutions which diverge from the continuous solution are not unique, and the checkerboard patterns which result are entirely dependent on the element pattern in the finite element mesh. If the approach of Mullander et aI.” is taken, the patterns of divergence from the continuous solution are interesting, but the continuous solution still exists. Thus, the trabeculating patterns are interesting as possible trabecular patterns, but the background density distributions (i.e. the density as averaged over several trabeculae) are probably close to the continuous (but unstable) solution. As detailed in a prior publication16 the optimization function is derived from the remodelling stimulus, so that the predictions should be well represented. Note, however, that the finite element remodelling simulations do not depend on the optimization approach, they are based solely on Equation 4. Thus, the general agreement between the optimization predictions and the simulations seemed to confirm the usefulness of the optimization approach. This optimization approach is related to, but distinct from, other optimization approaches in which (a) the arrangement of a fixed mass of bone tissue is optimized to minimize strain energy in the structure or (b) weight is minimized for a given structural stiffness”. In the optimization scheme used in this study, a weighted sum of overall potential energy and a measure of the mass in the structure is used as the optimization function. This is similar to the method used by Subbarayan and Bartel’“. The distinction is primarily in the parameters used to define the problem. In approach (a), mass is fixed, and a Lagrange multiplier is used which corresponds roughly to ‘16 in the related continuum-level remodelling rule. In approach (b), the compliance is fixed, and a Lagrange multiplier is also used which is related to Q(, as well. However, in both (a) and (b) the Lagrange multipliers are unknowns in the optimization process. That is, the value of the Lagrange multiplier in either (a) or (b) is a result of the optimization, and is not specified at the outset. This corresponds to leaving QO as an unknown in the remodelling rate equation, and finding a value which satisfies the optimization criterion. Although, for a fortuitous choice of parameters, similar structures can be found in these approaches, the related bone remodelling algorithms are not the same.
231
ACKNOWLEDGEMENTS This research was supported in part by Pittsburgh Supercomputing Center grant Number 1 P41 RR06009 from the NIH National Center for Research Resources, and by Grant Ricerca Finalizzata of the Italian Ministry of Health, del. CIPE 13/l@% “Studio dei fattori the governano il fenomeno di adattamento funzionale de1 tessuto osseo ai carichi meccanicci”.
REFERENCES 1
Currey
2
Harrigan TP, Hamilton JJ. Bone strain sensation via transmembrane potential changes in surface osteoblasts: loading rate and microstructural implications. 1 Biomech 1993; 26: 183-200. Martin RB, Burr DB. Structure, Function, and Adaptation of Compact Bone. New York: Raven Press, 1989. Beaupre GS, Orr TE, Carter DR. An approach for timedependent bone modeling and remodeling - theoretical development. 1 Orthop Res 1990: 8: 651-661. Carter DR, Orr TE, Fyhrie DP. Relationships between loading history and femoral cancellous bone architecture. J Biomech 1989; 22: 231-244. Harrigan TP, Hamilton JJ. An analytical and numerical study of the stability of bone remodelling theories: dependance on microstructural stimulus. J Biomech 1992; 25(5): 477-488 (Erratum 26: 365-366). Hart RT, Davy DT. Theories of bone modelling and remodelling. In: Cowin SC, ed. Bone Mechanics. Boca Raton, FL: CRC Press, 1989; 253-277. Huiskes R, Weinans H, Reitbergen BV. The relationship between stress shielding and bone resorption around total hip stems and the effects of flexible materials.
3
4
5
6
7
a
JO. The Mechanical Adaptations of Bones. Princeton, NJ: Princeton University Press, 1984.
Clin Orthop Rel Res 1992; 274: 124. 9
10
11
Orr TE, Beaupre GS, Carter DR, Schurman DJ. Computer predictions of bone remodeling around porous-coated implants. J Arthroplasty 1990; 5(3): 191-200. Harrigan TP, Hamilton JJ. Bone remodeling and structural optimization. J Biomech 1994; 27(3): 323328. Harrigan TP, Hamilton JJ. Necessary and sufficient conditions for global stability and uniqueness in finite element simulations of adaptive bone remodeling. Int J
Solids Struct 1994; 31(l): 97-107. 12
Huiskes R, Weinans H, Grootenboer HJ, Dalstra M, Fudala B, Sloof TJ. Adaptive bone remodeling theory applied to prosthetic design analysis. J Biomech 1987;
13
Weinans adaptive
14
Carter DR, Hayes WC. The compressive behavior of bone as a two-phase porous structure. J Bone Joint Surg 1977; 59A: 954-962. Rice JC, Cowin SC, Bowman JA. On the dependence of the elasticity and strength of cancellous bone on apparent density. J Biomech 1988; 21: 155-168. Harrigan TP, Hamilton JJ. Optimality conditions for finite element simulation of adaptive bone remodelling. Int J Solids Struct 1992; 29(23): 2897-2906. Harrigan TP, Hamilton JJ. Finite element simulation of adaptive bone remodelling: a stability criterion and a time stepping method. Int J Numer Meth Engng 1993: 36: 837-854.
20: 1135-1150. H, Huiskes R, Grootenboer HJ. The behavior of bone remodeling simulation models. JBiomech
1992; 25: 1425-1441.
15
16
17
Biomaterials
1996, Vol. 17 No. 2
232
18
19 20
Optimal Adina R&D Inc. ADINA-IN for ADINA User’s Manual (Report ARD 87-4). ADINA R&D Inc., Watertown, MA, USA, 1987. Bathe KJ. Finite Element Procedures in Engineering Analysis. Englewood Cliffs, NJ: Prentice-Hall, 1982. Harrigan TP, Hamilton JJ. Bone remodeling adjacent to total hip replacements: a naturally occurring material design problem. I Comput Aided Mater Des 1993; 1: 27-40.
Biomaterials
1996, Vol. 17 No. 2
remodelling
21
22 23
around
intramedullary
stems:
T.P. Harrigan
et al.
Crandall SH, Karnopp DC, Kurtz EF, PridmoreBrown DC. Dynamics of Mechanical and Electromechanical Systems. New York: McGraw-Hill, 1968. Strang G, Kohn RV. Optimal design in elasticity and plasticity. Int J Num Mefh Engng 1986; 22: 183-188. Subbarayan G, Bartel DL. Bone remodeling around a hole: a comparison of theoretical and experimental results. Trans Orthop Res Sot 1991; 16:425.