## Abstract

The Eshelby stress (static energy momentum) tensor is derived for bone modelled as an inhomogeneous piezoelectric and piezomagnetic Cosserat (micropolar) medium. The divergence of this tensor is the configurational force felt by material gradients and defects in the medium. Just as in inhomogeneous elastic media, this force is identified with the thermodynamic force for phase transformations, in bone it is the thermodynamic cause of structural transformations, i.e. remodelling and growth. The thermodynamic approach shows that some terms of driving force are proportional to the stress, and some acting on material inhomogeneities are quadratic in the stress—the latter outweigh by far the former. Since inertial forces due to acceleration enter the energy–momentum tensor, it follows that the rate of loading matters and that both tension and compression stimulate growth, which is favoured at heterogeneities.

## 1. Configurational forces in bone

From time to time, bone undergoes remodelling in a cycle of osteoclastic and osteoblastic activity. Such remodelling is stimulated by physical exercise, which translates into mechanical stress on the skeleton. For a long time (Wolff 1884, 1892, 1986), it has been known that bone growth is stimulated by an applied mechanical force. A modern version of these ideas is the ‘mechanostat’ of bone remodelling (Currey 2002).

Wolff's law (1884, 1892, 1986) is the first example of mechanotransduction, the effect of mechanical solicitation on biological activity. External influence creates a wondrous wealth of shapes (D'Arcy Thompson 1942). For recent reviews on the tissue, muscular and cellular levels of mechanotransduction, see Cowin (1990, 2004), Hamil & Martinac (2001) and Kjaer (2003). Though the fact that in principle mechanical loading stimulates bone growth and remodelling has been known for 120 years, quantitative progress has been slow—the most important discovery being that not only loading but also the loading rate has an influence (Mosley & Lanyon 1998; Goodship & Cunningham 2001).

The nature of the connection between mechanical solicitation and biological response is still unclear. This is due to the complexity of the problem, which can be divided into three stages as follows.

The external loading of the skeleton, be it static or dynamic in walking, jumping, standing still, is being transmitted through a complicated loading train via muscles, cartilage, macroscopic bone architecture and mesoscopic structure (cortical or trabecular) to the osteocytes. These, in turn, control the remodelling process that is carried out by bone absorbing osteoclasts and bone depositing osteoblasts (Jee 2001). The intermediate members of the loading train with their mechanical characteristics of damping, attenuation and time delay considerably modify the external signal. The local situation of stress and strain might well be very different from the external signal applied. We are not aware of any locally controlled or monitored experiment in the sense that stress and strain were locally followed and the ensuing biological activity monitored. In principle, the problem of transmission between macroscopic external loading and local stress–strain conditions can be solved: on the theoretical side with sophisticated numerical simulations, and on the experimental side by

*in situ*measurements.Assuming that these local mechanical parameters become known, the question arises: what is the thermodynamic driving force produced by stress and/or strain or their rates? The purpose of the present paper is to provide an answer to this question. At first sight, one might surmise that growth, like bone plasticity might be stress or strain driven (Nalla

*et al*. 2003*a*,*b*) in the sense that the stress tensor or the strain tensor has a biological effect. Analogy with solid mechanics indicates that the situation is more complicated: as we will show, a certain combination of the stress and strain tensors, called the energy–momentum tensor (EMT), and more precisely, the divergence of this quantity, a vector, provide the thermodynamic driving force, which is the stimulus for subsequent biological activity.The biological response to this driving force is still another matter. It belongs to the domain of nonlinear irreversible thermodynamics. Both theoretically and experimentally, this response is the real issue of mechanotransduction: how the thermodynamic stimulus, identified by us, conditions cellular activity. As physicists, we do not address this side of the issue. Our work merely identifies the proper input parameter for the necessary biological study. The response might well be rate controlled, whatever the stimulus is. For example, the hydrodynamic part of poroelasticity might be limiting the biological response (Cowin 1999).

Another difficulty is the eminently hierarchical structure of bone (Rho *et al*. 1998). Macroscopic loadings can be surveyed, but the actual stress distribution on the osteonic level is another matter. Even the osteons, of less than millimetre size, are inhomogeneous: they consist of 30–50 lamellae of a few micrometres thickness, made up from collagen fibrils. The appropriate level for discussing osteoclastic and osteoblastic activity is not the skeletal but the microscopic and micromechanical levels.

