Electro-mechanical remodelling of bones

0020-7!?5/83 53.0@+X4 % 1983PergamanPren Ltd

ELECTRO-MECHANICAL REMODELLING OF BONES HILMIDEMIRAY Marmara Research

Institute,

Division

of Applied

Mathematics,

Gebze,

Kocaeli,

Turkey

Abstract-This paper is mainly concerned with an electro-mechanical model for living bones with the hope of explaining the remodelling process in bone, which is vital for bone healing and regeneration. In d&eloping the model, it is assumed that a wet bone is composed of a charged fluid (electrolytes and bone salts) and the solid bone matrix which has the piezoelectric property. Based on this simplifying assumptions, the balance laws and constitutive relations are derived for such a mixture, and for the illustration of present model. an example problem concerning the longitudinal wave propagation in wet bones is studied and possible application of it in bone healing is indicated.

1. INTRODUCTION

THE WET bones, like many other tissues (i.e. tendons) have the capacity of regeneration and remodelling [ I], which is known as the Wolf’s law and controlled by a feedback mechanism. It is generally believed that the main cause of this feedback mechanism is the interaction of electro-mechanical and chemical fields in the bony tissues[2]. The internal forces caused by external effects generate electrical signals which carry information to bone cells to regulate their biological functions and controlling the orientation of their secreted macromolecules [3]. From the physical point of view. the structure of a wet bone tissue can be considered as the collection of a fluid phase (blood) and the solid matrix. In the process of remodelling, the bone salts in the blood are transferred to bone matrix and vice versa. On the other hand, due to leakage of monovalent ions (Na, Cl, CO,) into the charged parts of the bone matrix (polarization charges), the internal electrical potentials decay very rapidly[l]. Observing these two distinct phenomena taking place in living bones, one may be tempted to consider two fluid species in the blood which are blood salts (Ca,, PO&) and the monovalent ions (Na, Cl, HO,). Such a consideration makes it difficult to build up a tangible mechanical model for bones. Nevertheless, treating the whole blood as a single fluid with independent mass and charge densities might still give a satisfactory and yet mathematically much simpler result. Since the first discovery by YasudaI41. that the bone is a piezoelectric material, there appeared some models in literature for dry and and wet bones. Due to simplicity, the dry bones have been studied somewhat in detail by Gjelsvik[S], Williams and Braeger[6] and Giizelsu[7]. On the contrary, because of the mathematical and physical complexity of the problem, the progress made on wet bones is rather meagre. Among such studies it is worthy of mentioning the works by Petrov[8] who utilized the model presented for solid state plasma[9] and Demiray and Giizelsu[lO] who suggested a three component model for bones. Although the motivation in[8] is encouraging, due to misuse of the axioms of continuum mechanics, the constitutive relations rather seem to be heuristic. The three component model of [lo] is complete though, mathematical complexity of the equations makes it difficult to apply the result to any problem of practical interest. In this work we will present an electromechanical model for wet bones consisting of a piezoelectric solid bone matrix and a fluid (blood) with separate mass and charge densities. It is assumed that the electrochemical reactions take place between the solid and fluid continua without interacting with charge densities. The latter assumption holds true when the blood is supplied with sufficient charge through nourishment. For an electro-mechanical model for bones with these characteristics, the kinematics and balance Iaws are briefly summarized in Section 2. Section 3 is devoted to study of thermodynamics and nonlinear constitutive relations of bones, including the charge and conduction current densities (as well as usual variables) as independent constitutive functions. In the last section the linear constitutive relations and the solution of an example problem have been presented and the applicability of present model is thoroughly discussed. 1117

We consider a mixture consisting of a solid bone matrix and a charged fluid which occupies the Euclidean three-space V. A material point N,,,, (TV= I. 2) of the cuth species of the mixture is carried into a spatial position x, at time t, through its proper motions described by

