CC BY
25
Rutsian Journal of Biomechanics
www.biomech.ac.ru
SIMULATION OF THE BONE TISSUE ADAPTATION PROCESS UNDER THE PROPORTIONAL CHANGE OF THE LOAD
Department of Theoretical Mechanics, Perm State Technical University, 29a, Komsomolskii Prospect, 614600, Perm, Russia, e-mail: [email protected]
Abstract. A new mathematical model of internal adaptive remodeling of cortical and spongy human bone tissues is presented. This model allows predicting in real time scale the variations of bone tissue mechanical properties under proportional change of the load, i.e. without the change of symmetry of elastic characteristics anisotropy. In this case the character of the changed load remains to be close to the physiologic one. The field of application of this model is broad enough; it includes bone tissue adaptation problems in sport exercises, in development of new medical techniques of rehabilitation after the trauma or surgical intervention. The model uses the classical constitution relation of porous elastic media by Hegedus and Cowin, and a new tensor kinetic equation of internal remodeling that establishes a linear relationship between the rate of bone tissue stiffness tensor and a scalar deformational stimulus of the bone cells activation. In this equation the proportionality factor that characterizes the load sensitivity of the bone tissue, is the tensor of the fourth rank. In order to justify the new kinetic equation, the comparison of this new phenomenological equation of porous elastic adaptive media with the analogous equation that follows from the theory by Hegedus and Cowin has been carried out. It has been shown that under the assumption of linear relationship between the bone cells activation and the deformational stimulus, all the independent components of the sensitivity tensor may be determined in an analytical form as functions of the structure parameters. For the spongy bone tissue the analytical expressions of these functions have been derived. The test calculation of the adaptation process in a homogeneous porous bar under the increase of uniaxial compressive load is presented. The monotonous character of the time variations of all the parameters of the process qualitatively corroborates the validation of the developed model.
Key words: adaptive remodeling, real time scale, constitutive relations, kinetic equation, adaptive sensitivity tensor, cells activity, phenomenological equation
The nature of bone tissue adaptive remodeling is one of fundamental problems of biomechanics. As it can be seen from the overviews [1-3], development of the theory and simulation techniques of bone remodeling attracts the attention of many researches. The adaptive properties of the bone tissue govern the course of many clinically important processes of recovery of bone tissue mechanical characteristics, e.g. post-operational or posttraumatic rehabilitation. Developed by S.C. Cowin with co-authors in 1970th, the theory of elastic adaptive porous media [4-6] is the theoretical basis of all the models of bone tissue mechanical adaptation. These models usually assumed that the functional stimulus for bone cells adaptive activity is strain [5, 7, 8], stress [9, 10], or damage [11, 12], or strain and damage simultaneously [13]. However, the constitutive equations of theory by Cowin and
Yu.V. Akulich, R.M. Podgaets
Introduction
Hegedus involve the second rank tensor A(e) and the fourth rank tensor C(e). These tensors depend on the deviation of matrix volumetric fraction from its homeostatic value; therefore they are hard determinable and call for prolonged clinical experimental investigations. This is the reason of the fact that in literature up to now does not occur data about numerical values of the components of these tensors. In the present paper, in the development of idea on composition of the phenomenological equation for spongy bone tissue [14], we postulate the kinetic equation of internal remodeling for cortical and spongy bone tissues, and mathematically determine the form of its phenomenological equation. The used kinetic equation involves the fourth rank tensor as a function of the parameters of structure and activation of bone cells. This tensor characterizes the bone cells sensitivity to deformation of the bone element. The components of the sensitivity tensor, under some light enough assumptions, may be derived analytically. It takes an opportunity to simulate the bone adaptation in real tirpe scale, without long-duration clinical experiments.
Methods
Basic statements
The subject of the study is an adult human bone tissue, and the time variations of bone mechanical characteristics under the physiological load deviation from its homeostatic value are considered. The process of bone tissue continual renewal under the homeostatic physiological load does not included in the model. We assume that initial configuration is the configuration of homeostatic state, the strains in the bone tissue are small and follow the constitutive equation of the elastic adaptive media [5]:
cr = td(Z):e9 (1)
where a, s are stress and strain tensors, C is the fourth rank elasticity tensor, and £ is volumetric fraction of bone matrix at the current instant of time. During the process of bone adaptation strains and elastic characteristics are varying in time, and these variations are caused by the action of two different by nature factors: the mechanical factor from the external forces, and the biological factor from the activity of bone cells. Therefore the stress tensor in Eq. (1) is a function of two independent variables and s(t)9 so the material time derivative of the stress tensor is:
¿r = C' + :e, (2)
where the following designation is used
(3)
The elastic characteristics rate tensor C is dependent on the rates of matrix resorption and reposition by the bone cells, and on the magnitude of mechanical load stimulating their activity. As follows from the literature data on the mechanism of bone tissue adaptation [15-17], the bone cell responds to the stain variation in its neighbourhood, i.e. the deformational remodeling stimulus is valid. At the present time there are no data about the effect of the form or/and the volume on the bone cells activity, therefore we adopt the strain intensity si as a measure of cell neighbourhood strain. The kinetic equation of bone
remodeling may be expressed in the form that reflects the immediate effect of the remodeling stimulus on the tensor of elastic characteristics rate:
In Eq. (4) R is the fourth rank tensor of adaptive sensibility that characterizes the rate
of time variations of the bone tissue elastic characteristics in different directions, and s^10™ is
the homeostatic value of strain intensity. The substitution of Eq. (4) in Eq. (2) gives us the phenomenological equation of the linear elastic porous adaptive media
Sr^fa-e^Rie + C' (5)
The structure of the kinetic equation (4) follows from theory by Cowin and Hegedus [5], as the corresponding phenomenological equation that follows from this theory has the same form as Eq. (5):
r
o - (a(e) 4- A(e): <?)
C(e) + e™ 8e
:s+C(e):s , (6)
where a(e) + A(e): s = e is a scalar, e = is the deviation of matrix volumetric fraction
4 from its homeostatic value a(e) and A(e) are the experimentally determined functions.
The variation rate of bone tissue elastic characteristics is determined by the change rate of its structure that in turn depends on pores surface area with active bone cells. Therefore
it is naturally suppose that the tensor R has to be depend on the geometric parameters gi,
i - 157i characterizing the variable structure of the bone tissue. In this sense the tensor R(gi) is analogous to the tensor A(e) which has to be either specified or determined by prolonged
clinical experiments. The main advantage of the new kinetic equation (4) is the tensor R(gt)
under some assumptions may be determined analytically. The technique of such determination for the spongy bone tissue, and the test calculation are presented below.
The structure of the tensor R
The matrix form of Eq. (4) is ([C km] is symmetric matrix):
[Ckm ] = (st - e?om ],K,M = 16. (7)
As follows from Eq. (7), matrix [RKM ] is also symmetric one. Therefore the tensor R rosy involve no more than 21 independent components for general kind anisotropy of the ncniDdeled bone tissue. It is a fact (see [18]) the more is the number of material directions that ire equivalent by the elastic characteristics the larger is the number of zero elements in the
r^jix of elastic characteristics [CKM ] and hence in matrixes [C km j and [R KM ], and the less :;5 ne number of their independent components. Therefore the structure of the matrix [c km j ¿r i tfce number of its independent components in the current time are the characteristics of the roce tissue structure at the same instant of time. The main assumption used below is the s>Tnmetn of the elastic characteristics anisotropy in the bone tissue does not change during me remodeling process. Apparently the remodeling of such kind may take place under only quantitative change of the healthy bone tissue loading, when the character of the changed load remains close to physiological one, i.e. the loading is proportional or nearly proportional. The examples of such loading are the sport exercises, or the therapeutic limbs loading in order to recover the bone tissue mechanical properties after the enforced immobilization caused by rauma or surgical intervention.
Assuming the spongy bone tissue to be an isotropic linear elastic material, we can write the elastic characteristics rate tensor in a form that is well known for an isotropic elastic media with the Lame constants X and ju (see [18]):
X + 2ju X X X X + 2 ju X 0 X X X + 2 ju
[cKM\-
M 0 0
0 0 M 0
0 0 A
(8)
It is follows from Eqs. (7) and (8) that the matrix of sensitivity coefficients [RKM ] is of the same structure with three independent elements:
[Rrm]-
R\2 Rn
Rn Rn Ru 0
Rn Rn Ru
R 44 0 0
0 0 R 44 0
0 0 R 44
(9)
Hence, for the tensor R determination we have to specify its three components: i?nn, R\\22 9 ^1212 (here the known convention on correspondence between two-index and four-index component denotations of the fourth rank compliance tensor [19] is used).
Calculation of the tensor R components
The changeover to technical elastic moduli in Eq. (8), in accord with Eq. (7) gives us the following relationships:
d f M1-v ) 1
dt L(l + vXl-2 v)J
d f ÇEv \
dt Ul + vXl-2v)J
df ÇE ) = R
dt\\ + 2v j
n I hom \
= R\m\ei-ei j, 1212 (pi~si )>
(10) (11) (12)
where the bone tissue elastic modulus E and Poisson's ratio v at the continuous media level are the functions of the volumetric fraction of bone matrix c. In can be seen from Eq. (20) that the variable is explicitly involved as a multiplier in all the components of the elasticity tensor. Calculating of derivatives in Eqs. (10)-(12) results in the expressions that specify the components Run, Rn22, Run'-
R\\\\ =■
r-/1 x JdE \ _ôv E{ l-v) + 4 —(1 -v)-E—
(l-2vXl + v)+^^(l + 4vXl-v)
Si ~ £,
hom
[(l + vXl-2vff
,(13)
122 ~
Ev + Z
rdE
—v + E— ^
(l-2vXl + v)+^v|^(l+'4v)
hom
[(1 + KX1-2,)]=
^1212 ~
E + c —
(l-2v)-2 $E
8v
hom
(1 + 2 v)2
(14)
(15)
The model of the spongy bone tissue structure
Follow to paper [20], the pores of the structure are assumed to be spherical with radius r and the pores number n in the volume unit of bone tissue is constant in time. Therefore the
values r and n are the components of the parameters vector gi, i-\,n. The volumetric fraction of bone matrix (in 1 mm3 of bone tissue) is
£ = 1-4/3 nn?
(16)
As the first approximation we assume the first time derivative of the pore radius is a linear function of the adaptation stimulus
r = -a(£i-e?om), (17)
where a is a factor that depends on the bone cells activity [21], hence with Eq. (16) we have an expression for the first time derivative of the volumetric matrix content:
i = -47rnar2(£i-4om). (18)
Substituting of Eq. (18) in Eqs. (13)-(15) yields the analytical expressions of the
tensor R components that take into account the bone cells activation (the factor a) and the structure variation:
E{l-v)+4
R = 4nnr a-
^-(l-v)-E^-
\
J)
(l-2vXl + v)+£E—(l + 4vXl-v)
Ry 122 = 47tnr a ■
r
Ev + E,
dE „dv — v + E—
[(l + vXl-2v)r
\\
,(19)
dv
j j
(l-2vXl + v)+^v—(l + 4v)
i
[(l + vXl-2v)]=
E + Í
R\2\2 = 4xnr a-
5£
H.
(l-2v)-2
dv
dt
(20)
(21)
(1 + 2 v)2
The derivatives dE/d% and dv/d% in Eqs. (19)-(21) can be specified from the Ennerical experiments with 3D-specimens of porous elastic material.
Results
The verification of the presented model has been performed by numerical solution of the unidimensional adaptation problem for the homogeneous porous elastic medium under uniaxial compression stepwise increased by 5%. Such nearly uniaxial compression is typical e.g. for the spongy bone tissue in the subchondral zone of the human hip head. Mechanical characteristics of the medium assumed to be equal to spongy bone tissue ones. The pore
Fig. 1. Time variations of bone remodeling parameters.
surface is spherical. The constitutive relation (1), kinetic equation (4) and phenomenological equation (5) for this case may be written in the following form:
a = $Ee\ (22)
^ = (23)
dt
& = R(£-£hom)e + ZE£. (24)
In accordance with data obtained by the usage of the block model [22], the stiffness of the homogeneous porous medium may be approximated as a cubic parabola:
= (25)
where Em is the elastic modulus of bone matrix. The calculation of the derivative at the left hand side of Eq. (23), on account of Eqs. (18) and (25) gives us
R = \27znr2%2aEm. (26)
As during the stepwise load change & = 0, the following differential equation relative to strain yields from Eq. (24):
R(e~ehom)s + ^Ee = 0. (27)
This equation was numerically integrated by the Euler's scheme under the following initial conditions and process parameters:
r0 - 0.33 mm, s0 = 0.01395 , E0 = 791 MPa ,
w = 5, Em = 12660MPa, a -0.09 mm/day, ehom = 0.01329, where E0 is the initial Young's modulus value. The factor a value was determined according :: data of [23].
The results of the numerical integration are presented in Fig. 1. The strain-time plot is shown in Fig. la. Under 5 % load increase the compressive strain stepwise rises, and its recovering from the initial value £0 to the homeostatic one goes on 50-60 days.
Simultaneously the following changes of the structure and physical characteristics come about: the pore radius decreases (Fig. Id), the volumetric matrix content rises (Fig. le) and consequently, Young's modulus (Fig. lb) and stiffness (Fig. If) also are on the increase. These data correlate with present-day concepts of the character of bone tissue adaptation [2]. The curve in Fig. 1c shows that the sensitivity of the spongy bone tissue with respect to the vxam stimulus is not a constant value (as it is stated e.g. in [24]), but changes during the rem-c-cieling process. The monotonous character of changing the adaptation process parameters :: iis steady-state values qualitatively confirms the correctness of the presented model, as the ptotic solutions stability in the problems of adaptive poroelasticity is proved ±&9ceucaliy [6]. To estimate the accuracy of the model, it is necessary to compare its results - ± the experimental data.
Conclusions
The model presented above allows forecasting in real time scale the changes of the rone tissue mechanical characteristics during the internal adaptation process of spongy and cordcai human bone tissues. This model is applicable for investigation such adaptive processes where the character of changed load remains close to physiological one, i.e. the loading is proportional one or nearly to it. However the field of this model application is broad enough, it includes bone tissue adaptation problems in sport exercises, in development of new medical techniques of rehabilitation after the trauma or surgical intervention.
References
1. РЕГИРЕР С А, ШТЕЙН A.A., ЛОГВЕНКОВ C.A. Свойства и функции костных клеток: биомеханические аспекты. Современные проблемы биомеханики. Вып. 10. Механика роста и морфогенеза, Москва, Издательство Московского университета, 174-224, 2000 (in Russian).
2. TABER L.A. Biomechanics of growth, remodeling, and morphogenesis. ASME Appl Mech Rev, 48(8): 487-545, 1995.
3. COWIN S.C. Structural adaptation of bones. Appl Mech Rev, 43(5-2): S126-S133,1990.
4. COWIN S.C. and HEGEDUS D.H. Bone remodeling I: theory of adaptive elasticity. J Elasticity, 6(3): 313-326, 1076.
5. HEGEDUS D.H. and COWIN S.C. Bone remodeling II: small strain adaptive elasticity. J Elasticity, 6(4): 337-352, 1976.
6. COWIN S.C., NACHLINGER R.R. Bone remodeling III: uniqueness and stability in adaptive elasticity theory. J Elasticity, 8(3): 285-295,1978.
7. COWIN S.C., MOSS-SALENTIJN L., MOSS M. L. Candidates for the Mechanosensory System in Bone. Trans ASME. J Biomech Eng, 113(2): 191-197. 1991.
8. AKULICH Yu.V., DENISOV A.S., PODGAETS R.M., AKULICH A.Yu. Computer simulation of the elastic modulus of hip head bone tissue during the hip joint loading variations. Russian Journal of Biomechanics, 1: 53-61, 1999.
9. COWIN S.C. Bone stress adaptation models. Trans ASME. J Biomech Eng, 115(4): 528-533,
1993.
10. CARTER D.R. Mechanical loading histories and cortical bone remodeling. Calcif Tissue Int, 36: 19-24, 1984.
11. PRENDERGAST P.J., HUSKES R. Damage-adaptive bone remodeling: exploring the link between micro-damage and strain stimulus. 2nd World Congress Biomech, 2: 51, Amsterdam,
1994.
12. PRENDERGAST P.J., TAYLOR D. Prediction of bone adaptation using damage accumulation. J Biomech, 27(8): 1067-1076, 1994.
13. HAZELWOOD S.J., MARTIN R.B., RASHID M.M., RODRIGO J.J. A mechanistic model for internal bone remodeling exhibits different dynamic responses in disuse and overload. J Biomech, 34: 299-308, 2001.
14. AKULICH Yu.V., PODGAETS R.M. The phenomenological model of adaptable adult spongy bone tissue. Russian Journal of Biomechanics, 1-2: 53-57, 1998.
15. LANYON L. E. Functional strain as a determinant for bone remodeling, Calcif Tissue Int, 36: S56-S61, 1984.
16. LANYON L. E. Osteocites, strain detection, bone modeling and remodeling. Calcif Tissue Int, 53 (Supl.1): S102-S106, 1993.
17. RUBIN С. Т., LANYON L. E. Osteoregulatory nature of mechanical stimuli: function as a determinant for adaptive bone remodeling. J Orthop Res, 5: 300, 1987.
18. MASE G. E. Theory and Problems of Continuum Mechanics. McGraw-Hill, New York, 1970.
19. АШКЕНАЗИ E.K., ГАНОВ Э.В. Анизотропия конструкционных материалов. Справочник. Ленинград: Машиностроение, 1960 (in Russian).
20. MARTIN B.R., HAYNES R. R. The investigation of bones substructure using megahertz sound and a porous model. ASME Publ, N70-WA/BHF-11. 11 p, 1970.
21. MARTIN R. B. The effects of geometric feedback in the development of osteoporosis. J Biomech, 5:447-455, 1972.
22. McELHANEY J. H., ROBERTS V. L. Mechanical properties of cancellous bone. AIAA Paper, N 71-111. 8 p, 1971.
23. SOTIN A. V., AKULICH Yu. V., PODGAETS R.M. The model of cortical bone tissue adaptive remodeling. Russian Journal of Biomechanics, 5(1): 24-31, 2001.
24. ADACHI Т., TSUBOTA K., HOLLISTER S.J. Trabecular surface remodeling simulation for cancellous bone using microstructural Voxel finite element models. J Biomech Eng, 123: 403-409, 2001.
МОДЕЛИРОВАНИЕ ПРОЦЕССА АДАПТАЦИИ КОСТНОЙ ТКАНИ В РЕАЛЬНОМ МАСШТАБЕ ВРЕМЕНИ
Ю.В. Акулич, P.M. Подгаец (Пермь, Россия)
В статье представлена математическая модель процесса внутренней адаптационной перестройки спонгиозной и кортикальной костных тканей человека, позволяющая прогнозировать изменение механических свойств в реальном масштабе времени. Модель использует классическое определяющее соотношение пористой линейно - упругой среды, полученное Hegedus and Cowin [5], и новое тензорное кинетическое уравнение внутренней перестройки, устанавливающим линейную связь между скоростью изменения тензора жесткости костной ткани и скалярным деформационным стимулом активности костных клеток. В качестве коэффициента пропорциональности, характеризующего адаптационную чувствительность костной ткани к изменению физиологической нагрузки, в этом уравнении используется тензор четвертого ранга. Приводится обоснование структуры построенного кинетического уравнения путем сопоставления нового феноменологического уравнения пороупругой приспосабливающейся среды с аналогичным уравнением, следующим из теории Hegedus and Cowin [5]. Показано, что для процессов внутренней перестройки костной ткани, в которых не изменяется вид симметрии анизотропии ее упругих свойств f характер измененной нагрузки остается близким к физиологическому, то есть осуществляется пропорциональное нагружение или достаточно близкое к нему), и в предположении линейной зависимости активности костных клеток от деформационного стимула все независимые компоненты тензора чувствительности уоглт быть определены как функции параметров структуры. Для спонгиозной костной тузни получены аналитические выражения этих функций. Приводятся результаты тестового расчета процесса адаптации однородного пористого стержня в реальном масштабе времени при увеличении одноосного сжатия. Полученный монотонный характер установления всех параметров процесса адаптации является качественным подтверждением обоснованности разработанной модели. Библ. 24.
Ключевые слова: адаптационная перестройка, реальный масштаб времени, определяющие соотношения, кинетическое уравнение, тензор адаптивной чувствительности, активность клеток, феноменологическое уравнение
06 March 2002