Our objective is to treat bone lamellae within the framework of continuum mechanics. But what type of continuum mechanics? The stacking of the triple helices of the tropocollagen molecules is quite complicated (Fratzl 2003; essentially, a solidified two-dimensional liquid in the plane normal to their lengths and a one-dimensional crystal parallel to it), it is unlikely that it could be perfect; the existence of defects must be allowed for. In the continuum framework, this is taken care of by admitting the existence of a dislocation density; in other words, one allows for internal stresses in addition to the external ones from loading (Kirchner 2006). Plasticity within and between the collagen fibrils which constitute the lamellae is beyond the scope of this paper (Puxkandl *et al*. 2002; Kirchner 2006; Gupta*et al*., 2007), but the structural defects responsible for it might be relevant for remodelling. The low symmetry of the collagen assembly suggests that not only displacements but also changes of spin must be allowed for, the appropriate elasticity is micropolar (Cosserat) elasticity (Maugin 1998). For the same low symmetry reason, bone is also piezoelectric (Fukada & Yasuada 1957; Fukada 1968). Finally, bone is heterogeneous, even on the osteonic level owing to the orientation differences between lamellae, therefore, a variation of the constitutive relations with position must be allowed for. From a mechanical point of view, remodelling is the transformation of old into new material, akin to solid-state transformations between phases (Rice 1975). It is known that the influence of stress on such transformations is governed by configurational forces (Eshelby 1951; Rice 1975), which have also been called material forces (Maugin 1993; Fatemi *et al*. 2002). These, in terms, are derivable from the EMT. Our programme is, therefore, to derive explicitly the theory of the EMT for a dislocated heterogeneous elastic, piezoelectric, magnetoelectric micropolar medium thought to represent bone. For pedagogical reasons, we proceed step by step from less complicated to the most complicated medium.

The EMT is a central concept of field theory. A field theory starts with the idea that at any point in space there is a quantity, the ‘field’, scalar, vectorial or tensorial, which creates an energy density (the ‘Hamiltonian’). Minimizing the energy with respect to the field variables gives the Euler–Lagrange equations, called field or balance equations. These, usually differential, equations relate the field to its source.

By considering a slight shift of the field with respect to the material (taking a ‘variational derivative’) with the sources and the material kept fixed, one obtains new balance equations. These relate, usually in differential form, the EMT to configurational forces acting on the sources and material inhomogeneities (Wenzel 1949; Barut 1964; Maugin 1995).

Maxwell (1873) constructed the tensor in question for electromagnetic fields, Eshelby (1951) constructed it for elasticity, (Kluge 1969*a*,*b*) for Cosserat elasticity and Cherepanov (1974) for piezoelectricity. Other generalizations pertain to non-local micropolar (Lazar & Kirchner 2006), gradient (Lazar & Kirchner 2007*a*) and microcontinuum (Lazar & Anastassiadis 2006; Lazar & Maugin in press) elastic field theories.

The formal theory will result in expressions that give clues to what can cause remodelling. Our results will show that some of the configurational (thermodynamical, material) forces are proportional to the stress, some proportional to its square. The latter seem quantitatively more important than the former. It is hoped that, once medical experiments will have decided if remodelling responds to one or the other, conclusions might be drawn about the micromechanical mechanism of remodelling. For the time being, comparison with structures is restricted to simulated ones (Weinkamer *et al*. 2004).

## 2. Electrostatics

This is the paradigm of field theory, one considers the EMT in its simplest form. The energy is(2.1)where *E*_{i} is the electric field and *D*_{i} is the electric displacement. The integral is taken over the volume of the body. For simplicity, we assume a linear relationship between the electric displacement *D*_{i} and the electric field *E*_{i}(2.2)where *ϵ*_{ij} is the tensor of dielectric constants. Observe here that the anisotropy permitted by allowing *ϵ*_{ij} to be a symmetric (and not a diagonal) tensor has no bearing on the argument. The field equations are, with *q* being the electric charge,(2.3)(2.4)The latter equation implies that the electric field must be the gradient of an electric potential *V*(2.5)which, in turn, implies that there is a certain symmetry to the gradient of the electric field,(2.6)