Throughout this work. Greek indices enclosed within parenthesis label the constituents of the mixture. The mapping is assumed to be one-to-one and onto, thus the inverse motion exists and characterized by X(,,) = X(,,,(X. I).

(a = !,?I.

(2.2)

Existence of inverse motion implies that (2.1) and (2.2) possess continuous partial derivatives with respect to its arguments, except, possibly at some singular points, lines or surfaces. Thus the Jacobian of motions must not vanish, i.e. (2.3) The velocity and acceleration vectors of species are defined by (2.4) with

Here the summation convention applies on repeated indices. For its use in the subsequent analysis, it might be useful to report the following identity. (2.5)

Having this much about the kinematics of mixtures we may list the balance equations governing the species and the mixture as a whole. Following Demiray[9, 1t] the balance equations of charged mixture are listed below: (i) Conservation of mass

where p(a) is the mass density and c(,) is the rate of mass transfer from other species into ath component. Since there is a mass transfer among species, the density of solid continuum cannot be obtained through the relation pt,, = Jopt,,, where pFY,is the initial mass density of the solid bone matrix. Throughout this work the letters s and f enclosed within the parenthesis are used to label the solid and fluid phases, respectively.

Electra-mechanical

where t$’ is the the rate of linear semi-colon are Throughout this electromagnetic

1119

remodelling of bones

stress tensor, gj”’ is the electromagnetic body force per unit volume and RI”’ is momentum transfer from other species into ath species. The indices following used to denote the partial derivatives with respect to that coordinate. work the effect of mechanical body forces will be neglected. The expression of body forces is given by gi“ = EM&,

gy’ = qE,

(2.8)

where E and P are respectively the electric field density and polarization density of solid matrix, and q is the charge density of fluid continuum. Assuming that the process is slow, we have neglected the magnetic effects. (iii) Balance of angular

momentum

t& + ttfkj,+ P,kE,, = 0. Here the indices enclosed associated tensor, e.g.

in brackets

(2.9)

are used to denote

t;;/] = &l;;‘-

the skew-symmetric

$‘).

part of the

(2.10)

of energy

(iv) Conservation

p,,,i,,,

= @VI;)+

qti + h,,, - Rla)vla)+ P)(&&,

- E& + e(,,

(2.11)

with

T e(,, = 0.

(2.12)

where Ed,) is the internal energy density, qip’ is the heat flux, I+,, is the electromagnetic energy density per unit volume (thermal heat sources are neglected) and e,,, is the rate of energy transfer from other components into ath species. For its brevity, the over-dot is used to represent the material time derivative. The expression of electromagnetic energy is given by h(,, = E&, (v) Electromagnetic

h,, = 0,

+k = P, + Pkv;!,.

(2.13)

field equations Vx E = 0 0 .D = q

$

+J =0

2 + div(qv”‘) = 0 where D is the electric displacement

(Faraday’s

Law)

(2.14a)

(Gauss (Law)

(2.14b)

(Amp&C’s Law)

(2.14~)

(conservation

(2.14d)

of charge)

vector and J is the total electric current density defined by D = E,,E + P,

J = qo”.

Here e0 is the dielectric constant of vacuum. Again, assuming slow, the magnetic effects have been neglected.

(2.15) that the dynamical

process is

1110

(vi) Principle of entropy Although there are several views as to whether the entropy inequality should he formulated for each species in the mixture or a single inequality for the whole mixture. in this particular work we shall use a single entropy inequality as follows:

where q(o) is the entropy density of ath continuum and H is the common absolute temperature of the species. For multi-temperature distribution it might be pertinent to use separate entropy inequality for each member. Introducing the following Legendre transformation

and using the energy eqns (2.11) and (2.12) we have - &)(&)f

no,&-&&,,t

T)&+

t:‘cI.;‘+ t$‘ry:

- P&I

+(qV’+

4:“)

61,

7-

c(&,,-

ICltrJ

+ T,j, ~0.

(2.18)

where we have defined T,, j, and c

Here j, stands for the conduction current of fluid continuum relative to solid matrix. We require that the inequality (2.18) should be valid for all thermodynamiccrlly admissible processes.

3. CONSTITUTIVE

EQUATIONS

In this section we shall develop a set of constitutive equations of a mixture consisting of a solid bone matrix and a non-viscous fluid (blood).t For this purpose, the constitutive dependent variables will be selected to be

and the independent

constitutive

Here one should note that, since pI) cannot only be obtained from the deformation and the initial mass density of solid continuum, we have chosen it also as an independent constitutive variable. The general form of a constitutive dependent function should read

The constitutive functions are further restricted by the axiom of objectivity, material symmetry and the second law of thermodynamics. The first implications of the second law of thermodynamics may be stated as &,) =

‘i’t,#iK. Ptst, 0. Ek). (3.4)

&I, = ~(~)(~(~~. 4. 0). +Although it is well known shall neglect its effects.

that blood possess a considerable

viscosity.

for the sake of simplicity

in the derivations,

we

1121

Electra-mechanicalremodellingof bones Taking the material time derivatives equality (2.18) we obtain

of (3.4) and introducing

the result into the entropy

in

Since the inequality (3.5) is linear in the variables v w ,.k, &, i and &k, in order that (3.5) be valid k4 &, i and &k, the coefficients of these quantities must vanish, for all arbitrary variations of U,J, e.g.

a

!!b FkK, t:“I’= - p:s)apes, ak, + p(s) IK t$‘=-&jkI

rl(a)=-

The remaining

aF

Pk

=

-

p(s)

2

!!hfl Pcr'4+f

2 apv, 7T-= PW

(3.6)

+

(3.7)

alCr(ol, qt’ + qr’ = 0 (no net heat-flux). 3

(3.8)

a

a0

parts of the inequality

(3.5) take the following form (3.9)

It is interesting to note that for chemically reacting charged mixtures the solid component constitutionally has a pressure term, while the fluid part has a pressure contribution from the variation of charges of the fluid continuum. As a matter of fact, the latter observation implies that the gradient of charge density creates a driving force on the fluid particles. This point will be quite clear as we study the linear constitutive relations. The constitutive relations are further restricted by the axiom of objectivity which states that the constitutive functions should remain form invariant under the full group of rotation of the spatial frame of reference. The implication of this axiom may simply be stated as that the constitutive dependent variables should be a function of 0, pea,,q and the following set EKL

=

;(FkK&-

8KL), EK =

Thus, from (3.5), the stress tensor and the polarization t 8)= - &)

