K. Miller, K. Chinzei, "Modelling of Brain Tissue Mechanical Properties: Bi-phasic versus Single-phase Approach," in Proc. 3rd Intl Symp Comput Methods in Biomech & Biomed Eng, 1997. MODELLING OF BRAIN TISSUE MECHANICAL PROPERTIES: BI-PHASIC VERSUS SINGLE-PHASE APPROACH K. Miller1 and K. Chinzei2 1. ABSTRACT Recent developments in Robot(Aided Surgery ( in particular, the emergence of automatic surgical tools and robots ( as well as advances in Virtual Reality techniques, call for closer examination of the mechanical properties of brain tissue. The ultimate goal of our research is development of corresponding, realistic mathematical models. The paper discusses two candidates for tissue models: standard, non(linear, biphasic and single(phase, non(linear, viscoelastic. The mechanical behavior of brain tissue is highly non-linear. The stress(strain curves are concave upward containing no linear portion from which a meaningful elastic modulus might be determined. The tissue response stiffens as the loading speed increases, indicating a strong stress(strain rate dependence. The standard methods of modeling tissue as a biphasic continuum face serious problems: strong stress - strain-rate dependence can not be easily explained. According to our experiments, for brain tissue the stresses under fast loading can be six times higher than those under slow loading. Therefore, the use of the single(phase model is recommended. The non-linear, viscoelastic model, based on strain energy function with time dependent coefficients has been developed. The material constants for the brain tissue have been evaluated. Agreement between the proposed theoretical model and experiment is good for compression levels reaching 30% and for loading velocities varying over five orders of magnitude. One advantage of the proposed constitutive model is that it is not difficult to be employed in largere scale finite element computations. 2. INTRODUCTION Mechanical properties of living tissues form a central subject in Biomechanics. Keywords: Brain tissue, Mechanical Properties, Mathematical modelling In particular, the properties of the musculo-skeletal system, skin, lungs, blood and blood vessels have attracted much attention, for example, see [1, 2, 3, 4] and references cited therein. The properties of "very" soft tissues, which do not bear mechanical loads (such as brain, liver, kidney, etc.), have not been so thoroughly investigated.However, recent developments in robotics technology, in particular the emergence of automatic surgical tools and robots [5] as well as advances in virtual reality techniques [6], call for closer examination of the mechanical properties of these tissues. State of the art intra-operative imaging techniques can supply rich information of tissue deformation. Unfortunately, most medical imaging methods need tens of seconds to even tens of minutes to obtain a set of 3D images. Thus, by thinking of using an image measurement in the robot control system, one also has to consider a delay. The possible solution to this challenge is the prediction of the deformation based on the model [7, 8]. To improve the capabilities of surgical operation planning and surgeon training systems based on the virtual reality techniques ([6] and references cited therein), force feedback is needed. To achieve the above goals, the appropriate "very" soft tissue models are required. Knowledge of the mechanical properties of soft tissues and ultimately of their mathematical models is also required for registration. It encompasses matching images of different modality, such as MRI and Single Photon Emission Computed Tomography (SPECT), defining relations between coordinate systems (e.g., between a coordinate system associated with imaging equipment and those of robotic tools in an operating room), segmentation of reference features and defining disparity or similarity functions between extracted features [9]. Registration is a key technique for the computer-integrated surgery. Registration procedures involving rigid tissues are now well-established. If rigidity is assumed, it is sufficient to find several points such that their mappings between two coordinate systems are known. Registration of soft tissues is much more difficult because it requires a knowledge about local deformations. Here comes the place for accurate models of tissue deformation behavior. The reported experimental data on the mechanical properties of brain tissue are limited. Ommaya in [10] described mammalian brain as a "soft, yielding structure, not as stiff as a gel nor as plastic as a paste". Autors of [11] and [12] tried to establish elastic parameters of brain tissues by measuring induced changes in intra-cranial pressure. The experimental data, which might be used to determine constitutive relations for brain tissue, was reported in [13] and [14]. From the cited papers, it is difficult to extract precise information about the conditions under which the experiments were performed, especially concerning sample sizes and prevention of the adhesion between platens and sample surfaces. Based on these experimental results, nonlinear constitutive relations for human brain tissue were proposed [15, 16]. However both the experimental results and theoretical investigations concentrated on rapid loading conditions resulting in large strain-rates, typical for accidents and injury modelling, and can not find the direct application in surgical procedure simulations. 2. UNCONFINED COMPRESSION EXPERIMENT OF SWINE BRAIN TISSUE We decided to perform unconfined uniaxial compression experiments. An alternative was a confined compression experiment, used for example by Mow and coworkers [3] to validate biphasic models of soft cartilage tissues. The more detailed description of the in-vitro unconfined compression of fresh swine brain can be found in [17, 18]. Here we would like to summarise the results only. Cylindrical samples of diameter ~30 mm and height ~13 mm were cut. Four samples were taken from the frontal and posterial portions of the Sylvian fissure of each hemisphere for each swine brain. The ventricle surface and the arachnoid membrane formed the top and bottom faces of the sample cylinder. Thus the arachnoid membrane and the structure of the sulci remained as parts of each specimen. Uniaxial unconfined compression of swine brain tissue was performed in a testing stand, shown in Fig. 1. a) b) Fig. 1: Experiment setup a) general view with components: 1 - specimen and the loading platens, 2 - load cell to measure axial force, 3 - micrometer to measure axial displacement , 4 - laser to measure radial displacement. b) layout with coordinate axes The main testing apparatus was a UTM-10T (Orientec Co.) tensile stress machine. The vertical displacement (along z-axis in Fig. 1b) was measured by a micrometer with analog output. The experiment was documented by automatically taking CCD camera images. In the paper we will discuss the results obtained during the loading phase, for three loading velocities: · fast: 500 mm min-1 ¾ the fastest loading speed possible with our equipment), corresponding to the strain rate of about 0.64 s-1, · medium: 5 mm min-1, corresponding to the strain rate of about 0.64´10-2 s-1, and · slow: 0.005 mm min-1, corresponding to the strain rate of about 0.64´10-5 s-1. We performed 12 fast, 13 medium speed and 6 slow tests. The number of slow tests was limited because after each tissue delivery (usually 2 brains) we could perform only one overnight test. Figure 2 shows the Lagrange stress3 versus true strain (e = ln lz , where lz is a stretch in vertical direction, Fig. 1b) curves for three loading velocities. The standard deviation of the measurements and the theoretical predictions are indicated. The stress-strain curves are concave upward for all compression rates containing no linear portion from which a meaningful elastic modulus might be determined. The tissue response stiffened when the loading speed increased, indicating a strong stress-strain rate dependence. The results shown in Fig. 2 are in general agreement with those published in (Estes and McElhaney, 1970). a) b) c) Fig. 2. Lagrange stress - true strain relations for swine brain tissue, experimental and theoretical results. Error bars indicate standard deviation. Loading speeds: a) 500 mm min-1, corresponding to the stretch rate of about 0.64 s-1, b) 5 mm min-1, corresponding to the stretch rate of about 0.64´10-2 s-1, c) 0.005 mm min-1, corresponding to the stretch rate of about 0.64´10-5 s-1. 3. CANDIDATES FOR BRAIN TISSUE CONSTITUTIVE MODELS 3.1. Brain as a Biphasic Continuum The concept that the soft tissues can be treated as biphasic continua consisting of solid deformable porous matrix and penetrating fluid is widespread. It has been particularly useful in cartilage biomechanics, see for example [3]. The linear biphasic model of brain was proposed, for example, in [19, 20]. The governing equations of biphasic continuum are: Continuity: (1) Equilibrium: (2) where: fa - a phase content (f - fluid phase, s - solid phase), va - velocity of a phase. sa - a phase Cauchy stress tensor, Pa - diffusive momentum exchange between phases. When writing the equilibrium equation we neglected inertial body forces. If we make the additional assumption that the fluid is inviscid and that the diffusive momentum exchange is proportional to the relative velocity between phases, the constitutive equations are as follows: (3) (4) (5) where: p - apparent pressure, K - diffusive drag coefficient function, - Cauchy stress tensor of the solid phase. I - identity tensor (rank three), In general, the diffusive drag coefficient K is not constant. It is usually considered to be dependent (exponentially) on strain. The accepted way to relate the stresses to the deformation (in a solid matrix) is by means of the Helmholtz free energy: (6) where: S - Second Piola-Kirchoff stress of the solid phase (measured with respect to the original configuration) W - the Helmholtz free energy (per unit volume) function of the solid phase. This function depends on the current deformation only. E - Green's strain tensor (relative to original configuration) of the solid phase. This approach is based on the assumption that the solid phase stress depends only on the current deformation. Therefore, there is no energy dissipation in the solid, but the dissipation comes from interactions between phases only. It was shown analytically [21] for the linear biphasic model and confirmed numerically for the non-linear case [7, 8] that the ratio of the instantaneous stress (after sudden movement of the upper platen) to the equilibrium stress (after sufficiently long time following load application), as predicted by the biphasic theory, cannot be larger than , where n is the Poisson's ratio of the solid phase. This poses a severe limit on the stress dependence on loading velocity. The biphasic theory in its present form cannot be accepted for very soft tissues (e.g. brain), for which the stresses for the largest loading velocity in our experiments are about six times higher than for the smallest one (Fig. 2). 3.2. Brain as an Inelastic Single-Phase Continuum The reason for standard biphasic theory's inability to discribe strong stress-strain rate dependence is the underlying assumtion of solid matrix hyperelasticity. In the case of the confined compression experiment the relative velocity of the phases is almost equal to the velocity of the solid phase. Therefore the dissipation may be accounted for by choosing suitably large darg coefficient K, equation 5. In the unconfined experiment the velocity of solid phase is much larger than the relative velocity between phases. The results show that the solid phase is inherently dissipative. The dissipation in the system cannot be accounted for by adjusting the drag coefficient. The simple, phenomenological, modelling method introducing suitable dissipation into the system at the expense of single-phase description, is discussed below. 3.2.1. Modeling of finite deformation non-linear tissue behavior Let's start with the modeling of non-linear stress strain dependence using the strain energy function of the following form: (7) Where the strain invariants are: ; ; , (8) B is a left Cauchy-Green strain tensor. For infinitesimal strain conditions, the sum of constants C10 and C01 have a physical meaning of one half of the shear modulus: . (9) The energy dependence on strain invariants only comes from the assumption that brain tissue is initially isotropic. The assumption of tissue incompressibility results in setting the third strain invariant equal to one. The first two terms in (7) form a well known Mooney-Rivlin energy function, originally developed for incompressible rubbers (for discussion see [22]). In our experiment the deformation was orthogonal, and the left Cauchy-Green strain tensor had only diagonal components: , where is a stretch in vertical direction (Fig. 1b) (10) In such a situation, taking and , the only non-zero Lagrange stress components can be computed from the simple formula: . (11) 3.2.2. Modeling of the loading velocity dependent tissue behavior. To model the time dependent behavior of the tissue we propose to write the coefficients in the formula for energy function (7) in the form of a exponential series , (12) and the energy function in the form of convolution integral (13) From (11) we obtain: (14) Equation (14) served as a basis for the comparison of the theory and experiment. 3.2.3. Determination of material constants for swine brain tissue In the case of the compression with constant velocity, the integral (14) can be evaluated analytically4 . The result obtained with Mathematica [23] is long and is not presented here. It is clear from equation 14 that the expression for Lagrange stress is linear in material parameters Cij . Therefore it was easy to find them using least square method. To model accurately the tissue behavior for a wide range of loading velocities, we found it necessary to use two time-dependent terms in the Cij expansion (12) and to include second oredr terms in energy function (N=2 in eq. 13). A good agreement with experiment for all three loading velocities (Fig. 3) has been obtained for the following values of the parameters: time constants in eq. 12: t1=0.5 [s]; t2=50. [s]; equilibrium parameters: C10,inf = C01,inf =62.50 [Pa]; C11 = 0; C20,inf = 0.; C02,inf = 0.; parameters multiplying exponential with characteristic time t1: C101 = C011 = 39.66 [Pa]; C201 = C021 = 869.82 [Pa]; parameters multiplying exponential with characteristic time t2: C102 = C012 = 366.86 [Pa]; C202 = C022 = 310.65 [Pa]; 4. DISCUSSION AND CONCLUSIONS This study discusses two distinct mathematical model of brain tissue mechanical properties. The strong stress-strain rate dependence prohibits the use of standard biphasic models for brain tissue modeling. Therefore, the use of the single-phase, non-linear, viscoelastic model based on the concept of the strain energy function, written in the form of convolution integral with coefficient expressed in the form of exponential series, is advocated. One advantage of the proposed model is that the constitutive equation developed here is already available in ABAQUS [24] and can be used immediately for larger scale FEM computations. How to use the in vitro experimental results in the more realistic in vivo environment remains an open question. Further research is needed to determine brain tissue constitutive models, which would incorporate the influence of the blood and cerebrospinal fluid pressure and flow. ACKNOWLEDGMENTS The authors would like to thank most sincerely all the members of Biomechanics Division of Mechanical Engineering Laboratory, MITI and Dr. Tetsuya Tateishi not forgetting Dr. Yoshio Shirasaki's and Mr. Kazuya Machida, for their indispensable help. The financial support of The New Energy and Industrial Technology Development Organization (NEDO) and Department of Mechanical and Materials Engineering of The University of Western Australia is gratefully acknowledged. REFERENCES Fung Y.C., Biomechanics. Mechanical properties of Living Tissues. Springer-Verlag, New York, 1981. Gallagher, R.H., Simon, B.R., Johnson, P.C. and Gross, J,F. eds., Finite Elements in Biomechanics. John Wiley & Sons. New York, 1982. Mow, V.C., Ateshian, G.A. and Spilker R.L., Biomechanics of Diarthrodial Joints: A Review of Twenty Years of Progress. Trans. ASME, J. Biomech. Eng., 1993, 115, 460-467. Schmid-Schonbein, G.W., Woo, S.L-Y. and Zweifach, B.W. eds., Frontiers in Biomechanics. Springer, New York, 1986. Brett, P.N., Fraser, C.A., Henningan, M., Griffiths, M.V. and Kamel Y., Automatic Surgical Tools for Penetrating Flexible Tissues. IEEE Eng. Med. Biol., 1995, pp.264-270. Burdea, G., Force and Touch Feedback for Virtual Reality. Wiley. New York, 1996. Miller, K. and Chinzei, K., Modeling of Soft Tissues, Mechanical Engineering Laboratory News, 1995, 12, 5-7 (in Japanese). Miller, K. and Chinzei, K., Modeling of Soft Tissues Deformation, Journal Computer Aided Surgery, 1995, 1, Supl., Proc. of Second International Symposium on Computer Aided Surgery, Tokyo Women's Medical College, Shinjuku, Tokyo, 62-63. Lavallˇe, S., Registration for Computer Integrated Surgery: Methodology, State of the Art. Computer-Integrated Surgery, 1995, MIT Press, Cambridge Massachusetts, pp. 77-97. Ommaya, A.K., Mechanical Properties of Tissues of the Nervous System. J. Biomech., 1968, 1, 127-138. Walsh, E.K. and Schettini, A., Calculation of brain elastic parameters in vivo. Am. J. Physiol. 1984, 247, R637-R700. Sahay, K.B., Mehrotra, R., Sachdeva, U. and Banerji, A.K., Elastomechanical Characterization of Brain Tissues. J. Biomechanics. 1992, 25, 319-326. Estes, M.S. and McElhaney J.H., Response of Brain Tissue of Compressive Loading, ASME Paper No. 70-BHF-13, 1970. Galford, J.E. and McElhaney, J.H., A Viscoelastic Study of Scalp, Brain and Dura. J. Biomech., 1970, 3, 211-221. Pamidi, M.R. and Advani, S.H., Nonlinear Constitutive Relations for Human Brain Tissue. Trans. ASME, J. Biomech. Eng., 1978, 100, 44-48. Mendis, K.K., Stalnaker, R.L. and Advani S.H. (1995) A Constitutive Relationship for Large Deformation Finite Element Modeling of Brain Tissue. Trans. ASME, J. Biomech. Eng. 117, 279-285. Chinzei K., Miller K., "Experimental evaluation of compressive behavior of very soft tissues," Mechanical Engineering Laboratory News (ISSN 0286-2271), No. 11, 1996, pp. 1-3 (in Japanese). Miller, K. and Chinzei, K., Constitutive Modeling of Brain Tissue; Experiment and Theory, Journal of Biomechanics (submitted). Nagashima, T., Tamaki, N., Matsumoto, S., Horwitz, B. and Seguchi, Y. (1987) Biomechanics of Hydrocephalus: A New Theoretical Model. Neurosurgery. 21, No. 6, 898-903. Basser, P.J (1992) Interstitial Pressure, Volume, and Flow during Infusion into Brain Tissue. Microvasc. Res. 44, 143-165. Armstrong, C.G., Lai, W.M. and Mow, V.C. (1984) An Analysis of the Unconfined Compression of Articular Cartilage. Trans. ASME, J. Biomech. Eng. 106, 165-173. Rivlin, R.S. (1984) Forty Years of Nonlinear Continuum Mechanics. In Proceedings of the IX Int. Congress on Rheology, Mexico, pp. 1-29. Wolfram Res., Inc. (1994) Mathematica, Version 2.2, Wolfram Research Inc., Champaign, Illinois, USA. ABAQUS Theory Manual (1992) Version 5.2, Hibbit, Karlsson & Sorensen, Inc. 1 Lecturer, Department of Mechanical and Materials Engineering, University of Western Australia, Nedlands/Perth, WA 6907, Australia, email: kmiller@mech.uwa.edu.au 2 Researcher, Biomechanics Division, Mechanical Engineering Laboratory, MITI, Namiki 1-2, Tsukuba, Ibaraki 305, Japan, email: chin@mel.go.jp 3 To calculate Cauchy (true) stresses, the precise measurement of the cross-section area during loading or the assumption of tissue incompressibility is needed. 4 The result contains exponential integrals.