Now let us take an arbitrary infinitesimal functional derivative *δW* of the energy density. From equations (2.1) and (2.2) we get(2.7)In the following, we want to get a configurational force, therefore, we specify the functional derivative to be translational (no summation on *k*):(2.8)On the left-hand side of equation (2.7), we write(2.9)with the energy density(2.10)On the right-hand side of equation (2.7), we have, by using the symmetry (2.6),(2.11)The first term of equation (2.11) may be written with (2.3) according to(2.12)By equating equations (2.9) and (2.11), using equations (2.10) and (2.12), we obtain(2.13)The second integral contains the sources of the electric fields: the electric charge and the inhomogeneity of the material. The integrand of the first integral in equation (2.13) is the divergence of the EMT of electrostatics(2.14)In addition, the configurational force density is defined as the divergence of the EMT(2.15)with the configurational (energetic) force(2.16)Already the simple electrostatic example shows the characteristic feature that will reappear for all other fields we will consider (elasticity, piezoelasticity and Cosserat elasticity). If the material is not isotropic, the EMT is not symmetric. The first term is the divergence of the EMT, which is the configurational force, and is proportional to the field (electrostatically *E*_{k}); the second term is of quite another nature, being quadratic in the field but proportional to the gradient of the material constants (electrostatically the dielectric constants *ϵ*_{ij}). Even in simple electrostatics the EMT is quite different from general relativity, where it is symmetric and divergence-free (there is no anisotropic constitutive law and configurational forces are absent).

The divergence of the EMT can be integrated out with Gauss, and we find(2.17)If the variation of the dielectric constant is discontinuous, and does not allow to form the gradient , the usual argument of forming a penny-shaped volume across the surface between two media I and II gives as configurational force on the interface(2.18)Here and are the EMTs on the sides I and II of the interface, respectively.

## 3. Elastostatics

The introduction of the EMT to elasticity by Eshelby (1951) is considered as one of the most important advances of that field in the twentieth century. The concept was extensively applied to fracture mechanics by Cherepanov (1967) and Rice (1968). The elastic energy is given by(3.1)where *σ*_{ij} is the stress and *β*_{ij} denotes the elastic distortion (Kröner 1958). For simplicity, we assume a linear relationship between the two, but, as in electrostatics, that is not at all necessary. Hooke's law for full anisotropy reads(3.2)where *C*_{ijkl} is the tensor of elastic constants which possesses only the symmetry(3.3)Observe that the fields are now tensorial, not vectorial any more, and the constitutive law involves a tensor of fourth order and not of second order. It should be noted that one can think of going from the electric to the elastic case by adding on index in front: and . The field equation is(3.4)and the incompatibility condition reads(3.5)The latter equation implies that the dislocation density tensor *α*_{ij}, which can be considered as the source of the elastic distortion, must be divergence free in the second index(3.6)The concept of a dislocation density *α*_{ij} does not depend on the existence of a lattice, it is merely the curl of the, in principle, measurable distortion field *β*_{ij}. When one admits that bone is describable by continuum mechanics, one is forced to allow, at least formally, for the definition (3.5). Presumably, a dislocation density is to interpret as some fault in order, for example, in the highly correlated collagen sequences. By multiplying equation (3.5) by *ϵ*_{mnj}, one finds for the elastic distortion(3.7)If no dislocations are present, the elastic distortion is just the gradient of a displacement *u*_{i}: .

Now let us take an arbitrary infinitesimal functional derivative *δW* of the elastic energy density. From equations (3.1) and (3.2), we get(3.8)As before we specify the functional derivative to be and on the left-hand side of equation (3.8), we write again(3.9)with the energy density(3.10)On the right-hand side of equation (3.8), we have(3.11)where the second and third terms have been subtracted and added. The purpose is to obtain the square bracket with the meaning of equation (3.7). The third term may be written with (3.4) according to(3.12)By equating equations (3.9) and (3.11), using equations (3.10) and (3.12), we obtain the expression(3.13)The second integral contains the sources of the elastic fields: the dislocation density, the body force and the inhomogeneity of the material. The integrand of the first integral in equation (2.13) is the divergence of the EMT of elasticity (Eshelby stress tensor)(3.14)Again, the configurational force density is the divergence of the EMT(3.15)with(3.16)The first term is the so-called Peach–Koehler configurational force on a dislocation density *α*_{ij} (Peach & Koehler 1950). The second term is the Cherepanov (1967) configurational force on a body force *f*_{i}, and the last term, the Eshelby (1951) force on the elastic inhomogeneity . In §8 we will show that in bone, due to its eminently heterogeneous microstructure, the last term is by far the most important one.

