Optimality conditions for finite element simulation of adaptive bone remodeling

OPTIMALITY

CONDITIONS FOR FINITE SIMULATION OF ADAPTIVE BONE REMODELING T.

Department

ELEMENT

P. HARRIGAS und J. J. HAWLTOS

of Ortho&dic Surgery. University of Missouri at Kansas City Medial 301 Holmes Sweet. Kanscls City. MO 63lOY. U.S.A.

School.

Abstract-Bone

remodeling in vcrtebratcs is widely quoted 3s LL process which optimizes the use of structural material. subject IY mechanical rcquircments. In vertchrutes, bone remodeling continues throughout life and tends IV preserve the structure of B particulur bone over decades of life. This implies that the processes which remodel bone arc stable. at least over ;I time period spanning many years. Many recent numeric;tl simulaGons of bone remodeling have used rate equations which have not been carefully assessed for stability and their ability lo produce an optimal sIructure. In this study. we w-state the conditions ncccssury for stability of ;1 particular hone remodeling rate equation derived in ;I related study. and wc investigate whether the rate equation used can produce an optimal structure. Within the context of a tinitc element discrctlzation, we show that thls rate cquatwn dws INN product a structure op~imizcd with respect to density. By making 21 smrplc niodilic;ition to lhe st;thlc remodeling r;iIc cqu;lti
The structure of bona within

the skclctal system of vortcbratcs has often been cited as an

sxamplc of ii natural syslcm in which lhc adi~ptation

of a tissue

to mechanical rcquircmcnts

has been achicvcd. kmc rcmodcling is ;I gcncral term which dcscribcs the proccsscs which maintain bone throughout the bone microstructure. within

adult life. Tissue

is laid down and rcsorbcd at surfaces within

since the rigid nuturc of the cxtraccllular matrix prcvcnts growth

the volume of the tissue. The biological principals which govern bone rcmodcling

arc currently

under ;I great deal ofdcbatc.

and numerical simulations

of bone rcmodcliny

show that the inllucncc of mechanical ctli~ts in rcmodsling arc also dillicult ability of thcsc simulations

to correctly mimic the rcmodrling

to ;ISS~XS. The

processes in bone is still being

assessed. and there is no rcmodcling rate equation which has been shown to be dclinitivcly bsttcr than others. Although

bone is constantly

being laid down and removed in order to

compcnsatc for fatigue damugu, no current thcorics for remodeling can account for this base-line ccllulur activity. The term “bone romodcling” in the shape of bone and changes in the distribution continually

maintaining and adapting the structure

is usually used to describe changes

of bone density with time. Thus,

of bones,

rcmovc fatigued matcriul and lay down new material within mechanical sophistication Wolff’s

GIII

bone illustrate

the dcgrcc of

occur in nature.

law is the most widely known theory for bone adaptation to stress. WollT’s

law is un extrapolation through

which

by

the biological proccsscs which

made by observing the correlation

the canccllous bone in iI number of joints

bctwccn patterns seen in a section

and the lines drawn in an analysis 01

stress using the mcthod of chaructoristics [XC Wollr ( 15X36)].The theory is ;L qualitative one which states that the internal and cxtornal structure of bone changes to adapt to changes in the loads applied to it. WolH”s law is often intcrprctcd today as a postulate that bone is distributed and oricntcd so as to support the loads placed on it with the minimum amount of material. This mutcrial optimization theory is very widely acccptcd as ;1 simple explanation for a complex bony structure.

but ;I rigorous

proof of the principal has not been

arrived at. A number of numcricnl models have been cxploitcd to study specific charactcrizations of Wolfi’s law rcccntly. Thcsc stud& can qualitatively predict the density

and J. J

T. P. H.WUGAS

z!w

distribution

of bone near joints.

qualitative

postulates studies

measures).

of this postulate

use continuum-averaged

of the mechanical

In this study.

quantities

questions

regarding

can certain

remodeling for a structure

within

First.

on the outside Cortical

classified

lead to an optimal

structure.

of bones.

and cancellous

(Reilly

and Burstein,

which

is very variable.

extent,

and the elastic

and orientation

1974). Cancellous The elastic modulus