a

JPW

Sk,+

(3.10)

vector of solid bone matrix follow

!i!hFkKE,

FkK&_ + P(s) ‘(” ?bd aEKL

(3.11)

$& FkK.

(3.12)

pk = - pcs,

Equations which are solve any equations

EkFkK, IK = FkKjk.

aEK

(3.7), (3.8), (3.11) and (3.12) give the most general forms of constitutive relations anisotropic and nonlinear. However, with this generality, it is not an easy task to practical problem. Therefore, in what follows we will study the linear constitutive and hopefully see some of its applications

II”

H.DEMIRAY

_I

4.LINEAR CONSTITUTIVE RELATIONS In the linear constitutive theory the difference between the material and space coordinates may be disregarded. This is tantamount to stating that the displacement field is small. We further assume that the displacement gradients, electric field and the deviation of densities from their initial values are so small that the second order effects may be neglected. Thus, for the linear theory the strain energy of solid bone matrix may be expressed as

where i)(S)stands for the deviation of solid mass density from the initial density, and the coefficients have the following symmetry properties

Introducing (4.1) into (3.11) and (3.12) and noting the following relation FkK

=

SkK

(4.3)

-+- SkMhf.K

where && is the coordinate shifter and uK is the component of the in~nitesimal displacement vector; the stress and polarization vector of solid bone matrix read,

Here in obtaining these equations we have kept only the first order terms in the expansion and assumed that the stress and polarization vanish as Em,. EK and bcy,go to zero. The stress tensor of fluid continuum, in the linear theory, from (3.7), may be expressed as (4.6) where x0 is the hydrostatic pressure, p&, is the initial mass density, q. is the initial charge density and brn and 4 denote the deviations of mass and charge densities from their initial values. Now let us return to inequality (3.9). For the linear constitutive relations the proper form of T, and c should be as follows Tl= c = d,;,,,td&,+d,4

vj,

+ GE,

(4.7) + CKLEKL.

(4.8)

Here Y is of order of G/q:, where V is the collision frequency of sofid and fluid particles. Furthermore, a priori assumption has been made on the vanishing initial value of mass transfer term. Besides v z 0, substitution of (4.7) and (4.8) does not yield any simplification on the material coefficients di, CK and CKL. The constitutive equations are further restricted by the axiom of material symmetry conditions. The bone matrix is known to possess hexagonal polar (C,>structure (Anderson and Ericksson [ 121, Fukada and Yasuda[4]). Selecting the plane perpendicular to X, axis as the plane of rotational symmetry, the matrices of material moduli take the following forms

Electra-mechanicalremodellingof bones

1123

- B,d

(4.9)

Formal expressions may be given for other tensors associated with material coefficients. In the linear theory we disregard the differences between the material and spatial coordinates. Hence, the linearized constitutive equations take the following final form t:‘r’= - pl),&Z& - &E, + Bmnemn)&, + B&W - GmEm + Bk,mnemn

(4.10) (4.11)

(4.12)

Ri” = qovjr c =

d,fi,,,