As in electrostatics, the divergence of the EMT can be integrated out with Gauss, and formula (2.17), but now with *P*_{ik} from equation (3.14) is valid. If the variation of the elastic constant is discontinuous, and does not allow to form the gradient , formula (2.18) is applicable, but with and being the EMTs (3.14) on the sides I and II of the interface, respectively.

By now, the structure of the EMT has become clear: it is the energy density in the diagonal minus the direct product of the two (vector) electric fields, or minus the (tensorial product) of the elastic fields (Kirchner 1999).

## 4. Piezoelectric medium

A combination of electric and elastic effects is the situation of piezoelectricity. Since collagen and bone are known to be piezoelectric, this is of some relevance. In a piezoelectric medium, there exist not only the electric field *E*_{i} and *D*_{i}, and the mechanical fields *β*_{ij} and *σ*_{ij} independent of each other, but also the constitutive relations couple the electric and mechanical fields. The energy density is(4.1)where *e*_{ijk} is the tensor of piezoelectric constants. The constitutive relations are(4.2)(4.3)After algebraic manipulations, we obtain the EMT(4.4)and its divergence, the configurational force density(4.5)We note that due to the electric part in the stress (4.2), the Peach–Koehler force in a piezoelectric material (first term in equation (4.5)) has an electric part. Therefore, the Peach–Koehler force is caused by both stress and electric fields on dislocations in a dislocated elastic dielectric. This observation is in agreement with Maugin & Epstein (1991) and Maugin (1993).

## 5. The field equations for an inhomogeneous Cosserat continuum and electromagnetic fields

In classical elasticity, treated in §3, every point of the body is just a point that suffers a displacement upon deformation. From the point of view of solid mechanics, because the smallest volume element of bone is to be considered as a molecular crystal, each point is to be associated also with a direction (the orientation of the collagen). Upon deformation, each point not only undergoes a displacement, but the orientation vector changes as well. The appropriate formalism is the Cosserat elasticity, in which, besides the asymmetric force stress *σ*_{ij} caused by the micropolar distortion *γ*_{ij} there exists also a couple stress *μ*_{ij} caused by a wryness tensor . Owing to the piezoelectric properties of bone, one must also allow for the presence of electric fields, as in §§2 and 4. The inclusion of magnetic fields of possible therapeutic significance is straightforward.

The equilibrium equations for the elastic stresses (force stress *σ*_{ij}, couple stress *μ*_{ij}) are(5.1)(5.2)and the static Maxwell equations read (the electric displacement *D*_{j}, the electric field *E*_{j}, the magnetic field *H*_{j} and the magnetic induction *B*_{j})(5.3)(5.4)(5.5)(5.6)where *f _{i}*,

*l*,

_{i}*q*and

*j*

_{k}are the body force, body couple, electric charge density and current vector, respectively. For the current vector, we have the conservation law(5.7)It can be seen in equation (5.2) that, in addition to body couple, the skew-symmetric part of the force stress

*σ*

_{ij}gives a source term for the couple stress

*μ*

_{ij}in a Cosserat or micropolar elasticity. On the other hand, the force stress

*σ*

_{ij}has only the body force as source in equation (5.1). The additional source term in equation (5.2) will be important for the configurational forces in Cosserat elasticity.

In linear micropolar elasticity, the incompatibility equations are(5.8)(5.9)By differentiating, we obtain the conservation laws for the dislocation and disclination densities(5.10)(5.11)which mean that the disclination density tensor is divergence free in the second index and the divergence of the dislocation density tensor is determined by the skew-symmetric part of the disclination density tensor. If no dislocations and disclinations are present, the micropolar distortion is given in terms of the gradient of a displacement *u*_{i} and a micro-rotation vector *φ*_{i} and the wryness is just the gradient of the micro-rotation(5.12)(5.13)In addition, the electromagnetic fields have the potentials(5.14)(5.15)where *V* denotes the electrostatic potential and *A*_{l} is the magnetic vector potential.

In the case of material linearity, the energy density is given by(5.16)in terms of the fields: the elastic micropolar distortion *γ*_{ij}, the elastic wryness , the electric field *E*_{j} and the magnetic field *H*_{j}. Here the micropolar elastic tensors *C*_{ijkl}, *A*_{ijkl} and *B*_{ijkl}, the dielectric permittivity *ϵ*_{ij}, the magnetic permeability *ν*_{ij}, the tensor of magnetoelectric moduli , the tensor of piezoelectric moduli , the tensor of piezomagnetic moduli , the tensor of the piezowryness effect and the tensor of the magnetowryness effect depend on position. Only, the symmetries hold(5.17)The energy within a fixed volume of the body, over which the integrals will be taken is(5.18)Consequently, the constitutive equations where the material tensors vary with position read (see also Eringen (1999) for homogeneous media)(5.19)(5.20)(5.21)(5.22)

