CC BY
42
MECHANICAL IMPEDANCE OF BIOLOGICAL SOFT TISSUES: POSSIBLE
MODELS
E.M. Timanin, E.V. Eremin
Department of Radiophysics in Medicine, Institute of Applied Physics, Russian Academy of Sciences, 46, Ulyanov Street, 603600, Nizhny Novgorod, Russia
Abstract: Theoretical expressions for the layer impedance characteristics are written when the layer is rigidly connected with a base and when an indentor vibrating on its surface does not produce shear stresses but produces different profiles of normal pressure. By means of the computer-based experimental set-up the impedance characteristics spectra of the homogeneous gelatinous layer are registered as well as their changes during changing the indentor diameter. Calculations of impedance characteristics are conducted in the models with "uniform", "parabolic" and "hyperbolic" pressure profiles under the indentor and comparison of calculation results with experiments is performed. The most suitable model for experiment description appeared to be the model with the "uniform" profile.
Key words: biological soft tissues, mechanical impedance, layered systems modeling, hardware-software experimental set-up
Introduction
Impedance characteristics [1] of biological soft tissues are determined in experiments in which a small, hard, vibrating indentor is pressed in these tissues and indentor kinematics parameters (displacement (U), velocity (V) or acceleration (A)) and resistance force to tissue deforming (F) are measured. For the complete description of tissues behavior in these experiments one can use real and imaginary parts of any of three equal characteristics: complex stiffness K = -F/U, complex mechanical impedance Z = -F/V and complex inertia M = -F/A, or any pair of independent values, in particular, real parts of stiffness and impedance (ReK and ReZ). Studies of impedance characteristics of biologic soft tissues have been performed for a long time [2-10], but only recently they received a new impetus owing to the development of modern computer measurement facilities and data processing [11-14]. Such studies began as far back as 40-ies in connection with the problem of different contact sensors matching with the human body surface [2-4]. A bit later studies of relationship between impedance characteristics and tissue conditions as well as development of ways of evaluation of these conditions on the basis of impedance characteristics measurements began [5-14]. New lines within the framework of these studies are development of the method to reconstruct mechanical parameters of multi-layered tissues on the basis of spectral impedance measurements (that is to say on the basis of frequency dependencies of impedance characteristics) and development of the method of continuous monitoring of tissues viscoelastic parameters with high time resolving power on the basis of single-frequency impedance measurements [11-14]. These methods offer new opportunities of tracking the changes of tissues viscoelastic characteristics, primarily muscles, in the course of different physiological and pathological processes and in the course of the development of response to different test influences. Thereby, new opportunities are being offered for biomechanical and medical study of human neuromuscular system, for example, when studying the mechanism of the motor control or when studying an action of different drugs.
Development of mathematical models of impedance characteristics of biological soft tissues, having, obviously, an independent scientific interest, is the most important stage of development of the method to reconstruct mechanical parameters of multi-layered tissues and the method of continuous monitoring of tissues viscoelastic parameters with high time resolving power. At present several such models are known [8,11,14-17]. The most complete one is the three-layered model proposed on the basic of the experiments on human forearm tissues [11]. However this model can be hardly used for interpreting the results of experimental studies of impedance characteristics of other body areas, which are very diverse in the construction and mechanical characteristics, because of the difficulty of its parameters identification, requiring specialized software. It seems possible to use for this purpose other models ("models with the pressure source of vibrations") which are less strict, but greatly less labor-consuming in calculations, and based on the approach taken when solving the Lamb's problem [14,18-20]. These possibilities are demonstrated below by comparison of calculations in the single-layered model of such a type with the experiment data obtained for the homogeneous gelatinous layer.
Models
To construct models with the pressure source of vibrations, similarly to the model [11], the following approach is used [18]. Firstly, consideration is restricted to the axially symmetric case and the general solution of equations for the acoustic field in the linear elastic medium in Hunkel images is found. Secondly, the boundary conditions describing the considered layered object more exactly are assigned and unknown functions A and B entering the general solution are found. Thirdly, by means of the Hunkel inverse transform the stiffness K (impedance Z) of the considered object is found as the relation between the force P = - F applied to the indentor and its displacement U (velocity V). The specific feature of models with the pressure source of vibrations is the fact that on an external surface of the object the simplified boundary conditions are taken used for the first time by Lamb at the statement of the vibrating indentor problem on the surface of the homogeneous half-space [18]. They are the condition of normal pressure definiteness and the condition of shear stress absence on the whole surface, including the region under the indentor. The latter condition is interpreted as the condition of the indentor slippage and, in principle, its implementation can be provided in the experiment by special procedures. The pressure profile under the indentor p(r) in the models of this type must be chosen such a way to conform to experiments.
The general solutions for Hunkel images of components of the displacement vector and the stress tensor in the medium are got as follows [14]:
Uz0 (k, z) = -k 2 e K(ZA1 (k) - k 2 e ~K'zA2 (k)+ k ¡e K/ZB1 (k) -K/ e ~K/ZB2 (k), U} (k, z) = k k te KtZA1 (k) - k k te ~K'ZA2 (k) - ke K/Z B1 (k) -ke ~K'ZB2 (k), a 0Z (k, z) = -2|k 2 k te KtZA1(k) + 2|k 2 k t e"K t Z A2 (k)+
+ |(k 2 + k 2)e K lZB1(k ) + |(k 2 + k 2)e "K lZB2(k), a }Z (k, z) =|k (k 2 + k 2)e KtZA1(k) +|k (k 2 +k 2)e "K tZA2(k )-
- 2|kcle K l ZB1 (k ) + 2|k k le "K l ZB2 (k ).
9 9 9 9 9 9
Here k is a parameter of the Hunkel transform; parameters kt = k - kt and k / = k - ki are
2 2 I 2 2 2 I 2
defined by the wave numbers of shear and longitudinal waves kt = ra / ct , ki = ra / c/ , where co is the circular frequency of indentor vibration, ct = |/p and c2 = (A,+ 2|)/p are
velocities of shear and longitudinal waves, defined by density p and by the Lame constants A and f.
The system of equations for unknown functions A;, A2, B;, B2, corresponding to the single-layered object with the boundary conditions of complete adhesion on the lower base z = H and with the discussed approximate boundary conditions on the upper surface z = 0, is got from (1) in the following form:
- 2|k2 k tA1 + 2| k 2 k t A2 + |(k 2 + k 2) B1 + |(k 2 + k 2) B2 = - p(k),
2 2 2 2 |k (k + k 2)A1 + |k(k + k 2)A2 - 2|kK/B1 + 2| 1k /B2 = 0,
-k2 eKtHA1 -k2e"KtHA2 + k¡eK/HB1 -k¡e~k/hB2 = 0,
kkt eKtHA1 -kkte"KtH A2 -keKlHB1 -ke"KlHB2 = 0 . (2)
The function p(k) here is Hunkel image of the pressure profile p(r) on the external surface of the layer. When the flat round indentor of radius a vibrates on the layer surface in the area r < a, it is possible to assume [14,19] that distribution of pressure under the indentor is uniform
/2 /99
na , "parabolic" p(r) = 2[1 - (r/a) ]P / na or "hyperbolic"
p(r) = P/2naVa2 - r2 and that outside of the indentor there is no pressure on the surface in all cases.
Expressing the unknown functions, entering the linear system of algebraic equations (2), through its determinants, substituting these expressions in formula (1) for normal displacement of the layer upper surface, applying the inverse Hunkel transform and determining the indentor displacement through the displacement of the layer surface under the
a
indentor averaged over its area [19] U = -1r [uz (r,0)2nr dr, for the layer stiffness we
na 0
obtain the following formula
P 1
K = —=---. (3)
U °f K/ (D13 + D14) - k 2( D11 - Dn)R(k) dk 0 (k2 +k2)(D13 - D14) - 2k2K t(D11 + D12) Here determinants of the third order D1j- are algebraic complements of elements of the first line of the principal determinant of system (2). Function R(k) entering formula (3) is defined by pressure under the indentor and has the following form
R(k) = - (4)
kna |
if the pressure under the indentor is uniform;
R(k) = - 8J1(ka)Jfa), (5)
k n a |
if the pressure under the indentor is distributed according to the "parabolic" law;
R(k) = - J1(ka)s2n(ka), (6)
kn a |
if the pressure is distributed according to the "hyperbolic" law.
Expression (3), in which the determinants correspond to the system of equations (2) with various variants of function R(k) (4) - (6) will be used below for numerical calculations and for approximation of experimental data with the purpose to chose the most adequate model. The variants of models corresponding to the different functions R(k), will be called for compactness A - models, PA- models and GA-models, respectively.
Fig. 1. A scheme of software-hardware complex for the spectra study of impedance characteristics of
biological soft tissues and their phantoms.
1 - power amplifier (type 2707);
2 - vibration exciter (type 4801T);
3 - impedance head (type 8001);
4 - removable indentor;
5, 6 - amplifiers of signals of sensors (type 2626);
7 - computer with CT4170 soundcard;
8 - an object of studies.
Experiment
For the experimental research of frequency dependencies (spectra) of the impedance characteristics of biological tissues and their phantoms, the specialized hardware-software complex [14] was constructed which provides reception of spectra of the impedance characteristics in the electronic form, in which they can be easily used for the further processing, in particular, for identification of the model of the specific object. The experimental set-up based on vibration-exciting and vibration-measuring equipment of the Bruel & Kj? r is used as the basis of the complex [7,8]. In the new set-up (Fig. 1), the processing of signals is carried out in the computer with the help of the specialized software working in Windows 95/98. The input of signals is realized with the help of CT4170 soundcard of the Creative Labs. The program shell allows to determine, to measure and to save on the hard disk the impedance characteristics spectra of the researched object in a range up to 512 Hz. Time of reception of one spectrum is 1 second, frequency resolving power is 1.22 Hz. There is an opportunity of averaging of any number of received spectra. Compensation of mass attached to the force gauge is carried out in each experiment before measurements, that is compensation of the force gauge accelerometric sensitivity. To do this the signals from gauges corresponding to vibrations of the indentor in air are entered in memory of the computer and during measurements the appropriate amendments are done. Besides, the calibration of the system is carried out before measurements by placing of a load of the known weight on the indentor. The appropriate signals are also entered in memory of the computer and are used during measurements for scaling of the determined impedance characteristics. In the mode of measurements in windows of the program shell the frequency dependencies of real (ReM) and imaginary (ImM) parts of complex inertia in grams or frequency dependencies of real parts of complex stiffness (ReK) in N/m and complex impedance (ReZ) in N-s/m are displayed. These values can be saved on the hard disk and can be used for the further processing. Verification of work of the new complex was carried out in several experiments [14]. Firstly, the impedance characteristics corresponding to a testing
load attached to the indentor were registered. Secondly, the synchronous measurements of impedance characteristics of a human relaxed forearm were carried out by means of new complex and by means of the spectra analyzer of type 2034, which was connected in parallel to computer.
The special series of measurements on the homogeneous gelatinous layer of 30 mm thickness was conducted on the described complex. The values ReK and ReZ were registered by means of three indentors with diameters 6, 10 and 16 mm. Each measurement was conducted under steady-state pressing of the indentor in the object on 1 mm. Averaging of 20 spectra was conducted in the course of each pressing. The density of gelatinous sample r « 1008 kg/m3 was determined by measuring the mass of the sample and its volume as well as the velocity of longitudinal waves in the sample c/ « 1500 m/s was determined by measuring time of spreading the ultrasonic pulse from the surface up to the base and back. The registered experimental data were read in MathCAD-files for calculations of impedance characteristics. Looking over the rheological parameters of models was conducted there for best approximation of experimental data. Experimental curves will be given below together with the results of numerical calculations.
Numerical calculations
The numerical calculations in models were carried out by means of MathCAD 6.0 directly by the formula (3), taking determinants from the system (2). The account of viscous properties of the layer material was carried out by replacement of its elastic parameters by the complex operators corresponding to a viscoelasticity type that can be done when solving any problem on the stable vibrations of linear viscoelastic body [21]. As a model of viscoelastic behavior the elementary Foigt body was chosen [22]. According to this model the Lame constants should be taken as: | = |0 + zra^, A = A0 + zra£, where fi0 and A0 are the static
elasticity modules, and r and % are the viscosity modules. Just this expression for f was accepted initially when calculating ct and kt, which as a result appeared to be complex. When calculating 1 the real value c/ was accepted initially which was taken from the experiment c/ ~
1500 m/s. Analyzing the complex expression for cl = ^(A+ 2|)/p , it is possible to find out,
that with reduction of frequency its real part tends to the value c/ = -J(A0 + 2|0)/p , and its imaginary part tends to zero. The validity condition of this limiting process will be the condition ra«Qlim = (A0+ 2|0)/(^ + 2^), which probably should be valid at frequencies below 1 kHz, where the measurements were carried out.
Before calculation of integrals in MathCAD the research of integrands was carried out and the range was determined, in which they remain essentially distinct from zero. The upper limit of integration was chosen of the order 7500 ^ 10500 that lies outside this range. As integrands have, when k is small, a rather sharp splash, if the viscosity of the material is small, it is necessary to divide the interval of integration into two subintervals. The 1-st one is rather short (up to k = 500 ^ 2000) and contains the splash, the 2-nd one is longer and the functions slowly damp here. When calculating integrals it was verified whether the results dependent on the upper limit of integration and on the method of dividing the interval of integration on subintervals.
Approximation of the experimental data was carried out by variation of viscoelastic parameters of models and by visual comparison of calculation results and experimental curves displayed on one graph. In all cases, first of all, the elastic parameters were selected to fit the level of the low-frequency plateau of stiffness curve, and then the viscosity parameters were selected to fit the level of the impedance curve in the range of middle and high frequencies.
(a)
400
0
-400
E
§ -800
PR 1200
■1600
(c)
400 0
E -400 ^ -800 £-1200 -1600
(e)
400 0
E -400 ^ -800 £-1200 -1600
1
100
200 300 f, Hz
400 500
H= 5 kPa,
9 -
100
200 300 f, Hz
1 1
9 ^^
(b) 2 1,5
I 1
z
^ 0,5 &
0
400 500
H= 5 kPa,
0
r = 0.2 Pa-s.
<d) 2 1,5
E
"in 1 z
^"0,5
C£
0 0
r = 3 Pa-s. (f) 2
1,5
E
"in 1
"¡3 0,5 Pi
100
100
„1—
100 200 300 400 500
f' Hz /u = 4 kPa,
0 100
r = 0.8 Pa-s.
200 300 f, Hz
400 500
200 300
f, Hz
200 300
f, Hz
400 500
400 500
0
0
0
Fig. 2. Experimental (1) and calculated (2) impedance characteristics of the gelatinous layer. Graphs (a, b) correspond to A-model, (c, d) - to PA -model, (e, f) - to GA -model. Parameters of models, except the ones given on graphs, are d = 10 mm, H = 30 mm , p =1008 kg/m3, ct = 1500 m/s.
Results and discussion
The comparison of various models by opportunities of reproducing properties of the homogeneous gelatinous layer gives the following results. The best conformity of calculations and experiments is observed in the A - model (Fig. 2). The model reproduces the low-frequency plateau of the curve ReK(f), the high-frequency plateau of the curve ReZ(f) and the qualitative picture of the layer resonances. Moreover the reproduction of all these features of curves appeared to be valid with the fixed parameters of the model for various diameters of the indentor (Fig. 3). The high-frequency fall of ReK(f) curve in this model as well as in the other considered models, is reproduced more abrupt in comparison with the experiment. Probably, it is connected with the accepted approach "of the pressure source of vibrations".
100
200 300 f, Hz
400
400 0
£ -400 ^ -800 £-1200 -1600
600 -600 g-1800 ^-3000 £-4200 -5400
100
200 300 f, Hz
100
200 300 f, Hz
500
d = 10 mm
1 I Tvf^Ujju
d = 16 mm
2 1,5 1
400 500
"S 0,5 oí
0
2 1,5 1
«¡g 0,5 oí
0
400 500
100
100
100
200 300 f, Hz
400 500
200 300 f, Hz
400 500
200 300 f, Hz
400 500
Fig. 3. Experimental (1) and calculated (2) in A-model frequency dependencies of impedance characteristics of the gelatinous layer corresponding to different diameters of the indentor. Parameters of the model: H = 30 mm, p =1008 kg/m3, f = 5 kPa, ] = 0.2 Pa-s, c, = 1500 m/s.
0
0
0
0
0
0
The important feature of A - model is that the conformity to the experiment of the level of losses ReZ(f) in the middle and in the upper parts of the used frequency range turns out to be valid automatically after taking very small values of viscosity h and after selecting the elasticity module of the layer ¡ for reproduction of the level of the low-frequency plateau of stiffness ReK(f). Variation of viscosity in the range of values 0.1 ^ 1.0 Pa-s practically does not influence the level of losses, and determines only the form of resonances of layer modes. For obtaining qualitative conformity of these resonances to the experiment it is necessary to set n ~ 0.2 Pa-s. Thus, A-model represents the radiation losses in the gelatinous layer in the frequency range behind the resonances as mainly "elastic". Because the model reproduces the change of losses with the change of the indentor diameter (Fig. 3), such representation is thought to be close to reality.
PA-model and GA-model give worse conformity to experiments even for one diameter of the indentor (Fig. 2). The PA-model with small viscosity gives the radiation losses with underestimating and it is necessary to increase the value of viscosity r essentially for losses level reproduction. It worsens reproduction of the form of resonances of layer modes and provides the description of radiation losses in the gelatinous layer as the sum of "elastic" and "viscous" components comparable to each other. The change of losses in the model after
changing the indentor diameter, however, does not correspond to the experiment and for reproduction of these losses level appropriate to the new indentor diameter it is necessary to set a new viscosity. The GA-model with small viscosity gives the radiation losses correctly "on the average", but smooth variations around this average level are reproduced here. By means of viscosity increase it is possible to damp the variations on the curves corresponding to the indentors d = 6 mm and d = 10 mm, but they remain essential on the curve corresponding to the indentor d = 16 mm.
Conclusions
Thus, when describing the impedance characteristics of the homogeneous layer within the framework of the models with the pressure source of vibrations, the model with uniform distribution of pressure under the indentor appears to be most adequate to the experiment (except for the description of behavior of stiffness ReK at high frequencies). It is impossible to improve conformity by means of acceptance of "parabolic" or "hyperbolic" pressure profiles. On the basis of these findings we may recommend for description of impedance characteristics of biological tissues to use primary models with uniform distribution of pressure under the indentor.
References
1. Mechanical impedance. In: Shock and vibration handbook, 2nd Ed, Harris C.M. and Crede C.E. (Editors), McGraw-Hill Book Company, New York, 10, 1976.
2. FRANKE E. Mechanical impedance of the surface of the human body. J Appl Physiol, 3(1): 582590, 1951.
3. GOLIKOV V.A., ODINTSOV S.G. Influence of primary transducer on the human body mechanical impedance. Novosti Meditsinskoi Tekhniki, 3: 25-28, 1978 (in Russian).
4. VERMARIEN H., van VOLLENHOVEN E. The recording heart vibrations: a problem of vibration measurement on soft tissue. Med Biol Eng Comp, 22: 168-178, 1984.
5. MRIDHA M., ODMAN S. Characterization of subcutaneous edema by mechanical impedance measurements. J Invest Dermatol, 85(6): 575-578, 1985.
6. TIMANIN E.M. Prospects for measuring the rheologic characteristics of human soft tissues based on the recording of their transverse rigidity. Biofizika, 34(3): 512-516, 1989 (in Russian).
7. TIMANIN E.M. On contribution of shear waves into a transverse stiffness of soft biological tissues in vibrating indentor investigations. 13 International Congress on Acoustics, Belgrade, 4: 215-218, 1989.
8. TIMANIN E.M. A model of formation of impedance properties of soft biological tissues. In: Methods of vibrational diagnostics of rheological properties of soft materials and biological tissue, IAP RAS, Gorky, 75-91, 1989 (in Russian).
9. OKA H. Estimation of muscle fatigue by using EMG and muscle stiffness. Conference Proceedings of 1996 IEEE/EMB, 131-132, 1996.
10. OKA H., IRIE T. . Mechanical impedance of layered tissue. Med Prog Technol, 21(Suppl): 1-4, 1997.
11. SKOVORODA A.R., AGLYAMOV S.R. The reconstruction of mechanical properties of layered viscoelastic media based on impedance measurements. Biofizika, 43(2): 348-352, 1998 (in Russian).
12. TIMANIN E.M. Possibilities of myotonography. Meditsinskaia Tekhnika, 2: 39-41, 1998 (in Russian).
13. TIMANIN E.M., REYMAN A.M., EREMIN E.V. Software-hardware complexes for studying the impedance characteristics of biological soft tissues. The 2-nd Congress of Biophysicists of Russia, Moscow, 2: 626-627, 1999 (in Russian).
14. TIMANIN E.M. On the possibilities of impedance characteristics description of biological soft tissues in models with a pressure source of vibration, Preprint IAP RAS N 488, Nizhny Novgorod, 1999 (in Russian).
15. OESTRAEICHER H. Field and impedance of oscillating sphere in a viscoelastic medium with an application to biophysics. JASA, 23(6): 707-714, 1951.
16. ARVIN G.I. Matching of acoustical transducers with biological objects. Meditsinskaia Tekhnika, 3: 26-29, 1972 (in Russian).
17. GLUSHKOV E.V., GLUSHKOVA N.V., TIMANIN E.M. Impedance and waveguide properties of organic tissues. Akusticheskii Zhurnal, 39(6): 1043-1049, 1993 (in Russian).
18. NOVACKI W. Theory of elasticity, Mir, Moscow, 705, 1975 (in Russian).
19. GLUSHKOV E.V., GLUSHKOVA N.V. On the dynamic contact stiffness for an elastic layer. Prikladnaia Matematika i Mekhanika, 54(3): 474-479, 1990 (in Russian).
20. KLOCHKOV B.N., SOKOLOV A.V. Waves in a layer of soft tissue overlying a hard-tissue halfspace. Akusticheskii Zhurnal, 40(2): 270-274, 1994 (in Russian).
21. OGIBALOV P.M., LOMAKIN V.A., KISHKIN B.P. Mechanics of polymers, Moscow State University, Moscow, 158, 1975 (in Russian).
22. FUNG Y.C. Biomechanics. Mechanical properties of living tissues, New York - Heidelberg -Berlin, Springer - Verlag, 41, 1981.
МЕХАНИЧЕСКИЙ ИМПЕДАНС БИОЛОГИЧЕСКИХ МЯГКИХ ТКАНЕЙ:
ВОЗМОЖНЫЕ МОДЕЛИ
Е.М. Тиманин, Е.В. Ерёмин (Нижний Новгород, Россия)
Разработка математических моделей импедансных характеристик биологических мягких тканей кроме очевидного чисто научного интереса является актуальной в связи с развитием в последнее время способа непрерывного мониторинга механических параметров тканей с высоким временным разрешением по данным одночастотных импедансных измерений и способа реконструкции механических параметров слоистых тканей по данным спектральных импедансных измерений (то есть по частотным зависимостям импедансных характеристик). В данной работе для интерпретации импедансных характеристик биологических мягких тканей и их фантомов предлагается использовать "модели с силовым источником колебаний", основанные на приближениях, принимаемых при решении задачи Лэмба. Возможности этих моделей изучаются на примере сопоставления расчетов в однослойных моделях такого типа с данными экспериментов на однородном слое желатина, полученными средствами специализированного программно-аппаратного комплекса. Проведённый анализ позволяет заключить, что наиболее адекватной экспериментам оказывается модель с равномерным распределением давления под штампом (за исключением описания поведения действительной части комплексной жесткости на высоких частотах). Улучшить соответствие за счет принятия "параболического" или "гиперболического" профиля давления не удаётся. Библ. 22.
Ключевые слова: биологические мягкие ткани, механический импеданс, модели слоистых систем, компьютерные средства измерения
Received 30 November 1999