(4.13) +

where e,, is the infinitesimal

d&+

d34

+

CkEk

+

(4.14)

Cklekl

strain tensor defined by 2ekl

uk.1

=

+

(4.15)

U1.k

liL being the components of the displacement vector u. Furthermore, we have approximated terms of T, as given by eqn (4.13). The linearized form of the balance equations follows from (2.6)-(2.8) and (2.14).

R, in

(4.16)

t;{ik+ q,&, + GE: -

(5) _

0

auY’

(4.17)

R , - Pu, 7;

(4.18) where cp(x, t) is the quasi-static electric potential, ik is the external current density and 4 is the external free charge density: Substituting (4.10)-(4.15) into (4.16)-(4.18) a set of differential equations governing the fields &, vif’, pin), 4, jk and cp(totally 13 unknown) can be obtained. The number of differential equations is sufficient to determine the electromechanical field completely. Appropriate boundary and initial conditions may be posed for each specific problem of concern. 5.ONEDIMENSIONALPROBLEMSIN

BONES

In general it is very difficult to deal with 3-dimensional problems in bones, since it is anisotropic in its structure and contains other species interacting with it. It will be instructive to study first the l-dimensional problems. For instance, the waves propagating along a broken bone had been used by pediatricians to speed up the healing processes. Assuming that all the phenomena occurs along the bone axis (x3 = z axis), in l-dimensional case we have

(4.19)

II?4

H.DEMIRAY

and the remainings are all zero. Introducing (4.19) into (4.10)-(4.14) we obtain

where the prime denotes the differentiation with respect to space variable z. Introducing (4.20) into (4.16) and (4.17) the following patial differential equations ire obtained 2

(4.21)

(4.22)

(4.23)

(4.24) (4.25) (4.26)

(4.27) Here we have assumed that the external charge and the current densities vanish. Before proceeding further some interesting results can be deduced from eqn (4.22). Assuming that the inertia of the fluid and the mechanical pressure are small and they can be neglected, we have (4.28) Here u = l/v stands for the conductivity of the fluid medium. It is easily seen that the gradient of induced charge density contributes the induced conduction current. Further, assuming that density gradient of solid bone matrix is small, from (4.21), (4.25) and (4.27), one obtains