## 6. The configurational force for bone growth and remodelling

Now we construct the static EMT of a Cosserat medium, characterized by the field equations (5.1)–(5.9), (5.17), (5.18) and the constitutive relations (5.19)–(5.22). As in §§2–4 for electrostatic, elastic and piezoelectric media separately, both sides of the energy expression (5.18) with (5.16) are subjected to the variation . The result is (see also Lazar & Kirchner (2007*b*))(6.1)where and denotes the plastic part of the micropolar distortion tensor. It is seen that the EMT for the quite complicated Cosserat medium is a mere addition of the EMTs for an electrostatic (2.14) and elastic (3.14) medium, plus the magnetic and Cosserat terms. The complexity of the medium does not appear in the structure of the EMT, but in the coupling between the fields in the energy density (5.16), and the constitutive relations (5.19)–(5.22).

Now, more important, the divergence of it gives the configurational (material), thermodynamic forces on the sources and inhomogeneities. The result is(6.2)with(6.3)Equation (6.3) is a sum of configurational force densities: the Peach–Koehler force density on a dislocation density *α*_{il} in the presence of the force stress *σ*_{ik} (Peach & Koehler 1950), the Peach–Koehler force density on a disclination density in the presence of the couple stress *μ*_{ik} which is the Mathisson–Papapetrou type force density also known from gauge theory (Gairola 1981; Maugin 1993; Hehl *et al*. 1995), the Cherepanov force density (Cherepanov 1981) on a body force *f*_{i} in the presence of the micropolar distortion , the force density on a body couple *l*_{i} in the presence of the elastic wryness , the force density on the force stress *σ*_{ki} in the presence of the plastic wryness which is due to the fact that the skew-symmetric part of the force stress is a source term of the couple stress (see equation (5.2)), the electrostatic force density on the electric charge *q*, the Lorentz force density on a current *j*_{l} in the presence of a magnetic field *B*_{k} and 10 force densities on inhomogeneities like the force density on the elastic inhomogeneity derived by Eshelby (1951). Let us mention that the Peach–Koehler forces due to dislocations and disclinations and also the third part of equation (6.3) contain electric and magnetic parts in addition to mechanical parts. Thus, the Peach–Koehler forces are caused by stresses, electric and magnetic fields on dislocations and disclinations in a piezomagnetic, piezomagnetic and magnetoelectric Cosserat medium. These electromagnetic terms could drive the healing process and the remodelling of bones.

If there are no sources, *α*_{ij}=0, *Θ*_{ij}=0, *f*_{i}=0, *l*_{i}=0, *q*=0, *j*_{l}=0 and no plastic fields, , , the configurational force is felt only by the gradients of the elastic, micropolar, piezoelectric, piezomagnetic and magnetoelectric constants. In this case and without inhomogeneities, the form of the EMT in equation (6.1) is divergence free. Without electromagnetic fields, equation (6.1) agrees with the expression for the so-called ‘Maxwell's stress tensor of a Cosserat continuum’ given by Kluge (1969*a*). For compatible micropolar elasticity such an expression was derived by Lubarda & Markenscoff (2003). In addition, Pucci & Saccomandi (1990) deduced a similar expression for defect-free isotropic micropolar elasticity. Nikitin & Zubov (1998) derived such tensor for finite deformations of compatible micropolar elasticity (see also Maugin (1998)).

If interfaces (where these constants are discontinuous) are present, the gradient terms in (6.3) become singular. Formula (2.18) is valid, but now with and as the Eshelby stress tensors (6.1) of the phases I and II, respectively. This formula is the extension of Rice (1975) who derived it for phase transformations in materials science. There the force on the interface causes its movement during the phase transformation.

## 7. Application to bone remodelling

Both the mechanical stress *σ*_{ij} and the EMT *P*_{ij} have the dimension of an energy density (J m^{−3}) or force per unit area (N m^{−2}), which are formally the same, because J=Nm, but this observation is misleading. Since the enunciation of Wolff's law (1884) biologists, biophysicists and the biomechanics community were tempted to argue in terms of stress, and not the energy–momentum tensor. It is therefore important to clarify the conceptual difference between the two.

