JOURNAL OF THEORETICAL AND APPLIED MECHANICS 43, 3, pp. 555-574, Warsaw 2005 FUNCTIONAL ADAPTATION OF BONE AS AN OPTIMAL CONTROL PROBLEM Tomasz Lekszycki Institute of Fundamental Technological Research, Polish Academy of Sciences, Warsaw e-mail: tlekszyc@ippt.gov.pl The functional adaptation of bone is a process of bone tissue remodeling induced by variable in time mechanical demands that the skeleton has to satisfy. It is a very complex but highly organizedprocess composed of events atmicro-level (molecular and cellular) but having effects inmacro-scale (va- riation of bone internal structure and external shape).Mathematicalmodels of this phenomenon proposed in the literature represent formulas postulated on the basis of the results of medical observations or laboratory investiga- tions and describe locally the evolution of a material in space and time. In the present paper a use is made of the hypothesis of optimal response of bone, proposed earlier by the author,what enables derivation (instead of po- stulation) the remodeling rules from a very general and global assumption. It turns out that such a formulation has many similarities to engineering optimal control problems. The link between the postulated local adaptation rules and those derived from the global assumption is also discussed. Key words: adaptation, bone, objective functional, optimal control, remo- deling, tissue 1. Introduction Bones are complex biological structures which have the ability to adapt theirmicro-structures and external shapes to variable in time conditions. The- re are many factors of different nature, often related to each other already identified, as these which influence the remodeling process of bone tissue and evolution of bone shape. In spite of the fact that the complete understanding of mechanisms responsible for bone remodeling is not yet possible, it is well known that mechanical state of bone is one of the major factors controlling the changes in tissue structure. Such changes governed by time-variable me- chanical loads are usually called functional adaptation of bone. 556 T.Lekszycki Functional bone adaptation is a very complex process representing a chain of different interactions between the cells and intercellular matrix. Mathe- matical and computational modeling of this process is an important part of the research focused on investigation of the mechanisms responsible for tis- sue remodeling. It has also important practical aspects: in the future, reliable modelsmight be used in planning the orthopedic surgery and predicting bone reactions in different situations, for instance in the context of endoprosthe- sis durability, rehabilitation treatment or healing processes. It follows from recent investigations that bone remodeling is a highly organized process, the mechanisms of which are controlled at different levels. Many models have be- en proposed with postulated ”remodeling laws” – the mathematical relations describing evolution of bone micro-structure. On the other hand it can be expected that such an intelligent behavior of bone which is accomplished at micro-level (molecular and cellular) but has an effect at the macro-scale, sho- uld be ruled by some more general principle. Indeed, this is possible (under certain assumptions), to explain the adaptation of bone and the associated tissue remodeling as an issue that falls into a category of optimal control pro- blems. Moreover, some of the derived necessary optimality conditions can be interpreted as remodeling law. In fact in specific cases the derived formulas are similar to the already postulated remodeling laws. The fundamental as- sumption that is necessary to consider bone adaptation as an optimal control problem, is a hypothesis of optimal response proposed in Lekszycki (1999, 2002). According to this hypothesis, the bone reacts to variable in time and spacemechanical loads in an optimal way so that the rate of some global cost functional is extremized in order to assure the best possible improvements of bone configuration. In te following sections, the basicmechanisms responsible for bone remodeling are discussed and the ”optimal control approach” is pre- sented. To illustrate the links and differences between the ”local” postulated models of bone adaptation and these following from the ”global” approach presented here, a specific model is selected and briefly discussed. 2. Control of bone remodeling Themechano-sensory system inbone represents a chain of events, cell acti- vities andbiochemical processes related to each other and controlled according to some rule.Abasic understandingof the role of cells and thehypotheses con- cerning themechanisms responsible for functioning ofmechanosensory system in bone is necessary to propose the mathematical models of bone adaptation. Functional adaptation of bone... 557 2.1. Fundamental roles of cells responsible for functional boneadaptation One of the most important tasks of the skeleton is to support the mecha- nical loads associated with everyday activity and to protect internal organs. The ability to bear the extreme loads is to a large extent possible due to the bone ability to adapt to these functional demands by controlling its mass and morphology. The problem of bone adaptation has attracted researchers for more than one hundred years. Before the end of nineteenth century,Wolff for- mulated a statement generally knownas ”Wolff’s law”.According to it, ”every change in the form and function of bones or their function alone is followed by certain definite changes in their internal architecture, and equally definite secondary alterations in their external conformation, in accordance with ma- thematical laws”, see Wolff (1892). Since Wolff formulated this statement, a large numberof theoretical and experimentalworks have beenperformed.This includes observations at different levels ofmagnification, sometimes employing very sophisticatedmethods as well as more andmore advanced investigations of complex processes present in bone and responsible for its changes. Indeed, it follows from the results of the research that living bones are in continuous alteration – an activity which manifests itself in perpetual renovation of the bone ”material” and possibly, inmodification of itsmicro-structure and exter- nal shape.With improved experimental tools interesting resultswere achieved and new light was thrown into this subject. Nevertheless, in spite of all these efforts, not everything is understood completely and sometimes our knowledge is based on hypotheses that require more investigations. However some ideas have been already generally accepted. Bone tissue is a porous non-homogeneous and strongly anisotropic ma- terial undergoing continuous alteration due to complex biochemical proces- ses. It is composed of a solid matrix (build of mineral crystals and colla- gen fibers) and living cells – some of them buried in the matrix and others located in the pores filled with marrow. The turnover of bone is basical- ly associated with two simultaneous effects, bone formation and bone re- sorption implemented by specialized cells. These effects are closely coupled with each other in time and in space. They play a crucial role in mode- ling, maintenance, repair and aging of bones. One of the important factors at macro-level that contributes significantly to the local control of bone re- modeling is the variable in time mechanical loading determining strain di- stribution in bone. The process of bone remodeling is very slow compared to the mechanical loading changes. The remodeling effects can be observed at macro-level after days, weeks or evenmonths, while for the everyday activities and the associated mechanical loads, the time scale counted in seconds can 558 T.Lekszycki be used. Thus speaking about variable in time mechanical load we have in mind rather themean value (or amplitudes) and not the actual instantaneous values. This is a generally accepted concept that three families of cells are ma- inly involved in changes of bone micro-structure and evolution of the bone tissue itself. These cells are osteocytes, osteoblasts and osteoclasts. This fact has been already confirmed by many observations and the results of labora- tory investigations. The balance between the cells, their proliferation, diffe- rentiation and apoptosis are controlled both by local growth factors as well as by systemic hormones and plays a crucial role in the processes contribu- ting to bone turnover. This is more and more evident from the recent inve- stigations that osteocytes play a role of sensor cells controlling the process of bone formation and remodeling, while osteoblasts and osteoclasts are the ”actor” cells (Cowin and Moss-Salentijn, 1991; Ruimerman et al., 2005; Bur- ger and Klein-Nulend, 1999; Knothe et al., 2004; Klein-Nulend et al., 1995). Osteoblasts are the bone-forming cells. They produce new bonematrix in the regions where the additional material is needed to improve the bone perfor- mance. On the other hand, osteoclasts are responsible for tissue removal in the domains where material is not used. The balance between these two ef- fects, bone formation and bone resorption, is controlled by the complex and highly organized interactions between the cells and extracellular matrix. De- pending on the actual situation, the bone formation may be in balance, may exceed or be lower than the rate of bone resorption. Faster formation results in thicker and mechanically stronger bones. This may happen for example in response to increased mechanical loading. In another case, bone resorp- tion can be faster than bone formation what may be associated with bo- ne mechanical disuse or decease as osteoporosis and what results in loss of mass and deterioration of mechanical strength. The control of these effects is accomplished both via direct cell-to-cell signaling and via soluble mole- cules from the sensor cells (osteocytes) to the effector cells (osteoblasts or osteoclasts). The most numerous cells in a mature bone are osteocytes which play the role of sensors and represent probably the key element of the mechanosenso- ry system. They amount to about 90% of all cells in the bone tissue (Parfitt, 1977). Osteocytes have fixedpositions – they are buried in the bonematrix, do not divide and have long lifetime. Osteocytes have long processes (about 80) – ”fingers” which are connected with the processes of neighboring osteocytes. They form a network which is crucial in signal transmission to the actor cells. Osteoblasts are ”actor” cells which produce newbonematrix by collagen syn- Functional adaptation of bone... 559 thesis andmaking it calcify. After receiving the signals, osteoblasts located at the internal surface of the bone, that is the surface of the pores, build osteoid composed of collagen and other organic components. Osteoid is a ”brick” – an element of highly complex micro-architecture the tissue matrix is built of. The osteoblasts continue their activity and some of them, previously atta- ched to pore surface, become surrounded by the newly created material and transform into the osteocytes buried in ahardmatrix (Dudley andSpiro, 1961; Baud, 1968; Palumbo et al., 1990). During rapid growth, the proliferation of progenitor cells fromthemarrowpresent in the pores assures a necessary num- ber of new osteoblasts which replace the ones already buried in the matrix. But at a certain stage of this process, proper signaling slows down the pro- genitor cells proliferation and the remaining osteoblasts stop their production of osteoid while mineralization of the matrix continues. While the osteoblast transforms into osteocyte in the space called lacunae, its cellular volume de- creases and the collagen synthesis decays (Nefussi et al., 1991). Simultaneously the development of long cell processes with gap junctions starts (Doty, 1981). They are placed in channels called canaliculi. Thematrix around the osteocyte cells and their processes is not calcified, somechanically it ismoreflexible com- pared to the rigid calcified regions. This is an important observation in one of the proposed concepts trying to explain the mechanosensory functions of this complex system. It was alreadymentioned that the osteocytes are distributed in the three-dimensional space occupied bymatrix and that they form a com- plex network – they are connected with neighbors by cell processes and joined at gap junctions at their ends (Doty, 1981). Some of the osteocytes remain also in direct contact with other cells – osteoblasts and with the internal surface of the bone. This network can possibly form a system of mechanosensing and mechanotransduction in bone (Burger and Klein-Nulend, 1999; Klein-Nulend et al., 1995; Cowin and Moss-Salentijn, 1991; Weinbaum et al., 1994; Cowin et al., 1995). The third family of cells involved in bone remodeling are osteoclasts. They are responsible for bone resorption, at certain moments the signals appear in regions where the tissue should be removed and attracted osteoclasts are attached to the bone internal surface and dissolve the matrix. The relation between osteoclasts and osteoblasts activities determines the balance between the bone resorption and production. Therefore this is natural to expect that these processes are not independent of each other and are somehow controlled (among the others – by osteocytes activities). The ”coupling factor” discussed in Frost (1964) between osteoclast resorption and osteoblast formation has probably a mechanical nature (Rodan, 1991). 560 T.Lekszycki 2.2. Mechanosensory system in bones In modeling the bone adaptation, understanding of the mechanosensory systemandthemechanisms responsible for tissuevariations is of crucial impor- tance. In spite of numerous experimental data, the precisemechano-biological pathways of strain-induced bonemetabolism are not really known. Hence, the proposed theories are partly based on assumptions, see e.g. Ruimerman et al. (2005). However, some of the concepts win more and more prevalence in re- cent years and experimental data confirm their correctness. Since the Wolff law was formulated, the fact that bones are dependent onmechanical state is accepted and confirmed by many research results. Later investigations have proved that osteocytes – the sensor cells – are sensitive tomechanical loading and osteoblasts build the tissue while osteoclasts remove it.What exactly the osteocytes are sensitive to and how they ”fill” the mechanical situation and transformthe received information into the signals directed to osteoblasts, and osteoclasts is not clear as yet. Therefore inmanyworks an assumption ismade that the stimuluswhich excites osteocytes is proportional to a specificmeasure of strain or stress fields in the bone.One of themost commonly usedmeasures is the density of strain energy.Based on this assumption, severalmathematical and computational models of bone remodeling have been proposed. More than ten years ago, a hypothesis was presented by Cowin and co- workers according to which in intact bone the osteocytes are mechanically activated by flow of interstitial fluid through the lacunocanalicular system (Cowin andMoss-Salentijn, 1991; Weinbaum et al., 1994; Cowin et al., 1995). If this assumption is correct, the main stimulus for bone adaptation is the strain-driven motion of interstitial fluid through the canaliculi and along the osteocyte processes, which is somehow sensed and transduced by the osteocy- tes. Osteocytes then send signals to the actor cells, osteoblasts and osteoclasts and control their activities. The interstitial fluid flow in canaliculi can be in- deed the representative of gradients of the local strain fields since other voids in the tissue are about 103 times larger and their pressure is more uniform and almost equal to the blood pressure (Cowin and Weinbaum, 1995). On the other hand, according to one of the concepts, the osteoclasts are assumed to be recruited by osteocyte apoptosis due to micro-damage or cracks that damage the canalicular system important for nutrient transportation (Bronc- kers et al., 1996; Noble et al., 1997; Verborgt et al., 2000). They remove the cracked matrix while osteoblasts, due to increased level of strains associated with resorbing activity of osteoclasts, build the new one of improvedmechani- cal characteristics. This way the tissue is renewed and changes its anisotropic characteristics according to the mechanical state. Functional adaptation of bone... 561 3. Models of bone adaptation and the hypothesis of optimal response Since Wolff had formulated his famous statement, many attempts have been done to proposemathematical formulas – this ”mathematical law”men- tioned inhis statement –according towhich the configurations of bones evolve. Despite intensive research on this subject, there was no unanimity for many years concerning such problems as, for instance, what are the most impor- tant effects responsible for bone remodeling, what are the mechanosensory mechanisms including sensing of different signals and transmitting them to the effector cells, what are the mechanisms of bone maintenance, deposition and resorption and others, see e.g. Burr andMartin (1992), Cowin andMoss- Salentijn (1991), Harringan andHamilton (1993).Manymathematical models of bone adaptation based on different assumptions and taking into considera- tion diverse mechanical and non-mechanical effects have been proposed, see e.g. Carter and Orr (1992); Cowin (1995); Cowin and Hegedus (1976); Hart and Davy (1989); Hegedus and Cowin (1996); Levenston and Carter (1998); Luo et al. (1995); Prendergast et al. (1997); Prendergast and Taylor (1994); Taber (1995). Generally speaking, the models can be classified into three groups: bio- mechanical models, those based on structural optimization methods, and the models derived with the use of optimal response hypothesis. The first of the mentioned groups is the largest one. In most cases, the biological and medical observations and results of experiments are used to advance the hypothesis concerning possible causes of bone variations, theme- chanisms of stimulus sensing and signal transferring to the effector cells and the essence of tissue remodeling. Based on these hypotheses and the theoreti- cal investigations, the mathematical description of the adaptation process is postulated. Usually suchmodels are of a local nature in the sense that, accor- ding to an information about the actual local mechanical state of a bone, the local rate of tissue remodeling is calculated at the actual time. Such models can be verified using numerical computations and the results of clinical and experimental investigations. Some of the more recent works take into acco- unt the results discussed in the previous sections concerning the nature of the mechanosensory system in bone at the cellular level, see e.g. Ruimerman et al. (2005), Mullender and Huiskes (1995), Tezuka et al. (2005), Doblaré and Garcia (2002), Burger and Klein-Nulend (2003), Lemairea et al. (2004). The approaches based on the assumption that bones can be considered as optimal structures fall into the second group, see e.g. Bendsøe and Kikuchi (1988), Fernandes et al. (2000), Rodrigues et al. (1999), Folgado et al. (2004), 562 T.Lekszycki Tanaka andAdachi (1999). This assumption is considered controversial.More- over, such an approach does not enable any analysis of the alterations in bone due to conditions variable in time.On the other hand, it provides an asympto- tic solution, a possible bone configuration in equilibrium state i.e. under exter- nal loading constant in time, after a sufficiently long period of time. However one should be aware of the fact that the bone structure, even in equilibrium state, is not always optimal. The choice of the objective functional is arbitrary and represents an important step in this approach. It might be considered as a weak point of the procedure. Difficulties in including non-mechanical effects in the formulation represent an additional drawback. The last approachmentioned herewas proposed byLekszycki (1999, 2002) and is based on the formulated hypothesis of optimal response of a bone. It will be discussed in the next section. Many of the proposedmodels fall into the category for which the adapta- tion law can be symbolically written in the form dM dt =A(S(x, t)−S0(x, t)) (3.1) In the above relation M represents the bone mass in point x and time t, density of the material, Young’s modulus or any other parameter that is re- sponsible for local material characteristics. Of course, in the case when M represents mass, an additional relation is formulated – the relation between the local mass and the elastic characteristics of a material. S(x, t) is a sti- mulus – the signal that drives bone remodeling. It is often assumed that the stimulus is proportional to the density of strain energy, but also other solu- tions have been proposed depending of themechanosensory effects included in the analysis. S0(x, t) is a reference value of a stimulus, that is the value for which the bone is in a remodeling equilibrium state. This function should be assumed on the basis of other investigations (what is a week point of such an approach). In the following sections the relation (3.1) is compared with the derived adaptation law. It follows that it can be interpreted as a specific case of the relations obtained from the ”optimal control” approach based on the hypothesis of optimal response of the bone. 3.1. Hypothesis of optimal response and its relationwith optimal control It follows from the considerations concerningmechanosensory functions in bone presented in the previous sections that the adaptation process of bone can be considered as a class of optimal control problem. In spite of the fact that the formulation discussed here is not a typical approach known in the Functional adaptation of bone... 563 optimal control theory, many similarities appear what shows that the bone adaptation process can be indeed considered as a class of optimal control problem. To formulate the problem in mathematical terms, the hypothesis of optimal response of bone, formulated by Lekszycki (1999, 2002) is used. In many works an assumption was made that the bone represents an optimal configuration. Of course, if something is optimal or not is to a large extent a matter of optimization criterion. According to one of them, the object can represent an optimal solution, according to other it could be even the worst one. This assumption is based on the observation that an internal structure of bone is similar to optimal engineering structures, especially when some measures of strains or stresses are taken as the criterion with the constraint imposed for the overall mass, or opposite – when the mass is minimized with the constraints for maximum level of stress or strain measures. Thus it is probably not baseless. However, there are two important points that should be mentioned in this context. As it was already discussed in the previous sections, bone remodeling is an extremely complex phenomenon that depends on the processes of very diffe- rent levels of size, starting deep at molecular level with the results observable at macro-level. Additional effects of different nature such as biochemical, me- chanical, electrical and others are involved in its control and accomplishment. Many of these effects are closely related to each other and represent a com- plex control scheme, other are independent and work in parallel. Thus the ideal model should enable consideration of these effects and their possible interrelations. Anothermore important point is the fact that the optimal con- figuration – assuming that the criterion was correctly selected – represents some asymptotic, steady solution whichmight be achieved under the assump- tion that external and internal conditions that stimulate the changes in a bone are constant and do no vary in time. Of course, this is not the case in a real situation. We already know that a bone is exposed to conditions variable in time, bothmechanical and biological, which determine the processes involved in control andmaintenance. Therefore the optimal solutionmay provide some theoretical state which in fact can never be reached in practice. Nevertheless, in many cases the differences between this theoretical solution and the actual bone configurationmay be small or even negligible. Unfortunately, the ”opti- mization approach” provides only the final state under stimulation constant in time and do not enable to follow the remodeling process in time. These observationswere amotivation to propose a new approachwhich enables deri- vation (instead of postulation) the remodeling formulas including time effects (Lekszycki, 1999, 2002). This approach makes it possible to include in the formulation different effects as we learn more about the subject. 564 T.Lekszycki The starting point in the following considerations is the hypothesis of opti- mal response. According to it, the bone is not optimal but it reacts optimally to conditions variable in time. In other words, the bone attempts tomake the changes in its configuration within actual constraints to approach the best solution which is never achievable since it varies due to conditions variable in time. This way the bone is still in a state of pursuit of the optimal configura- tion. The general points of the approach are discussed below. Basic assumptions. Anassumption ismade that the considered effects are slow and the inertia terms are negligible. In addition, it is assumed that the theory of small displacements and velocities is valid. In order to describe the variations of bone, the control functions characterizing its structure should be selected µ(x, t). The derived remodeling law relates the velocities µ̇(x, t) to variable in time states of the bone. Criterion. In order to compare different bone structures, the functional G(µ(x, t)) is defined. It depends on a set of time-variable control functions determining the bone configuration. Greater/smaller value of this criterion means a better bone structure. The hypothesis of optimal response of a bone. According to this hy- pothesis, the bone reacts at each instant in an optimal way: that is, the rates µ̇(x, t) should assure the extremum of the objective functional. Theobjective functional. Theobjective functional results fromthe choice of criterion and the hypothesis of optimal response of the bone. It is assumed that it is represented by the rate of the criterion G. The global and local constraints. The constraints should be defined so as to take into account the important issues affecting the remodeling pro- cess. This is a crucial point in this formulation since different mechanical and non-mechanical effects can be included. This way, with growing knowledge concerning the mechanisms responsible for bone remodeling, an extension of actual models will be possible in future. Theadaptation law. Fromthe stationarity condition of the objective func- tional, with constraints attached to it by means of Lagrange multipliers, the optimality conditions follow. Some of them can be interpreted as the remode- ling law. Functional adaptation of bone... 565 4. Application of hypothesis of optimal response in derivation of bone adaptation law 4.1. General considerations A scheme of derivation of governing formulas is presented in this section. Let us introduce the following notation. C(µ)—tensor of material parame- ters where µ(x, t) is a control function defining the components of material tensor C (e.g. Young’s modulus in case of isotropic material or density ofma- terial in a more general case) and t denotes time (t is treated as a parameter – we consider only slow variations in time and inertia effects are neglected). As a result of this derivation, appropriate formulas are obtained for evolution in time of the function µ(x, t) following external conditions variable in time (e.g. mechanical loading or boundary conditions). Let U,W and V represent the set of kinematically admissible displacement fields, the set of kinemati- cally admissible variations of displacement fields and the set of kinematically admissible variations of velocities of displacement fields, respectively. Let us introduce also the following definitions a(u,v)= ∫ Ω C∇u ·∇v dΩ a′(u,v)= ∫ Ω Ċ∇u ·∇v dΩ l(v)= ∫ Ω b ·v dΩ+ ∫ Γf f ·v dΓ l′(v)= ∫ Ω ḃ ·v dΩ+ ∫ Γf ḟ ·v dΓ (4.1) where Ċ = dC dµ µ̇ ḃ= ∂b(x, t) ∂t ḟ = ∂f(x, t) ∂t µ̇= ∂µ(x, t) ∂t and Ωdenotes adomainoccupiedby thebodywhile Γf is apart of aboundary surface where the loading is defined.We can now express the potential energy as Π(u)= 1 2 a(u,u)− l(u) u∈U (4.2) and its time derivative as Π̇(u,u̇)= dΠ dt = 1 2 a′(u,u)+a(u, u̇)− l(u̇)− l′(u) (4.3) It is easy to check that the stationarity conditions of the functional (4.3) with respect to independent variations of u and u̇ are satisfied δu̇Π̇(u,u̇)= 0 δuΠ̇(u,u̇)= 0 (4.4) 566 T.Lekszycki Thus the weak formulation of the analysis problem is a(u,δu̇)− l(δu̇)= 0 ∀ δu̇∈V a′(u,δu)+a(u̇,δu)− l′(δu)= 0 ∀ δu∈W (4.5) Let us now define the comparison functional which represents a measure ne- eded to compare different systems (bones) G= ∫ Ω S(u,µ) dΩ (4.6) According to the hypothesis of optimal response, the cost functional is defined as (we assume that the domain Ω does not evolve in time), Ψ = dG dt = ∫ Ω Ṡ dΩ (4.7) Let us define now the optimization problem min µ̇ Ψ(u,u̇, µ̇) (4.8) with additional global and local constraints applied a(u,δu̇)− l(δu̇)= 0 ∀ δu̇∈V a′(u,δu)+a(u̇,δu)− l′(δu)= 0 ∀ δu∈W ∫ Ω hi(µ̇(x, t)) dΩ−Ai(t)= 0 i=1, . . . ,Ng gi(µ̇(x, t))­ 0 i=1, . . . ,Nl (4.9) where the following notation has been introduced: gi(·) and hi(·) – local and global constraints defined for the control function µ̇, Ai is a limit imposed globally on the control function, Nl and Ng denote the numbers of local and global constraints, respectively. Let us build an extended cost functional bymeans of the Lagrange multi- pliers λ1, λ2, ρi, ηi, and slack variables αi L(u,u̇,ua1,u a 2, µ̇,ρi,ηi,αi)= =Ψ(u,u̇, µ̇)−a(u,ua2)+ l(u a 2)−a ′(u,ua1)−a(u̇,u a 1)+ l ′(ua1)+ (4.10) + Ng∑ i=1 ρi(t) [∫ Ω hi(µ̇(x, t)) dΩ−Ai(t) ] + Nl∑ i=1 ∫ Ω ηi(x, t) [ gi(µ̇(x, t))−α 2 i(x, t) ] dΩ Functional adaptation of bone... 567 where additional functions ua1(x, t) and u a 2(x, t) are defined using Lagrange multipliers λ1 and λ2. They represent state variables of the so-called adjoint system ua1 =−λ1δu u a 2 =−λ2δu̇ (4.11) Comment. In a general case of an arbitrary comparison functional, an ad- ditional system called ”adjoint system” appears. This results in the necessity of analysis of this system since the adaptation relations are expressed in terms of state variables of both the primary and adjoint systems. 4.2. Examplary derivation of the adaptation law for a specific case In a specific case when the comparison functional G represents the global measure of stiffness, the situation is simplerbecause both systems, theprimary and the adjoint ones are equal to each other S = 1 2 C∇u ·∇u G= ∫ Ω S dΩ (4.12) Ψ = 1 2 ∫ Ω (Ċ∇u ·∇u+2C∇u̇ ·∇u) dΩ Let us apply the specific constraints to derive a sample adaptation law a(u,δu̇)− l(δu̇)= 0 ∀ δu̇∈V a′(u,δu)+a(u̇,δu)− l′(δu)= 0 ∀ δu∈W (4.13) ∫ Ω µ̇(x, t) dΩ−A0(t)=h1(µ̇)= 0 (4.14) ∫ Ω µ̇2(x, t) dΩ−B0(t)=h2(µ̇)= 0 µ̇(x, t)− µ̇max(x, t)= g1(µ̇,x)¬ 0 −µ̇(x, t)+ µ̇min(x, t)= g2(µ̇,x)¬ 0 (4.15) −µ̇(x, t)H(µmin(x, t)+θ−µ(x, t))= g3(µ̇,x)¬ 0 µ̇(x, t)H(µ(x, t)−µmax(x, t)+θ)= g4(µ̇,x)¬ 0 568 T.Lekszycki where the following notation hasbeenused, µmin,µmax, µ̇min, µ̇max –minimal andmaximal values of the function µ and its velocity, respectively, A0, B0 – global constraints imposedon µ̇.These functions shouldbedefinedon thebasis of experimental results and clinical observations. H(·) denotes Heaviside’s function. θ represents small neighbourhood of the limit values. According to the last two constraints, the function µ in the neighbourhood close to µmin can not decrease and when it is close to µmax, it can not grow. Let us build an extended cost functional by means of Lagrange multipliers λ1, λ2, ρ1, ρ2, η1, η2, η3, η4 and slack variables α1, α2, β1, β2. Then the cost functional has a form L(u,u̇,ua1,u a 2, µ̇,ρ1,ρ2,η1,η2,η3,η4,α1,α2,β1,β2)= = 1 2 a′(u,u)+a(u, u̇)−a(u,ua2)+ l(u a 2)−a ′(u,ua1)−a(u̇,u a 1)+ l ′(ua1)+ +ρ1(t) [∫ Ω µ̇(x, t) dΩ−A0(t) ] +ρ2(t) [∫ Ω µ̇2(x, t) dΩ−B0(t) ] + + ∫ Ω η1(x, t) [ µ̇(x, t)− µ̂max(x, t)+α 2 1(x, t) ] dΩ+ (4.16) + ∫ Ω η2(x, t) [ µ̇(x, t)− µ̂min(x, t)−α 2 2(x, t) ] dΩ+ + ∫ Ω η3(x, t) [ µ̇(x, t)H(µmin(x, t)+θ−µ(x, t))−β 2 1(x, t) ] dΩ+ ∫ Ω η4(x, t) [ µ̇(x, t)H(µ(x, t)−µmax(x, t)+θ)−β 2 2(x, t) ] dΩ From the stationarity condition of the cost functional we obtain: • state equations for the primary system, • state equations for the adjoint system, • set of applied constraints, • set of equations for Lagrange multipliers and slack variables, • adaptation rule. For the specific case considered here (global compliance as a comparison func- tional) we have ua1(x, t)=u(x, t) u a 2(x, t)= u̇(x, t) (4.17) Functional adaptation of bone... 569 and the remodeling equations, representing one of the necessary conditions for stationarity of the cost functional can be expressed as follows µ̇(x, t) = − 1 2ρ2(t) [1 2 ∂Cijkl ∂µ eijekl+ρ1(t) ] − η1(x, t) 2ρ2(t) − η2(x, t) 2ρ2(t) + (4.18) − η3(x, t)H(x, t) 2ρ2(t) − η4(x, t)H(x, t) 2ρ2(t) Lagrangemultiplier ρ1 canbe interpreted as a variable in time reference value, often used in the ”postulated” models. A numerical example of application of this adaptation law was calcula- ted. In Fig.1 the results of simulation are presented. At the initial state a homogeneous material was used. After application of mechanical load the mi- crostructure has developed. This microstructure was rearranged due to remo- deling after endoprosthesis implantation. The micro-structures presented in this figure are similar to the clinical observations of the real bones. 4.3. Postulated and derived models In the previous section, a simple adaptation model was mentioned (3.1). This formulahasa similar formto thederivedone (4.18). It canbenoticed that remodeling according to (4.18) is also proportional to the difference between some stimulus and its reference value, similarly as it was assumed in (3.1). But in the present case the reference value is not assumed and is not constant. It follows from the analysis and is dependent on the Lagrange multipliers that have to be determined during the analysis. Therefore such a model is more realistic. On the other hand, the ”optimal control” approach enables the control of the total amount of material via global and local constraints, what can be useful in many situations as, for example, in cases of osteoporosis or bonegrowth.Thepostulatedmodel enables also the control of thematerial but indirectly – bymodification of the reference value of the stimulus. In practical considerations it is very difficult since the reference value is not known. 5. Conclusions In the present paper the idea has been presented that the functional ad- aptation of bone can be considered as a modified optimal control problem. Such a formulation is associated with significant advantages. However, the 570 T.Lekszycki Fig. 1. The effect of a numerical simulation of bone adaptation before and after endoprosthesis implantation (from the left to right: initial state, adapted configuration after application of the mechanical load, remodeled configuration after endoprosthesis implantation). The structure of a real bone is presented below formulation should be very carefully proposed to enable consideration of im- portant biological effects that are involved in the control process. Among the others, the interactions between different families of cells should be included. An attempt to do so was alreadymade by Lekszycki (2002), but an improved formulation is still necessary. Acknowledgement This work was supported by the Polish State Committee for Scientific Research, grant KBN 3T11F00727. Functional adaptation of bone... 571 References 1. BaudC.A., 1968, Submicroscopic structure and functional aspects of the oste- ocyte,Clinical Orthopaedics and Related Research, 56, 227-236 2. Bendsøe M.P., Kikuchi N., 1988,Generating optimal topologies in structu- ral designusing ahomogenisationmethod,Comput.Methods Appl.Mech. Eng., 71, 197-224 3. Bronckers A.L.J.J., Goei W., Luo G., Karsenty G., D’Souza R.N., Lyaruu D.M., Burger E.H., 1996, DNA fragmentation during bone forma- tion in neonatal rodents assessed by transferasemediated end labeling, Journal of Bone and Mineral Research, 11, 1281-1291 4. Burr D.B., Martin R.B., 1992,Mechanisms of bone adaptation to the me- chanical environment,Triangle, 31, 2/3, 59-76 5. Burger E.H., Klein-Nulend J., 1999,Mechanotransduction in bone – role of the lacuno-canalicular network,FASEB J., 13S, S101-112 6. Burger E.H., Klein-Nulend J., Smit T.H., 2003, Strain-derived canalicu- lar fluid flow regulates osteoclast activity in a remodelling osteon – a proposal, J. Biomechanics, 36, 1453-1459 7. Carter D.R., Orr T.E., 1992, Skeletal development and bone functional adaptation, J. of Bone and Mineral Research, 7, S389-S395 8. Cowin S.C., 1995,On theminimization andmaximization of the strain energy density in cortical bone tissue, J. of Biomech., 28, 4, 445-447 9. Cowin S.C., Hegedus D.H., 1976, Bone remodeling i: theory of adaptive elasticity, J. Elascity, 6, 3, 313-326 10. Cowin S.C.,Moss-Salentijn L.,MossM.L., 1991,Candidates for theme- chanosensory system in bone, J. Biomech. Eng., 113, 191-197 11. Cowin S.C., Weinbaum S., 1998, Strain amplification in the bonemechano- sensory system,Am. J. Med. Sci., 316, 184-188 12. Cowin S.C.,Weinbaum S., ZengY., 1995,A case for bone canaliculi as the anatomical site of strain generated potentials, J. Biomech., 28, 1281-1296 13. Doblaré M., Garcia J.M., 2002, Anisotropic bone remodelingmodel based on a continuum damage-repair theory, J. Biomechanics, 35, 1-17 14. Doty S.B., 1981,Morphological evidence of gap junctions between bone cells, Calcif. Tissue Int., 33, 509-512 572 T.Lekszycki 15. Dudley H.R., Spiro D., 1961, The fine structure of bone cells,The Journal of Biophysical and Biochemical Cytology, 11, 627-649 16. Fernandes P., Rodrigues H.C., Jacobs C.R., 2000, Amodel of bone ad- aptationusingaglobaloptimization criterionbasedon the trajectorial theoryof Wolff,Mechanics in Biology, ASME 2000, AMD-Vol.242/BED-Vol.46, 173-184 17. Folgado J., FernandesP.R.,Guedes J.M.,RodriguesH.C., 2004,Eva- luation of osteoporotic bone quality by a computational model for bone remo- deling,Computers and Structures, 82, 1381-1388 18. FrostH.M., 1964,Dynamics of bone remodeling, In: FrostH.M. (Edit.),Bone Biodynamics, Little, Brown, Boston,MA, 315-333 19. HarringanT.P.,Hamilton J.J., 1993,Bone strain sensationvia transmem- brane potential changes in surface osteoblasts: loading rate andmicrostructural implications, J. of Biomech., 26, 183-200 20. Hart R.T., DavyD.T., 1989, Theories of bonemodeling and remodeling, In S.C. Cowin (Edit.),Bone Mechanics, CRCPress, Boca Raton, FL, 253-277 21. Hegedus D.H., Cowin S.C., 1996, Bone remodeling ii: small strain adaptive elasticity, J. of Elasticity, 6, 337-355 22. Klein-Nulend J., Van der Plas A., Semeins C.M., Ajubi N.E., Fran- gos J.A., Nijweide P.J., Burger E.H., 1995, Sensitivity of osteocytes to biomechanical stress in vitro, FASEB J., 9, 441-445 23. KnotheT.M.L.,Adamson J.R.,TamiA.E.,BauerT.W., 2004,The oste- ocyte,The International Journal of Biochemistry and Cell Biology, 36, 1-8 24. Lekszycki T., 1999, Optimality conditions in modeling of bone adaptation phenomenon, J. Theoret. Appl. Mech., 37, 3, 607-623 25. Lekszycki T., 2002, Modelling of bone adaptation based on an optimal re- sponse hypothesis,Meccanica, 37, 343-354 26. Lemairea V., Tobina F.L., Grellera L.D., Choa C.R., Suvab L.J., 2004,Modeling the interactions between osteoblast and osteoclast activities in bone remodeling, J. of Theoretical Biology, 229, 293-309 27. Levenston M.E., Carter D.R., 1998, An energy dissipation-based model for damage stimulated bone adaptation, J. of Biomech., 31, 7, 579-586 28. LuoG.,Cowin S.C., SadeghA.M.,ArramonY.P., 1995, Implementation of strain rate as a bone remodeling stimulus, J. of Biomechanical Engineering, 117, 3, 329-338 Functional adaptation of bone... 573 29. Mullender M.G., Huiskes R., 1995, A proposal for the regulatory mecha- nism ofWolff’s law, J. of Orthopaedic Research, 13, 503-512 30. Nefussi J.R., SautierJ.M.,NicolasV.,ForestN., 1991,Howosteoblasts become osteocytes: a decreasing matrix forming process, J. Biol. Buccale, 19, 75-82 31. Noble B.S., Stevens H., Loveridge N., Reeve J., 1997, Identification of apoptotic changes in osteocytes in normal and pathological human bone,Bone, 20, 182-273 32. PalumboC.S.,Palazzini,ZaffeD.,MarottiG., 1990,Osteocytedifferen- tiation in the tibia of newborn rabbit: an ultrastructural study of the formation of cytoplasmic processes,Acta Anat. (Basel), 137, 350-358 33. Parfitt A.M., 1977, The cellular basis of bone turnover and bone loss: a rebuttal of the osteocytic resorption –Bone flow theory,Clin. Orthop., 236-247 34. Prendergast P.J., Huiskes R., Søballe K., 1997, ESB Research Award 1996. Biophysical stimuli on cells during tissue differentiation at implant inter- faces, J. of Biomechanics, 30, 6, 539-548 35. Prendergast P.J., Taylor D., 1994, Prediction of bone adaptation using damage accumulation, J. Biomech., 27, 1067-1076 36. Rodan G.A., 1991, Mechanical loading, estrogen deficiency, and coupling of bone formation to bone resorption, Journal of Bone and Mineral Research, 6, 527-530 37. RodriguesH.C., JacobsC.R.,GuedesJ.M.,BendsøeM.P., 1999,Global and localmaterial optimizationmodels applied to anisotropic bone adaptation, In P. Pedersen and M.P. Bendsøe (Edit.), Synthesis in Bio Solid Mechanics, Kluwer Academic Publishers, 209-220 38. Ruimerman R., Hilbers P., Van Rietbergen B., Huiskes R., 2005, A theoretical framework for strain-related trabecular bone maintenance and ad- aptation, J. of Biomechanics, 38, 931-941 39. Taber L.A., 1995, Biomechanics of growth, remodeling, and morphogenesis, Appl. Mech. Rev., 48, 8, 487-545 40. Tanaka N., Adachi T., 1999, Lattice continuum model for bone remode- ling considering microstructural optimality of trabecular architecture, In P. Pedersen andM.P. Bendsøe (Edit.), Synthesis in Bio Solid Mechanics, Kluwer Academic Publishers, 43-54 41. Tezuka K., Wada Y., Takahashi A., Kikuchi M., 2005, Computer- simulated bone architecture in a simple bone-remodeling model based on a reaction-diffusion system, J. Bone Miner. Metab., 23, 1-7 574 T.Lekszycki 42. Verborgt O.G.J., Gibson, Schaffler M.B., 2000, Loss of osteocyte in- tegrity in association with microdamage and bone remodeling after fatigue in vivo, J. of Bone and Mineral Research, 15, 60-67 43. Weinbaum S., Cowin S.C., Zeng Y., 1994, A model for the excitation of osteocytesbymechanical loading-inducedbonefluid shear stresses,J.Biomech., 27, 339-360 44. Wolff J., 1892, Das Gesetz der Transformation der Knochen, A. Hirchwild, Berlin (Translated by Maquet P., Furlong R., The Law of Bone Remodeling, 1986, Springer, Berlin) Funkcjonalna adaptacja kości jako zagadnienie optymalnego sterowania Streszczenie Funkcjonalna adaptacja kości jest procesem polegającym na przebudowie tkanki kostnejwywołanej zmieniającymi sięwczasiewymaganiamimechanicz- nymi, jakiemusi spełniać szkielet kostny. Proces ten jest niezwykle złożony, ale doskonale zorganizowany i składa się z szeregu zjawisk zachodzących na po- ziomie mikro (molekularnym i komórkowym) lecz mających efek na poziomie makro (zmiana zewnętrznego kształtu kości oraz jej struktury wewnętrznej). Matematycznemodele tego zjawiska, postulowanewoparciu o obserwacjeme- dyczne i badania laboratoryjne, opisują lokalną ewolucję materiału w czasie i przestrzeni. W tej pracy zastosowano hipotezę optymalnej odpowiedzi kości zaproponowanąwcześniej przez autoraw celuwyprowadzenia (zamiast postu- lowania) związków rządzących przebudową kości w oparciu o bardzo ogólne założenia. Okazuje się, że takie sformułowanie ma wiele wspólnego z zagad- nieniami optymalnego sterowania. W pracy zaprezntowano przykład zasto- sowania omawianego podejścia oraz przeprowadzono krótką dyskusję na te- mat związkówmiędzy postulowanymimodelami i wyprowadzonymiw oparciu o przyjętą hipotezę. Manuscript received May 16, 2005; accepted for print June 1, 2005