- (cotA&"

t C3+” - cj= 0

(4.30) (4.31)

Here we have made use of the eqn (4.28). This gives us three coupled equations to be solved. Knowing cp,u and $, from the remaining eqns (4.22)-(4.24) and (4.26) fiCaJ and u can be determined. Solving 4 from (4.30) and substituting it into (4.29) and (4.31) we obtain

Electra-mechanical

remodelling

1125

of bones

Here f(t) is an arbitrary function of time-t-which may be taken to be zero at the outset. Now let us seek a harmonic wave type of solution to eqns (4.32) and (4.33). For this, we set (4.34) where ((I+,, U,) are the complex amplitudes, k is the wave number and w is the angular frequency of the wave. Introducing (4.34) into eqns (4.32) and (4.33) we obtain the following algebraic equations [-

ik’a2(Eo+ Ad + (P:,& - G3)kZ-

~~qd~o+ [P&o2+ i(u$J,,,k3 + (pI)&3 - &dk21 U, = 0 (4.35)

- [ik3Cu2(e,,+ A??) + kwq,v(e,, + A3J + ikq&,+

[ik’azC?33 + kwqovCq?x- ioq,‘,vl fJ,, = 0. (4.36)

In order to have a non-zero solution for a0 and U,,, the determinant from (4.35) and (4.36) must vanish

of coefficients matrix obtained

q,dvt’,,w’+ ipl’,,(a+ a&k2b2+ qd(ab + ~:+43C333 - C&3 - w2b)k2 - id’+@ -

dlo

+ i[dab + ~&A3C333 - C$,)k4+ qoak21= 0

(4.37)

B 3333,

(4.38)

where we have defined a =

,4$33

-

b =

EO+

Au.

Before studying the general dispersion equation (6.37) it might be illuminating to examine some special cases. (i) High conductivity. In this limiting case u is vanishingly small, so that we neglect terms in (4.37) including v. Thus, d&o

+ Gk’b2

+

[dab

+ ~&43C333

-

Cf3,)k4

f

wJd = 0

or o2 = WCC&~

- ~&A~&33 pcl,,(qo

-

abW4

-

wok2

+ &k’>

Naturally, the resulting wave is dispersive. Moreover, in this limiting case there is only one wave propagating in the medium. (ii) Low conducticity. Poor conductivity implies very large values of v. The terms in eqn (4.37) with v therefore, may be considered to be large as compared to others. Thus, we have

w7 + T&

PCs+

[tab

+ P&A~C&

-

C:33 - q,a,b)k’-

iqop&A3k - q$o = 0.

(4.40)

This again gives us a single and dispersive wave. Further, we note that the existence of initial charge density leads to dissipation of the wave, i.e. the frequency has an imaginary part. On the other hand, vanishing initial charge density of the fluid continuum gives non-dispersive wave for both limiting cases (see eqns (4.39 and (4.40)). (iii) Vanishing initial charge density. If the initial configuration of the fluid continuum is charge free (q,, = 0), then from (4.37) one obtains w2 + 4 (ab + p&W333 pc.4

- C:,,)k’ = 0.

(4.41)

II26

H. DEMTRAY

This is again anon-dispersive wave solution. It is also interesting to note that the phase velocit) . f~ this particular case, does not depend on the coefficient (Y?which represents the variation of fluid pressure with electrical charges. Some other intermediate cases may be studied as well. For instance. for moderately high conductive cases (small v), one can expand the frequency as a series of v and obtain further approximations for w. It is being a straightforward manipuIation~ the result wiil not be reported here. Since, there is no complete experiment giving the numerical values of material coefficients, it becomes rather difficuft to study further the general dispersion equation. One thing we know about the result is that the wave propagating in such a medium is dispersive. Considering the fluid body in wet bones as a source providing mass and charge exchanges to solid bone matrix, one can neglect the effect of mass density of fluid component. Furthermore, in bone remodelling processes, it is well-known that the electro-chemical reactions are predominant and, therefore. the effect of density of solid bone matrix on mass transfer may also be neglected. Thus, from (4.20) we have c = dxq - C,cp’+ C& or, in terms of applied electric potential, one has c = G&z, t). c&z, t) = (Doexp[i(kz - tot)]

(4.42)

where the coefficient C, is defined by c

0

=

WC33

-

~2&&3)[i~3dJ

+ (c333

-

p&&)k’+

F(+~ + iG333k + (p&13? - B&k~

iaM+ (bd&’ _ ikC ) 3

(4.43)

As one sees, given cp(z,t) as an external source to start the remodelling process. the mass production term c depends on cp(z,t) and the coefficient C,,. On the other hand, C, depends on the frequency of the wave, and the unknown coefficients d?, C3 and C,-,.+ By performing three experiments with three different frequencies and measuring Co, one can calculate the coeficients d,, C3 and Cj3. Then it becomes possible to judge on which of these coefficients are important.

CONCLUSION

Starting from a two component mixture model, the balance equations and constitutive relations of wet bones are obtained. In the model we proposed, the effects of electro-chemical reactions with mechanical fields have been taken into account. Finally, solving an example problem on a long bone, some remarks made on measuring the material coefficients. More il~ustratingproblems are to follow in separate publications. REFERENCES [l] C. A. L. BASSET, The Biochemistry and Physiology of Bone, Vol. 31, pp. l-69. Academic Press, New York (1971). [2] E. FUKADA andI. YASUDA, J. Pfzys. Sot. (Japan) 12, 1158 (1957), [3] C. A. L. BASSET and R. D. BECKER, Science 137, 1603 (1962). [4] I. YASUDA, J. KYOTO Pref. Univ. Med. (in Japanese) 53, 352 (1953). {S] A. GJELSVIK, J. ~io~ec~anics 6,69 (1974). (61 W. S. WILLIAMS and L. BREGER, Ann. of New York Acad. of Sci. 238, 121 (1974). [7] N. GIUZELSU, An Inoesfigation ofBoneTissues as a LineorPiezoelectric Material, Habilitation Thesis. Istanbul (1976). [8] N. PETROV, Bulgarian Acad. of Sci., Biomechanics 2, 31 (1975). [9] H. DEMIRAY and A. C. ERINGEN, Plasma Phys. IS, 903 (1973). [lo] H. DEMIRAY and N. GUZELSU, Int. J. Engng Sci. 15, 707 (1977). [I l] H. DEMIRAY, J. Tech. Phys. 19, 267 (1978). [l?] 1. C. ANDERSON and C. ERICKSON, Nature 227,491 (1970). (Receiued 16 Februnry

1982)

tThis statement is valid if all the material coefficients appearing in polarization density and stress are determined by some other experiments.