The divergence of the mechanical stress is related to an applied body force by the field equation (3.4), while according to (3.15) with (3.16) the divergence of the EMT gives the configurational force density on sources and inhomogeneities present. Both equations have the same structure and dimension, but their meaning is different. In the field equation (3.4), the source of the stress field, the body force density *f*_{i} is prescribed (imposed) and the differential equation must be solved to obtain the stress field *σ*_{ij} everywhere. Once this has been done, the EMT *P*_{ij} is constructed according to (3.14), and its divergence, according to (3.15) and (3.16), gives the configurational force felt by the force *f*_{i} due to the stress field present. More generally, the line of argument goes from sources of the field (body forces, dislocations, disclinations, charges and currents) to the fields, which are solutions of the field equations. From these solutions, the EMT is constructed everywhere, and its divergence, in turn determines the configurational forces on field sources and material inhomogeneities. Fracture mechanics is a simple illustration of this remark: the loadings are prescribed, the stress fields computed and the crack extension force is obtained via the two-dimensional version of the EMT (Cherepanov 1967; Rice 1968).

Under a given loading of a body, say by a point force or surface loading, in linear elasticity the stress is proportional to the loading, and, due to Hooke's law, the energy density (and the other terms in the EMT as well) are quadratic in the loading.

Inspection of the resulting configurational forces (6.3) is very revealing: the first terms, the configurational forces on the field sources, in bone that would be dislocations *α*_{ij}, disclinations and electric charges *q*, are proportional on any external loading or applied electric and magnetic field. On the other hand, the last 10 terms, the ‘material’ forces acting on inhomogeneities (where the elastic, electric, piezoelectric and Cosserat constants vary with position) are quadratic in the external loading.

At present, a concrete mechanism how dislocations and disclination could act in remodelling has not yet been envisaged, but if they play a role, their action should be proportional to the loading. The other way round, if physiologically stress is considered as important, that indicates that defects are essential for remodelling. More likely seems the hypothesis that deterioration of bone quality causes a material change with ensuring changes in the constitutive laws, thereby producing inhomogeneities. If remodelling is connected with these, their action should be proportional to the square of the stress (the energy). One concludes that, if defects matter, the thermodynamic driving force is linear in the applied loads (or electric and magnetic fields); if inhomogeneities matter, the thermodynamic driving force is quadratic in the applied loads (or electric and magnetic fields).

In the context of bone, the configurational (material) forces on interfaces between two microstructural entities (osteons and cement lines, etc.) cause an interfacial shift. The analogon of a phase transformation in materials is growth (or remodelling) in bone. Since the material constants in bone are not known, even not order of magnitude wise, it is difficult to judge the magnitude of the forces driving growth. It is, however, clear that the configurational forces acting on interfaces are quadratic in the fields present.

## 8. A quantitative estimate: inhomogeneity and rate effects

As simple, but instructive example we compare the magnitude of the Eshelby and the Cherepanov terms in (3.16) for an approximately modelled femur (see figure 1). Since the piezoelectric and Cosserat coupling constants have never been measured, we restrict the argument to classical elasticity of §3. Because real strain distribution is very complicated indeed, we assume that the diameter of the femur is small compared with its length. All three- and two-dimensional effects are neglected, only the normal stress *σ*_{33} and distortion *β*_{33} in the vertical direction need to be considered and we can suppress the indices. In the inferior epiphysis of height *H*, Young's modulus is supposed to vary sinusoidally with an amplitude *A* and a wavelength *λ* like . The maximum gradient of the Young modulus is .

In the diaphysis, there are volume force densities *f* due to gravity () and dynamical ones, . In typical movement, the acceleration *a* is approximately 10 times greater than *g*. The load transmitted from the diaphysis of length *L* is transmitted by a traction , where *L* must be taken longer than the actual femur to account for the weight of the anatomy above. Within the metaphysis, the distortion *β*(*t*_{a}) due to the traction *t*_{a} is of the order of *t*_{a}/*E* or . The distortion *β*(*f*) due to the body force density *f* is of the order of *σ*(*f*)/*E*, if *σ*(*f*) is the stress. On the average, the stress balances the total body force density, *σ*(*f*)=*fH*, which gives, order of magnitude wise .

According to equation (3.16), the Cherepanov force is and the Eshelby force is . The ratio of the Eshelby configurational and the Cherepanov configurational forces is therefore: . If one takes, reasonably enough, the amplitude to be , the Eshelby configurational force is a factor larger than the Cherepanov one. Certainly *L*>*H*, and can be taken of the order of 1 m to account for the rest of the skeleton. The inhomogeneity scale *λ* is the lamellar thickness within an osteon, 3 μm (Gupta *et al*. 2006), or the thickness of the cement lines, which is about the same. One reaches the conclusion that the Eshelby force on inhomogeneities is a factor *L*/*λ*, or 300 000 times larger than the Cherepanov force.