(Carter

bone and cancellous elastic modulus

throughout

the formulation Poisson

this study

clTccts of’gcomctric

with

modulus

on rcmodoling

the information

avvilablc

nutrient

the concerted significant

action

but in this paper application equally

WC will

applied

to the entire

region

The form of the rrmodcling ofchange

of bone density

density

(Huiskes state”

remodeling

and

infbrmation

and Williams,

(such structure

as oxygen,

of a

loaded,

position

1973). Thus the local

hormonal

or

bone is produced

act microscopically.

Clearly,

yet containing factors.

This

by

there are

dense bone), restricts

the

efrccts arc either small or are

for bone thus far is fairly

is proportional

the bone density

to a “stimulus” at that

point.

CI ~1.. 1987; Orr er 01.. 1989) use a stimulus

can bc taken

the tissue. responding

being mod&d.

by bone density.

is a value of stimulus

stimulus

with a constant

of the tissue. WC arc

or physiological

to the mechanical

rules proposed

or a mcasurc of stress, divided

This “null

and

1977). In order to keep

bone seems to bc primarily

information

where non-mechanical

at a specilic point

to the stress state of the bone investigators

rcmodcl

that the overall

attention

to regions

cortical

density

bone. and WC have not asscsscd the

1987 ; Warwick

processes which

restrict

consider

between

bone as isotropic,

