CC BY
30
FORMULATION OF INITIAL BOUNDARY-VALUE PROBLEM AND CONSTRUCTION OF COMPUTATIONAL ALGORITHM IN SIMULATION OF
GROWING BONE TISSUE
A.G. Masich*, S.A. Chernopazov*, Y.I. Nyashin*, E.Yu. Simanovskaya**
* Perm State Technical University, 29a, Komsomolsky Prospect, 614600, Perm, Russia ** Perm State Medical Academy, 39, Kuybishev Street, 614000, Perm, Russia
Abstract: A long-standing experience of clinical observations and methods of orthopedic remodeling developed at Perm School of children's stomatologists specialized in correction of congential hard palate defects have provided solid grounds for mathematical interpretation of this pathology. The objective of the present work is to develop a mathematical formulation and computational algorithm for a initial boundary-value problem of growing bone tissue. The growth processes taking place in bone tissues are the basis for the constitutive relations of the mathematical model. The conjecture that their evolution is governed by mechanical factors has been supported by numerous empirical observations. The initial differential formulation of the problem is changed to a variational one to facilitate its numerical realization. The two-dimensional problem of growth deformation is solved numerically for rectangular region using the finite-element method.
Key words: growing bone tissue, growth strain, differential statement, variational statement, finite element method
A bone has a wonderful ability for reconstruction and transformation. It belongs to functionally adjustable and highly regenerative organs.
S. A. Reinberg
Introduction
Mathematical modeling of biological growth and development of living tissues is actually related to a construction of new continuum models. The interpretation of numerous empirical data led scientists to a general conclusion that tissue behavior including the rate and peculiarities of growth and structure formation is highly sensitive to the action of different forces [1-4].
The first attempts to reveal mechanisms of this process date back to the previous century. Among them there is the world-known Guter-Folkman law, which relates the delay in tissue growth to an increase of compressive load on the growing elements, for example, human's backbones [2]. It was just this law that offered an explanation for the genesis of some forms of scoliosis [5]. A greater part of empirical data is concerned with the action of the force factors on the growth and structure of the bone tissue, which experiences continuous changes and transformations during the whole man's life [2].
The advancements of the continuum mechanics and evolution of its methods [6] have paved the way to an adequate mathematical formalization of the advanced concepts and deep understanding of various complex biological systems whose behavior is governed by interaction of mechanical and involved physicochemical processes. This approach has been
effectively used to develop a number of mechanochemical models for various biological tissues including growing continuum models [1].
In order to construct a mathematical model of tissue growth and evolution there is a need to take into account the observed and experimentally detected dependence of the growth rate on the mechanical stresses [1-3]. The stress dependence of the growth rate can be governed by different mechanisms, which reflect eventually the contribution of the stressstrain state to cellular activity. Serious efforts are currently made to develop physicochemical theories for interpretation of such interactions involving a complicated system of intermediate agents [7].
However, a purely mechanical deformation of biological materials is an essential constituent of growth. For example, a key factor in the growth of the plant cellular wall is the extension due to intercellular pressure, which manifests itself in rearrangement of cellulose microfibrils and partial breaking of their bonds [8].
The fact that the growth process involves mechanical factor suggests that the growth deformation can be controlled by application of external loads. A striking example of how we can artificially stimulate the growth process is the method of destructive elongation, in which significant mechanical stresses are generated in bone tissue with the help of different apparatus designed by G.A. Ilisarov. A high level of tensile stress, which is significantly greater than the normal physiological loading, is a necessary condition for a growth process, which after removal of tensile loads ends with complete calcinosis of the growth zone.
Another attempt to build a general mathematical model and formulate the problem on uniaxial tension of the growing callus has been reported in [9]. However the stated problem proved to be rather complicated for the analysis since it involved many parameters, which at present cannot be defined with sufficient reliability.
The objective of the authors' research is to investigate growth processes in maxilla tissues with pathological changes caused by children's congenital palate cleft, and in particular, the response of separated palatal fragments to the load generated by orthopedic apparatus [10]. The first attempt to develop a simple mathematical model of cleft palate behavior taking into account the growth strains has been made by A.G. Masich, Y.I. Nyashin, and E.Yu. Simanovskaya [11,12]. The computational algorithms have been constructed to determine the desired material constants (growth parameters) from the results of experiments and the optimal force needed to bring down the palate fragments from the nasal to oral cavity [12].
In this study, we present the results of further developments of the mathematical theory of growth strains [11,12] enabling us to estimate with more accuracy the effect of the mechanical factor on the growth deformation.
Mathematical statement of the problem
Here we consider the process of deformation in a growing elastic body, which
occupies the domain Q in R3 with the boundary S, Q = Q U S, S = Sv U Sa . The Sv is the part of the boundary under the prescribed kinematic boundary conditions and Sa is the part of the boundary under the action of the specified forces.
A key factor to modeling the growing living tissue is to ensure a highly rational description of its kinematics. This still unresolved problem is concerned with the development of the large deformation theory for the examined class of phenomena [1].
In the context of the proposed approach we restrict ourselves to a consideration of geometrically linear theory for which we admit the following system of constitutive relations [13]:
Se = N :S,fg = A +M : S,f = fe + fg. (1)
Here ~e is the tensor of elastic strain, f,fe,fg are the tensors of rates of total, elastic and
inelastic (growth) strains, respectively, S is the stress tensor, N is the tensor of elastic parameters, the tensor A characterizes the proper growth of the material (in the absence of
stresses), and the tensor M is responsible for the stress influence on the growth deformation.
Assuming that fe = d ~e we can write the following differential relation: dt
~ = A + M : S + d(N: S). (2)
dt
The rheological relation (2) is nearly coincident with the conventional relation for the viscoelastic Maxwell material differing from the former by nonzero term A, which is responsible for deformation in the absence of stresses.
Equation (2) together with the equilibrium equation
V• S + F = 0, Vr eQ, (3)
geometrical relation
f = 1 (w + VvT ) Vr eQ, (4)
boundary conditions
v = 0, Vr e Sv
S V (5)
n • s = p, Vr e SS
and initial condition S(r,0) = 0 generate a closed system of equations for v,f ,S . In equation (4), v is the vector of the medium growth rate.
Computational algorithm
Let us consider an isotropic growing material with the time-independent growth parameter. The problem is solved under the simplifying assumption for material constants involved in constitutive relation (1). In this case, the rate of growth strain can be written as
f = A~ + MS, (6)
where gS is the metric tensor, and A and M are the growth parameters. The elastic behavior of the medium is represented by
S = , ^ ^ (7)
where n is time independent tensor of elastic properties and N = n-1.
Substituting the equality fe = f-fg in relation (7) and taking into account (6) we
obtain
~ = n : (S - AS - MS). (8)
Define the space V by
V = {v e (C2 (Q))31 v = 0, Vr e Sv}
and suppose that u e V is kinematically admissible velocity field at some time step At.
Differential relation (8) can be represented by two finite-difference approximate expressions [14] in explicit form
st+At =~ t +At ft: (S - Ag - MS1) (9)
and in implicit form
~t+At = g t +At f:(g_ a~ — Mgt+At). (10)
To proceed further we must perform some preliminary transformations. First we express tensor gt+At in expression (10). From the tensor analysis [15] it directly follows that there exists some isotropic tensor of the 4th order CIII = rsg rs such that
CIII : B = B : CIII = B where B is the arbitrary tensor of the second order. Then expression (10) takes the form
A: gt+At = gt + At fl: (g — A~),
~ ~ ~ (11) where A = (CIII + At M ff),
from which we can express g t+At as
~ t+At = a-1 : (~t + At ff : (g _ Ag)). (12)
In order to find approximate solutions to differential problems (2)-(5), (9) and (2)-(5), (12) using available numerical techniques (for example, the Ritz, Bubnov- Galerkin and other methods) it is necessary to construct the variational statement to the initial differential statement.
Using the Jourdain's principle [16], constitutive relations (9) and (12), and extending the space V to
K = (v e (W12(Q))3\v = 0, Vr e Sv}
we arrive at the following variational problems. Problem 1. Find v e K :
At J g (u): ff : g dQ = J G : g (u)dQ. + J F • u dQ + J p • u dSa,
Q Q Q Sa Vu e K, (13)
where G = At ff : (Ag + Mgt) — gt. In the context of implicit scheme (12) the variational problem takes another form. Problem 2. Find v e K :
At J g (u): IT* :g dQ= J G* : g(u) dQ+J F • u dQ+ J p • u dSg,
Q Q Q Sg Vu e K, (14)
Jg
'* _ A * . Ag g . gt
where n* =A : n, G = Atn* : Ag — A : gt.
The value of the function on the boundary in definition of the set K has the meaning of a trace [17].
t T*
In relations (13) and (14), G and G represent some imaginary forces linked with the growth deformation.
Then we find approximate solutions at each time step by applying the finite element method (FEM).
Results of numerical realization and comparison with the beam theory solution
A numerical investigation of the algorithm convergence is made in the context of the model problem on the growth deformation of arbitrarily shaped cantilever beam, defined and solved in our preceding studies [12, 13]. In these works, the configuration of the growing beam loaded at its free end was determined using the beam theory and neglecting elastic strain.
In the present work we concentrate on the linearly elastic growth deformation of the isotropic rectangular domain at plane stress state (PSS). Fig. 1 schematically shows the fixation of the computational domain and its dimensions. The force p(t)=F, t e[0, tend] is applied to the point B. The domain is decomposed into the triangular finite elements and contains 675 nodal points.
Table 1. FEM with consideration for growth (PSS)
Extension Number of time steps u (mm)
explicit scheme 700 0.30195
1000 0.30200
implicit scheme 50 0.30198
100 0.30200
Bending
explicit scheme 700 3.75763
1000 3.83766
implicit scheme 50 3.85875
100 3.86823
Table 2. Beam theory
u (mm)
Extension 0.30223
Bending 4.01881
The calculations have been performed for a material with earlier defined growth
parameters A=0.01 1/ month, M=0.01 mm2/(g • month) [12]. The elastic modulus
6 I 2
A=10 g/mm and the Poisson's ratio V =0.3 [18].
Two loading schemes have been examined: extension under the action of the force F=1g and bending under the action of the concentrated force F=0.1g at the free end. The time of loading is tend=1 month .
The results of numerical calculations are summarized in Tables 1 and 2. Here u is the displacement of the characteristic point A shown in Fig. 1. The results of FEM calculations agree fairly well with the model calculations [14], in which the elastic strains are ignored. Under tensile load the displacement discrepancy at the point A (Fig. 1) is 0.07%-0.09% and under bending load - 3.9%-7%. The disagreement obtained in bending problems is probably due to the shape of the domain: in plane problems the hypothesis for plane sections does not hold in the region of force application.
Y, mm
20
X, mm
Fig. 1. Schematic representation of the computational domain and its dimensions
F
B
A
It should be noted that incorporation of implicit schemes in computations on PC offers a marked advantage of numerical stability over calculations based on the explicit schemes, which prove to be unstable at small number of time steps (<=700).
Thus it is seen that the benefits of plane-stress definition of the problem are manifold: it can treat regions with actual configurations, which is hardly possible in the case of the beam model; it provides a more adequate modeling of the boundary conditions; and, finally, it allows us to take into account the local contact stresses generated by orthopedic apparatus and their influence on the growth deformations.
Conclusions
The present work introduces the results of the latest investigations, which allows us to describe more fully the effect of mechanical factors on the growth deformation. The general mathematical formulation of the growth deformation and the algorithm of its solution have been developed. Two variational statements of the differential problem have been proposed for two schemes of integrating the constitutive relations with respect to time. The two-dimensional problem of growth deformation has been numerically solved for a rectangular region using the finite element approach. The comparison of the obtained results with the results based on beam theory highlights the advantages of the plane problem for our purposes. The results of this investigation suggest that the next possible step to the model refinement is consideration of geometrical nonlinearity.
References
1. РЕГИРЕР С.А., ШТЕЙН А.А. Механохимические модели морфогенеза. Теоретические и математические аспекты морфогенеза, Москва, Наука: 151-161, 1987 (in Russian).
2. РЕГИРЕР С.А., ШТЕЙН А.А. Механические аспекты процессов роста, развития и перестройки биологических тканей. Итоги науки и техники. Комплексные и специальные разделы механики, т. 1. Москва, ВИНИТИ: 3-142, 1985 (in Russian).
3. РЕГИРЕР С.А., ШТЕЙН А.А. Методы механики сплошной среды в применении к задачам роста и развития биологических тканей. Современные проблемы биомеханики, Рига, Зинатне, 2: 5-37, 1985 (in Russian).
4. TABER L.A. Biomechanics of growth, remodelling and morphogenesis. Appl Mech Rev, 48(8): 487-545, 1995.
5. ЕНТОВ В.М. О механической модели сколиоза. Изв АН СССР, Механика твердого тела, 4: 201-208, 1983 (in Russian).
6. СЕДОВ Л.И. Механика сплошной среды, т. 1-2, Москва, Наука, 1984 (in Russian).
7. WEINMBAUM S., COWIN S.C, ZENG Y. A model for excitation of osteocytes by mechanical loading-induced bone fluid shear-stresses. J Biomech, 27: 339-360, 1994.
8. COSGROVE D.J. Biophysical control of plant cell growth. Ann Rev Plant Physiol, 37: 377-405, 1986.
9. ЛОГВЕНКОВ С.А., ШАФИТ С.Е. Математическая модель ткани костного регенерата, полученного в результате дистракционного остеосинтеза. Биофизика, 41(6): 1332-1335, 1996 (in Russian).
10. ШАРОВА Т.В., СИМАНОВСКАЯ Е.Ю. Ортопедический способ устранения врожденного дефекта твердого и мягкого неба у детей с одно- и двусторонней расщелиной.
Методические рекомендации. Пермь, 1983 (in Russian).
11. МАСИЧ А.Г., СИМАНОВСКАЯ Е.Ю., НЯШИН Ю.И., БОЛОТОВА М.Ф. Биомеханическая модель процесса предоперационной ортопедической реконструкции твердого неба у детей при врожденной расщелине. Пермский медицинский журнал, 12(4): 28-31, 1998 (in Russian).
12. MASICH A.G., NYASHIN Y.I. Mathematical modelling of orthopedic reconstruction of children's congenital maxillary anomaly. Russian Journal of Biomechanics, 1: 101-109, 1999.
13. ШТЕЙН А.А. Деформирование стержня из растущего биологического материала, подвергнутого продольному сжатию. Изв АН СССР, Прикладная математика и механика, 59, (1): 149-157, 1995 (in Russian).
14. САМАРСКИЙ А.А., НИКОЛАЕВ Е.С. Методы решения сеточных уравнений. Москва, Наука, 1978 (in Russian).
15. ПОБЕДРЯ Б.Е. Лекции по тензорному анализу. Москва, Изд-во Моск. ун-та, 1974 (in Russian).
16. ГУН Г.Я. Теоретические основы обработки металлов давлением. (Теория пластичности). Москва., Металлургия, 1980 (in Russian).
17. САНЕС-ПАЛЕНСИЯ Э. Неоднородные среды и теория колебаний. Москва, Мир, 1984 (in Russian).
18. JONES M.L., MIDDLETON J., HICKMAN J., VOLP C., KNOX J. The development of a validated model of orthodontic movement of the maxillary central incisor in the human subject. Russian Journal of Biomechanics, 1-2: 36-44, 1998.
ПОСТАНОВКА И АЛГОРИТМ РЕШЕНИЯ НАЧАЛЬНО- КРАЕВОЙ ЗАДАЧИ ДЛЯ РАСТУЩЕЙ КОСТНОЙ ТКАНИ
А.Г. Масич, С.А. Чернопазов, Ю.И. Няшин, Е.Ю. Симановская (Пермь, Россия)
Актуальность настоящих исследований связана с потребностью клиники "объективизировать" накопленные многолетней медицинской практикой способы ортопедической реконструкции врожденной порочно развитой верхней челюсти у детей с целью их усовершенствования. Расщелина твердого неба является объектом биомеханического изучения. Основываясь на том эмпирическом факте, что на скорость роста живой ткани оказывает влияние механический фактор, предложено рассматривать твердое небо как растущую среду, что особенно характерно в детском возрасте. Приводится обзор работ по математическому моделированию роста живой ткани. Ортопедическая аппаратура, используемая при исправлении врожденного дефекта в твердом нёбе, является мощным механическим стимулом ростовых процессов в костной основе нёба. Ранее твердое небо рассматривалось как консольная изгибаемая балка. Предложены алгоритмы определения параметров роста, входящих в определяющие соотношения по результатам экспериментальных данных, и оптимальной силы на незащемленном конце для получения заданной конфигурации неба. В настоящей работе излагается дальнейшее развитие математической теории ростовых деформаций. Приводятся постановка и алгоритм решения краевой задачи для растущей живой ткани. Для исходной дифференциальной постановки строятся вариационные аналоги с целью их дальнейшего применения для отыскания приближенного решения задачи. Методом конечных элементов численно реализована задача для случая линейно-упругого ростового деформирования изотропной области прямоугольной формы в условиях плосконапряженного состояния. Полученные результаты сравнивались с расчетами по балочной теории. В дальнейшем предполагается провести расчеты для реальной геометрии сечения небного отростка в рамках предложенной модели, а также учесть геометрическую нелинейность. Библ. 18.
Ключевые слова: растущая костная ткань, ростовые деформации, дифференциальная постановка, вариационные аналоги, метод конечных элементов
Received 1 September 1999