The absolute value of the thermodynamic driving force acting at material inhomogeneities is, therefore, of the order of . The stress *σ*_{a} is approximately 100 MPa, Young's modulus approximately 10 GPa and, with *λ*=3 μm one obtains a value of , independent of the sign of the stress present. This is much larger than the force density due to gravity, which is for a bone density *ρ*=1, or 10 times as much during physical exercise.1

The expression (3.16) for the EMT, derived for static fields, remains valid for the dynamic situation where inertia effects exist. It suffices to add the inertial force to the body force as it occurs in the Cherepanov term, as the ‘static’ form of the EMT remains valid also for inertial forces. Since inertia forces are proportional to acceleration, their preponderance in the Cherepanov part of the configurational force imply, and allow for, the existence of rate effects as observed by Mosley & Lanyon (1998). They demonstrated that not only the time average of the load matters but also the actual functional time dependence. Contrary to superficial appearance, our final result for growth and remodelling rates, equation (9.2), contains such rate effects; they enter as inertial forces proportional to acceleration, in the thermodynamic driving force, the square of these accelerations appears and is then time averaged through the biological response function. Not the load, but the stimuli are time averaged.

## 9. Discussion

Analytic computation of the fields, like stress *σ*_{ij}, distortion *β*_{ij}, electric field *E*_{j}, electric displacement *D*_{j} and magnetic field *B*_{j} as function of the sources, like body loading *f*_{j}, body moment *l*_{j} and electric charge is practically impossible for any geometry of reasonable complexity and interest. Numerical methods must be used to compute Green's functions that relate the fields to the sources. A few conclusions can be drawn without explicit calculation. In spite of the same dimensions, mechanical stress (Nm^{−2}), energy density and EMT (Nm m^{−3}), behave qualitatively different. In linear elasticity, the realm of the present investigation, stress is linear in the applied load, while the energy density and the EMT are quadratic in the load.

This remark pertains to the growth model transcribed from plasticity theory (Epstein & Maugin 2000), where a ‘transplant’ tensor ^{−1}(Cowin *et al*. 1992; Fung *et al*. 1993; Rodriguez *et al*. 1994; Imatamin & Maugin 2001) describes the change from (1+*β*) to (1+*β*). An evolution law was formally suggested.

In simpler terms, growth can be described by a spontaneously created distortion rate . According to our results, it is not , but its divergence, the thermodynamic force, that controls its evolution. The evolution law for must be of div =_{thermo} over time. In linearized form that becomes(9.1)The biological response function derives from irreversible thermodynamics. If explicit age effects are neglected, so that bone development does not vary with the age of the subject, the response becomes a convolution of time:(9.2)Neither the time scale *τ* is known, nor the functional dependence of *G*_{bio} on ; it might be exponential. Biological response is slow compared with rate effects in loading. This biological response is the problem of mechanotransduction on the cellular level and beyond the scope of this paper. It might, however, be surmised that for small stimuli the response between stimulus and effect should be linear. This makes a square root biological response to the stimulus, which is quadratic in the load, unlikely and rules out linearity between applied force and growth. The essential point is that already the thermodynamic stimulus is quadratic in the load. Equation (9.2) resolves an apparent contradiction. Presumably, the response function *G*_{bio} is slow compared with the time scale of loading rates, amounting to a time average over _{thermo}, which itself is proportional to the square of the loading accelerations. It follows that contains rate effects, in agreement with experimental findings (Mosley & Lanyon 1998).

## 10. Conclusion

In the expressions for the thermodynamical driving force that triggers growth and remodelling, there are terms linear, and terms quadratic in the applied load. The latter, which dominate the former, act at material and structural inhomogeneities, like interfaces between lamellae or cement lines. The biological response to these is unknown, but it can be surmised that the formation of old with new bone starts at these inhomogeneities. Nucleation is heterogeneous and not randomly homogeneous anywhere in the material.