efrects. as in the skull (lightly

of this research

to some

of bone is abIt: to scnsc its anatomic

biochemical

of remodeling

non-mechanical

porosity

on the tissue porosity

to the density

is. no anatomic

(Netter,

This intlicatcs

joints.

well known

rcsponsc.

That

to the cells which

the available

concentrations).

heavily

and Hayes.

in cancellous

around

and has an elastic modulus

bone seem to act on a local basis within

lo that information

tissue stress and

found

law relationship

is related

which suggests that a small region

and act according

rate

these require-

is reasonably

bone varies with

(Carter

changes on the rcmoticlling

The cells which rcmodcl

porous.

canccllous

which

to local cuts from the micro-cnvironnicnl. is av;til;tblc

a power

by certain

rules for bone.

material

which

bone depends

WC have idealized

and an elastic

concentrating

the porous modulus

the range ofdensity

tractable,

ratio

We believe

1977 ; Rice cl trl.. 1988). Some authors

together.

to assess two stable’.’ Second,

on a remodeling

remodeling

of cortical

of cancellous

stress

bone. the dense layer of material

bone is highly

modulus

and Hayes, bone

bone,

and an elastic

and

as proposed

some requirements

into two types--cortical

of

these

the microstructure.

are these simulations

satisfies these two requirements.

bone has low porosity,

density

rate equations

ments can serve as checks on specific biologically-oriented Bone is generally

energy

remodeling

law? We have identified which

and a number

processes per se. they provide

which are important

these models.

rate postulates

Wolff’s

(strain

of the remodeling

we are using continuum-level

qualitative

interpretationsof

factors

Thus the original

have been confirmed.

are being assessed at this time. Although

and do not assess the details

an indication

equation

and the shape of long bone diaphyses.

made by the early anutomists

specific implementations structural

HNLTO~

at which

simple:

which is related

For example, which

and subtracted

;IS iI constant

value

throughout

some

is strain

energy

from a “null

state”.

the rate of change of density

a position-dcpcndcnt stimulus. A position-dcpcndent stimulus mechanical cffccts. An adjustable “null state” must bc carefully

the rate

is zero. This

the structure,

or as

is possible, due to nonused. however, since the

liberal USCof an adjustable parameter such as this can rcndcr the analysis meaningless. Careful study of the proccsscs used to simulate rcmodcling in these studies offers an opportunity temporal

to make the qualitative stability

of the numerical

postulates

of the early

bone remodeling

anatomists

simulations

more specific. The

has been discussed

in a

number of recent studies. Carter EI ai. (1989) discussed convergence issues, as did Fyhrie and Hollister (1990). Wcinans et al. (1989, 1990), Weinans (1991) and Huiskes et al. (1991) have explored two discrete

the stability remodeling

of bone remodeling elements,

Weinans

simulations

et al. (1990)

directly. indicated

Using a model that

the origin

with of the

unstable behavior was the remodeling rule. rather than the details of the finite element approximations. Huiskes ef al. (1991) have shown that the density distribution in bone can

Simulation be

simulated.

if the strain energy density

dependent “attractor

and Hamilton

the stability

from an intact bone is used as the position-

state” in a fairly coarse finite element mesh. The stability of this model

was considered tenuous. and refinement Hurrigan

2xw

of bone remodrling

of the finite element mesh was not attempted.

(1992a) have indicated from a continuous

of bone remodeling simulations

analytical model that

depends on the form of the remodeling rule.

The analytical model resulted in the derivation

of a stability

limit

for mechanical bone

remodeling postulates, and showed that the unstable behavior which resulted from failure to satisfy the limit was not similar

to the behavior seen itz ciao. Harripan

and Hamilton

(1992b) have developed a stability analysis for bone remodeling using finite element techniques. and have shown

that the requirement

for stability

as derived in the previous

analytical model is directly retlected in a finite element model. The stability of the algorithms used to simulate the temporal evolution of the density distribution to an effective Euler backward time-stepping procedure. This

were also studied, leading

study builds on the previous

finite element work. and develops ways in which the remodeling process developed can be considered as an optimal process. By optimal. we mean a process which manipulates state variables within

the structure in order to minimize a function of those state variables.

ANALYSIS

In order to assess the ways in which the remodeling rule used here can be regarded as an optimizing

process, we will first review the development of the stability conditions

finite element analysis. The conditions

using

for an optimal structure which are derived here are

based on a particular stimulus for bone remodeling. using strain energy density. Other bone remodeling theories, based on stress and strain measures, have been dcvelopcd. We chose to conccntratc on a rcmodcling stimulus cmaticillly simpler than stress

An incrcmcntal structure stability

based on strain cncrgy density. since it is math-

IllCilSlllW.

analysis has bcon clcvclopocl in or&r

being invcstigntrtl

to small perturbations

analysis. the structure

king

to study the rcuction of the

in density.

analyscd is nssumcd

b’or the purposes of the

to bc at a stats of rcmodcling

squilibrium

: i.c. WC:assume the clcnsity distribution has cvolvcd to a steady state situation. in this analysis are dt~clop~cl using it standard mcthocl for study of it dynamic system. with the rvontual tests for stability carried out using numerical methods. Stability

conditions

In this dcvclopnicnt. and the dcvclopmcnt to follow regarding optimalily

conditions,

the clcmcnt matrices arc written in global coordinates. The two indicts arc assumed to span the range

ofglobal

coordinate indicts, with non-zero entries only in the positions in which

the global coordinates couple to the clcmcnt. Thus we ncvcr refer to local coordinates for the clement stitTnrss matrices wc USC: hcrc. Also, the impliccl summation over repeated indices is only used when the rcpcatccl indices arc subscripts which follow displacements or stitfncss matrices. When summation over clemcnts is implied. it is explicitly cxprcssed using a summation sign, Wr will also use the subscripts /: g, p and (1to denote particular (global) degrees of freedom. The subscripts and Icft superscripts

c and s refer to element number.

When a quantity is a right superscript (used hcrc only on density). exponentiation Thus

WC have restricted our notation

in order to maintain

is meant.

clarity in the operations

we

intend. WC assume a material Poisson ratio which is indcpcndcnt of density. and it power law relationship

bctwcen density, (/I. and elastic modulus (Carter and Hayes. 1977). i.c.

E = E,rp. This

(1)

results in a material property matrix [C] which can be expressed as

[Cl = [C,]$‘“‘. with 4 the volumetric

(3

density, n the material exponent, and [C,] the material property

‘W)

T

matrix

for unit volumetric

element formulation

with the superscript

density (Bathe.

for element stiffness

19X2). Using

The average strain

matrix

matrix and

the element, “-f’, for a given displacement

energy density within

II is

with the exponent 111;t material constant. This continuum-lcvcl

strcsscs,

the microstructure. stimulus

in a standard finite

for unit density.

with 1:. the element volume. We assume the tissue stimulus

structure

this matrix

matrices yields

e denoting which element is used. IL,,, the element stiffness

*Klu the element stillness solution

P. H.\KKIGAN and J. J. H,\SIIL~OX

arc hi&r

strains

for bone remodeling is

equation is meant to retlcct the fact that the

and strain energy density are typici~lly magnified within

so that the stresscs,

StIXiIls

illld

than the contiliuuni-lcvcl

strain energy density within the micro-

variables. lkxupr~

ct trl. ( I990) proposed ;I

which accounts for the filet that the strength ofcn~~cllous

the square of the donsity of the tissue. so that the stress within

bone is proportional

to

the tissue is relittcd to the

ovcr:~ll stress in the invcrsc fashion. Since WCarc using strain energy density in this formuIiltion, antI sincc WCarc not making spcoific assumptions for rcmotlcling.

WC chose this relationship

ahout the spccilic physicill stimulus

with the saponcnt 111Icft ;IS it pitramctcr in the

relationship. Using this stimulus.

the overall rcniodcling rate equation is

l?iP*. = Y’, - ‘P,,. t?t Here the symbol Y,, denotes ;I “null point”, where the rcmodcling rate vanishes. In general. thcrc is a rate coclficirnt for simplicity

multiplying

the stimulus

(the right-hand sids of the equation), itnd

we have chosen to SGIIC time SO that the rate cocliicicnt is unity. The overill

rate equation can thus be written

as

(7fp< fp;” -‘%, “K,,,I~,, --. = __._.-.-._.._ ct 1V

yJ ,,.

(7)

“

The

finite clement equations to bc solved arc I<,,,,lr, = R,,. with K,,,. the ovcrilll stifTness R, the load vector for the matrix. formed by a sum of the slcmcnt stilrncss matrices. .lnd L linitc clcmcnt mesh. Taking perturbations of the linitc clcmcnt equations gives li,,,,Alr,, = - A$,,,LI,,,

with AK,,,, ths change in stilfncss for ;I perturbation with AcP, the change in donsity. ;IS

This

equation comes from taking derivitives

in density. This

(8) matrix can be written,

of the total stiffness matrix, written as a sum

Simulation

of individual

2901

of bone remodflinp

element stiffness matrices. with respect to density in each element. Thus,

change in displacements for a change in density distribution

the

can be written as

(10) Perturbations

in the rate equation. assuming a converged state. can be written as

= (d:“-“‘+(n-MI)&;-“‘-

“A~,)(~,‘K,,~,+ZII,.~K~~A\I~~+AII~K~~AL(,)

_ ~

zr; where we have used the fact that the element stiffness terms, neglecting higher order terms in perturbations.

0,

(II)

matrices are symmetric. Canceling and substituting

for Au gives

with

(13)

an d

(II - tll)ll,‘~K,,,ll* I),, = IS,.,

0,

zL’--- .’ 4,.

m

I)

(14)

9

c

with S,, equal to I il’c = s, and zero otherwise. HWXI CICC~Y

on ;I solution

which is cxponcntial

with time, and the equilibrium

in time, lhc solution

of qn

(12) will D,,

state arrived iit will be stiible, if the matrix M,,+

is negative dclinite. A matrix is negative definite if itll the cigenvalucs of the mntrix negative real parts. Since the rmttrix derived here is non-symmetric,

have

the eigenvalucs will be

complex numbers, ilnd special techniques for eigenvalue extraction are needed to compute them. Also. since the equations in perturbations separately. This

19921) is prcscrvecl. The ilniIlyticiIl the co&icients Harrigan

requirement

requirement that II be

(Harrigan

IL’SS thn

and Hamilton,

is necessary

III

to ktxp

in this diitgon:ll matrix Icss than zero. ilnd Hamilton

remodeling simulation non-symmetric

are linear, the matrix D, can be considered

shows that the i~nillyticitl stability

matrix.

(I99Zb)

have shown iI method whereby the stability

c;ln be itssessed without

of the

actually computing the eigenvalues of this

The procedure consists of factorizing

the symmetric

part of the

matrix into an LDLr form and assessing whether the diagonal matrix D contains negative elements (Crandall,

1955; Lancilstcr and Tismcnetsky.

1985; Wang.

1986).

An optimal structures approach to bone is useful for theoretical and practical reasons. The proposition that bone is iln optimal structure is one of the corncrstoncs of medical thought regarding bone. This proposition fits closely with other biological theories regarding natural selection for iidvuntagcous traits. If simulations of the natural rcmodcling situation which closely mimic itr rim

behavior arc found to optimize a given quantity

certain state variable. the theoretical framework

for these optimization

relative to a

postulates will be

strengthened. Also, if we can identify a particular state variable as being manipulated to optimize the structures studied. WC can indicate ways to frame or interpret studies of the

xl~

T

P.

HAKKIGAS

and

physical process within the microstructure. which optimizes

a function

J. J.

ttAMILT(IX

Practically.

a formulation

for bone remodeling

relative to a state variable is far easier to use than one which

does not. since the large literature practical and the theoretical

on optimization

can be used to sood effect. Both the

issues regarding bone remodeling are important.

conceptual issues are much debated. and simulations

since the

of bone remodeling often encounter

numerical problems. In order to assess whether these simulations some functional

equations. For example. Hamilton’s by minimizing

result

in a structure

which minimizes

on density, one can make an analogy to the vveak solution principle as applied to elasticity

an energy functional one can solve the elasticity field equations. Thus solving

the elasticity equations corresponds to finding a displacement distribution a functional

density distribution

minimizes

a functional on density if the density distribution

assessed whether a density distribution If a functional

which minimizes a functional

exists which is minimized

with respect to arbitrary

can be found by minimizing

admissible variations That

is, the first variation

equations. and the second variation

equilibrium

equations) should bc symmetric v&on

and results

function

then the rcmotlcling

with rcspcct to the two variations problem.

the i&as

ViIllIC.

ttwn

used.

can bc made mom

which mcasurcs some property

in a spccilic numerical

ilnd il’thc rcmoctcling rule used mininiizcs

in the functional will yield the

in the functional (the first variation in the

of the rcmodoling

concrctc. If thcrc is an indicator

the functional

in density. Also. the second variation in this

equilibrium

In ;I discrctizcd

a strain energy based “remodeling

on density.

at steady state by the remodeling rule. then

principal, a solution

functional should bc symmetric.

is given by

of the remodeling rate equation. That is. we have

can be considered to minimize

by analogy to Hamilton’s

distribution.

which minimizes

(stored elastic energy). In this study. we make an assessment of whether the

a long-term (i.e. “steady state”) solution stimulus”

of differential

problems shows that

of the density

wc can write that function as

(or lintls slablc points of) this indicator function,

rule can bc given by

(16) Thus

we can show that

JF(fp,) SCP, Thus,

a rcmodrling

(!?ff(tp,,

. . . , fp,,,

rl?f/(fp,. . . . .cP,,)

ctp, 2tp,

ccp,f?fP,

f'f-((P,) c'tp, .

rule can be dcrivcd from an overall indication

(17)

function on density if

and only if the dcrivativcs of that rule follow this equation. In our stability

stud&.

we have used a particular

rcmodcling rule in which

Ap, --,-. = I;‘(rp,).

(18)

i?r

and WC:have calculated the derivatives in cqn (I 7) for that rule, arriving at the matrices M,, and D,,. Since &I?, is not equal to M,,.. the dcrivativcs arc not equal for that particular rule. Thus, the mathematical form of that role is not the result of optimizing an indicator function on density. Careful study of eqn (I 3) shows that one possible way in which the matrix M,, can be taken as symmetric would be to take m = I. This would indeed satisfy the conditions for optimality, but as was shown in our previous studies (Harrigan and Hamilton, l992a, b), this choice would result in an unstable remodeling stimulus.

Thus.

it seems as though this

Simulation

of bone remodeling

1903

choice for remodeling stimulus can produce either a stable structure or an optimal structure. but not both. There is a way to arrive at a simulation

of bone remodeling

and which satisfies the mixed partial derivative condition.

which uses the above rule

however.

It involves modifying

the rate equation in such a way as to change the variable which is manipulated

to optimize

the function. Consider

the following

changes in the remodeling

rate equation.

Instead of using

density as our state variable. we can choose density to some positive power. as in

i’

=

(#I’.

(19)

with c a positive exponent, which we will arrive at later in this discussion. Note that when 7 is equal to zero. 4 is equal to zero. and when 7 is equal to one. 4 is equal to one. If we re-cast the rate equation as follows :

we will have the same equilibrium

states as before. since

(21) and the right-hand

side ofcqn

(30) is just the right-hand

side ofcqn

(7). multiplied

by each

clcmcnt volume. The situation

when (br approaches

zero needs consideration.

dcrivativc of?* will bc zero cvcn if the clcrivativc

ofC/J,. is not.

since at that point, the

From physical considerations,

this seems like a strength, rather than a watkncss. since the volumetric density in a particular region cannot bc less than zero. In gcncral. if the volumetric

density approaches zero in a

steady-state distribution.

approaches zero. Since II-~

the matrix product in cqn (20)

is negative for stability, density

goes

to

both 4: -“I and ry

aIs0

““.’ will approach

inlinity.

zero, thcsc two factors approach zero individually,

In a region where

and they

do

SO

in such

a way that their product is equal to ‘I-‘,,. This remodeling

rule will have ditferent temporal

behavior,

but the same equilibrium

states as eqn (7). as shown below. Since no current bone remodeling

rule purports to closely

mimic transient remodeling processes, we do not consider our modification

to be a serious

problem for the theoretical development. The left-hand side of eqn (20) could br troubling, WC have manipulated

since element volume appears there.

element volume here to remove it from the first term on the right-

hand side of eqn (20). Since the equilibrium

state for bone remodeling

has the right-hand

side of eqn (20) as zero within every element of the finite element model, multiplying right-hand

side ofeqn

(20) by element volume will not alter the equilibrium

since the stimulus is calculated the rate calculation

by element volume, and

divides this stimulus by the same clcmcnt volume, the effect of the finite

clcmcnt mesh paramctcrs expected

as the actual stimulus multiplied

will not intlucnce the rate calculation

from a typical dynamic

finite elcmcnt

model.

Thus.

any more than would be the change in temporal

behavior of the model between cqns (7) and (20) is related to calculating change of y. instead of the rate of change of I#. Although this remodeling rule seems artificial, minimizing

the

states. Also,

a functional

the the rate of

it can be shown to be the result of

on 7 if the exponent z is chosen equal to m. Since the remodeling

stimulus is the same in this equation as in the other equations shown, and since y is directly related to 4, we believe this result can have important implications for bone remodeling. In order to show that choosing z = m will satisfy the mixed partial derivative conditions, we must first assess the changes in displacements due to changes in density. The

T. P.

3wu

HAKKIGAN ;Ind J. 1. HASIILTOS

tinite element equations for equilibrium ~~lobai matrix 5

can be written as KC* =

R.and

if we expand the

K. we can write

or. written in terms of 7.

K,, = i

\vith “K,y the part of the stiffness simulation

i’:=eK,,+‘K,,.

(‘3)

I

f’

matrix

(e.g. an orthopedic implant).

which remains constant during Thus

the remodeling

we can take the derivative of KU

=

R

with

respect to ;‘,. and arrive at

(‘4) In order to choose z properly, we write the rate equation as

and wc take tlcrivativcs

as

(26)

Substituting

from cqn (14). we have

(27) The term 3;,, in this equation is equal to I if e = s and zero otherwise, so the first term ot the right of this equation is symmetric right-hand

in e and s by definition.

The second term on the

side of this equation contains a rather complicated matrix

be shown to bc symmetric

in the indices e and s, since

global linitc t’lcmcnt matrix.

'-'is A,,,

product which can

the inverse of the symmetric

In order to make the exponents on the terms yc and y, equal,

wc need only choose : = ttt. This

choice will

make eqn (27) symmetric

with respect to

indices c and s. and it indicates that eqn (17) is satisfied for this rate equation if we take dcriViltivcS with rcspcct

to y

instead Of 4.

DISCUSSION 13~ studying

the mathematical character of the remodeling

rate equations, we have

hecn able to assess two issues regarding bone remodeling as characterized by a remodeling rule based on strain cnorgy density : stability and optimality. In a related study (Harrigan and Hamilton. IO92b). WC used a finite clement discrctization to show an interrelationship bctwccn the bone rcmodcling process on a microscopic scale and the stability of the overall structure of bone on a global scale. This study confirmed stability requirements derived in ;I previous analytical study [Harrigan and Hamilton (1992a). and see after eqn (l4)]. We also arrived at a method to simulate temporal changes in bone density distributions within a tinitc clcmcnt analysis.

In this study, we build on the finite element study of remodeling

Simulation

1905

of bone remodeling

stability to study conditions for the presence of an optimal structure in bone. More specifically, we have arrived at a set of state variables, derived from a specific bone remodeling postulate. which can be considered to extremize a function of these variables over the entire structure. Although

we expect that the arguments put forward

finite element discretitation.

here do not depend on the

we have kept our development

within

this framework

for

simplicity in notation. A more general statement of the problem and results using functional analysis is probably possible. but is not attempted From this analysis of optimal choice for the remodeling optimizing

procedure.

been accomplished. and optimizing

conditions

state variable

A derivation

here.

for bone remodeling

which

of the actual functional

but an algorithm

to test candidate

rule part of an

which is minimized

bone remodeling

has not

rules for stability

behavior has been arrived at.

The function which the remodeling rateequation by integrating

we have arrived at a

makes the remodeling

the rate equation

appropriately,

We have found a method to differentiate on the density distribution,

minimizescan

theoretically

be derived

but we have been unable to do this as yet.

the complicated

dependence of the displacements

but we have not been able to integrate

the relationship

in a

sensible fashion as yet. The specific form of the rate equation in 7 which will result in an optimal structure can be related back to the rate equations in 4 by substituting

for 4’ for ;’ and by expanding

the

derivative with respect to time. This results in

(‘8)

which can bc written as

The form of this equation can bc seen as very close

to the form ofeqn (7), and this ccluation

is also very close to those studied in the litcraturc I991 ; Wrinans clement

illlOWS

(Huiskcs

rf (11.. 1987, 1991 ; Wcinans.

CI (11.. 1989, 1990). Since taking the time dcrivativc equal to zero for each Ollt2

10

factor

out tt1 ild

Cp“‘‘-” , the stable long-term

states for bone

rcmodcling will be the same in this analysis as with a strain energy density formulation as in Harrigan

and Hamilton

discussed aftrreqn

(199la).

The complication

such

as 4c becomes small. which was

(2 I), is still evident, only in eqn (29), the right-hand side is further divided

by 4’” - “q so that thr conditions necessary for the simulation to approach equilibrium regions with very low densities arc somewhat stricter than in eqn (7). The requircmcnt equilibrium

is now that the entire right-hand

fitShY than +L” _ I’. However.

this rcquircment

for for

side of eqn (20) approaches zero (with time) can be eliminated

by simply using eqn (20).

CONCLUSIONS

The optimizing standpoint.

nature of bone remodeling is important

Thcorctically.

bone has been postulated

from a thcorctical and practical

as an optimal

structure

century. and much ofthc biological and medical knowlcdgc ofboncadaptation rclativc to this assumption.

Practically,

an optimization

for over a

is intcrprstcd

proccdurc is much casicr to apply

than a point-by-point time stepping proccdurc. In a previous study (Harrigan and Hamilton. 19921). WC have shown a modification to existing bone remodeling algorithms which is ncccssary for stable time-stepping of thcsc equations.

In this study. WC show that by modifying

the state variable used to charactcrizc

the temporal rcsponsc of bone to strain energy density, we can product a simulation

which

satisfies the conditions for the existence of an overall cost function which is cxtrcmizcrl. We have derived a stable simulation of bone remodeling which is optimal in some sense. and which is similar to those in the literature

which appear to correctly mimic bone

2906

T. P. HARRIGAN and J. J. H~titrros

density distributions in a qualitative fashion. These simulations satisfy two qualitative mathematical requirements (stability and optimality) and are similar to simulations which qualitatively mimic bone density distributions. Further numerical analyses will be needed to assess whether the modified simuIations we arrive at here will correctly mimic the natural bone remodelting in a quantitative sense. REFERENCES Bathe. K. J. (t982). Fit&e Eiemetzr Procedwes In Engitreering Ann&six Prentice-Hall.Englewood Cliffs. NJ. Beat@, G. S...Orr. T. E. and Carter, D. R. (1990). An approach for time-de~ndent bone modeling and remodeling-theoretical development. 1. Orthopedic Res. 8.65 f -661. Carter. D. R. (1987).Mechanical loading history and skeletal biology. J. Biochem. 20. 1095-I 109. Carter, D. R. and Hayes. W. C. (1977). The compressive behavior of bone as a two-phase porous structure. J. Bone Joint Surgery S9A, 9511962. Carter. D. R.. Orr. T. E. and Fyhrie. D. P. (1989). Relationships between loading history and femoral cancellous bone architecture. J. &o&em. 22.3 1-244. Cowin, S. C. (1984). Mechanical m~eiling of the stress adaptation process in bone. C&if, 7%~. Inr. 36 (Suppl. If. S98-s 103. Crandall. S. H. (1955). Engineering Anulnis. McGraw-Hill. New York. Fyhrie. D. P. and Hollister. S. J. (1990). Comparison ofa trabecular tissue strain remodeling theory to experimental results. Tram. Orrhupedic Res. Sot. IS, 107. Harrigan. T. P. and Hamilton. J. J. (199 la). An analytical and numerical study of the stability of bone remodeling theories: Dependence on microstructural stimulus. 1. Bionrcch. ZS(S). 377488. Harrigan. T. P. and Hamilton. 3. J. (1991 bf. Fin&c clement simulation of adaptive bone remodeling: stability condi~i~~nsand a time stcppinp method. iptr. J. cVtmr. &fc*tlt. Enyit