The analogy with phase transformations in solids is valid insofar as one structure is replaced by another one. In solids, this happens by chemical reaction, diffusionless (martensitic) or diffusive transformations; the formation of a new phase might proceed via a sharp well-defined interface, or by slow interpenetration and the formation of diffuse interfacial regions (like in spinodal decomposition). In solid-state physics, a wide variety of transport processes responds to a thermodynamic driving force. Similarly, in bone we consider the thermodynamic driving force as cause of growth (via the detour of cartilage and subsequent ossification) and remodelling (the exchange of old against new bone). Consequently, they need not involve a phase transformation front in the form of a sharp interface, but proceed by replacement.

One day, physiological test might be able to distinguish if bone growth is linear or quadratic (as we predict) in the loading. As of now, no such experimental results are available. At present, one can only compare with mechanical simulations, as those of Huiskes *et al*. (2000) and Weinkamer *et al*. (2004). These model growth to be proportional to the strain energy density and strain, with specified response *G*_{bio}—as in (9.2). However, deconvolution of the causal and the response parts is notoriously difficult in such simulations, where both are varied as long as reasonable and agreeable similarity between calculated and natural structure is obtained. One might assume the pessimistic attitude, that from such simulations the validity of our hypothesis can be neither confirmed nor refuted convincingly.

In conclusion, in accordance with Wolff's law (1884, 1862, 1986), there exists a well-defined thermodynamic driving force, a vector, which we believe to cause biological activity. We derived explicitly the form of this vector in terms of stress and strain (and higher-order mechanical quantities for more sophisticated constitutive laws). It is the divergence of the EMT familiar to physicists.

This is a configurational mechanical force vector acting on bone defects and inhomogeneities. It is not the mechanical stress tensor itself, nor the divergence of the stress (which is the body force). Clearly, therefore, the body force density is not directly causing remodelling, but acts through the divergence of the EMT (a quadratic form in stress or strain), a concept familiar from solid mechanics (Eshelby 1951; Rice 1968) and physical field theories in general (Morse & Feshbach 1953; Barut 1964).

Quantitative comparisons strongly suggest that the driving force is greatest at bone inhomogeneities. In accordance with experimental evidence, not only loading but also the rate of loading has an influence on the driving force. This is simply due to the fact that the driving force is proportional to the sum of static loading plus inertial loading, which, for example in walking, is approximately five times larger. Surprisingly, this simple explanation that inertial loading implies rate effects has never been mentioned in the literature. The quadratic dependence means that tension and compression, and loading and unloading produce the same stimulus. Since the thermodynamic driving force is a mechanical quantity, it is in principle locally measurable and controllable. Our hope is, that, with the characteristic features of the thermodynamic driving force once identified, more specific experiments on the biological response can be envisaged. What needs to be measured is not the externally applied force, but the local stresses and strains *in situ*.

We summarize our three main results as follows.

Growth and remodelling in Wolff's sense respond not directly to the mechanical loading, but to the thermodynamic driving force (the divergence of the EMT).

The preponderant term in this thermodynamic driving force, quadratic in the loading, stems from acceleration (initial) and not static loading. This observation explains the positive rate effects actually observed (Mosley & Lanyon 1998). It is the time average of the thermodynamic driving force and not the time average of the actual loading that matters.

This thermodynamic driving force is likely to be greatest at structural inhomogeneities in bone, where the material constants change. Thus, biological activity is expected to be non-uniform.

## Acknowledgments

M.L. has been supported by an Emmy–Noether grant of the Deutsche Forschungsgemeinschaft (grant no. La1974/1-2).

## Footnotes

↵Instructive is the example of a bilayer plate, with an Eshelby force on the interface. The strains are

*β*(*i*)=*σ*_{a}/*E*(*i*),*i*=1, 2, the only nonzero component of the EMT is*P*_{33}(*i*)=(*σ*_{a})^{2}/*E*(*i*) and its discontinuity across the the interface is [*P*_{33}]=(*σ*_{a})^{2}[1/*E*(1)−1/*E*(2)]=−2(σ_{a})^{2}/*E*(1) if*E*(2)=*E*(2)/2. Compared with the applied stress*σ*_{a}, the Eshelby force on the interface is 2*σ*_{a}/*E*smaller, but quadratic in the stress. Independent of the sign of the stress, the interface is always driven in the direction of the softer medium, in order to minimize the total stored energy.- Received May 7, 2007.
- Accepted July 6, 2007.

- © 2007 The Royal Society