Human temporomandibular joint disc cartilage as a poroelastic material

Clinical Biomechanics 18 (2003) 69–76 www.elsevier.com/locate/clinbiomech

Human temporomandibular joint disc cartilage as a poroelastic material M. Beek, J.H. Koolstra*, T.M.G.J. van Eijden Department of Functional Anatomy, Academic Center for Dentistry Amsterdam (ACTA), Meibergdreef 15, 1105 AZ Amsterdam, The Netherlands Received 18 October 2001; accepted 10 September 2002

Abstract Objective. The hypothesis was tested that a poroelastic material model is potentially able to describe the mechanical behavior of cartilaginous tissues in dynamic indentation experiments. Design. This hypothesis was tested by comparing the results from model predictions with results obtained from cyclic indentation experiments. Background. The characteristics of cartilaginous tissues in general and of the temporomandibular joint disc in particular are generally identified by static confined or unconfined indentation experiments, while under physiologic circumstances these tissues are mostly loaded dynamically. Methods. Dynamic indentation experiments were simulated using an axisymmetric finite element model. The results from the simulations were qualitatively compared with the experiments. Results. The simulations showed several similarities with the experiments when the solid matrix was assumed to be hyperelastic. Both the maximum stress and the amount of energy dissipated decreased in each subsequent cycle. Furthermore, a similar dependency on the indentation frequency and amplitude was found. Conclusions. This qualitative study showed that a poroelastic material model can describe the dynamic behavior of the temporomandibular joint disc, provided that the solid matrix is modeled as hyperelastic. Relevance Temporomandibular disorders are presumably related to joint load distributions. Besides large static, dynamic loads are considered as a risk factor for cartilaginous wear. Dynamical loads, however, are also considered to stimulate the biosynthetic activity of cartilaginous tissues. Biomechanical analysis can be applied to estimate nonmeasurable joint loads. This enables to understand the underlying mechanisms of temporomandibular disorders, necessary to develop methods to prevent, diagnose and cure joint disorders. The present study shows that a poroelastic material model can be applied successfully to model the dynamical behavior of the temporomandibular joint disc in such analyses. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Temporomandibular joint disc; Poroelastic; Indentation experiment; Finite element method

1. Introduction Clinical problems concerning the temporomandibular joints are quite common among western populations [1]. The mechanical loads acting on the joints during habitual functions like chewing and talking, and during dysfunction like bruxism (grinding), are assumed to play

*

Corresponding author. E-mail address: [email protected] (J.H. Koolstra).

a role in the development of these problems. Under the influence of loads and deformations, the cartilaginous intra-articular disc of the joints may wear and become thinner or even perforated. This may lead to wear of the articular surfaces, which hampers the movements of the lower jaw. Assessment of mechanical loads in the joint, therefore, is crucial to understand its normal and pathological behavior. The mechanical behavior of the disc in the joint has been investigated by means of finite element modeling [2–7]. The results of these studies suggested that the disc plays an important role in distributing and absorbing

0268-0033/03/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 2 6 8 - 0 0 3 3 ( 0 2 ) 0 0 1 3 5 - 3

70

M. Beek et al. / Clinical Biomechanics 18 (2003) 69–76

loads acting on the joint. The applied models, however, were (quasi-)static and, therefore, the applicability of their results is limited to situations where the jaw hardly moves (e.g., clenching). The mechanical behavior of the temporomandibular joint disc has also been investigated in vitro by means of experimental studies [8–12]. The results of these studies indicate that its mechanical behavior is nonlinear, anisotropic and time-dependent. These studies were of static or quasi-static nature. However, the behavior of cartilaginous tissues during habitual dynamic loading differs from this (quasi-)static behavior [13]. Recently, sinusoidal indentation experiments were performed to investigate the behavior of the disc under more physiologic loading conditions [14]. It was found that the maximum reaction force and the amount of dissipated energy decreased in time (Fig. 1). Furthermore, these characteristics appeared to be dependent on location, indentation amplitude, and excitation frequency. This specific behavior can probably be attributed to mechanical interactions between the solid matrix (mainly collagen and proteoglycans) and the interstitial fluid (mainly water) of its cartilage [15]. In order to incorporate these interactions in biomechanical models, a biphasic theory has been developed [16], which assumed that a large part of the loads acting on cartilaginous structures is carried by interstitial fluid pressurization [17]. Furthermore, it has been demonstrated that a poroelastic material model is capable of modeling the behavior of cartilaginous structures similar to this biphasic model [18]. In contrast to the biphasic theory, the poroelastic theory is implemented in the majority of commercially available finite element software packages. Consequently, this material model seems to be advantageous for application in dynamic finite element models involving cartilage in general and the temporomandibular joint disc in particular. Although this material model has been used in (quasi-)static simulations of confined and unconfined compression experiments [18–20], it has not been investigated for its ability to adequately describe the specific characteristics found in the recent experiments, yet. The purpose of the present study was to test the hypothesis that the poroelastic material model is applicable to describe the dynamic behavior of the cartilaginous disc of the human temporomandibular joint. Therefore, a finite element model was developed and subjected to cyclic dynamical indentations similar to the experiments with regard to human temporomandibular joint discs as described by Beek et al. [14]. The focus of the simulations was directed to the ability to adequately model the specific qualitative characteristics found in the experiments, namely the decrease in time of both the maximum force and amount of dissipation, and the dependency on both the strain amplitude and the excitation frequency.

Fig. 1. (A) Experimental set-up of the dynamical indentation experiments with human temporomandibular joint discs performed by Beek et al. [14]. During the experiments, an intact disc was submerged in saline that was kept at a temperature of about 37 °C. (B) Typical time series of the measurements at an indentation frequency of 0.05 Hz. On top, the sinusoidal displacement of the top indenter. Below, the reaction force of the disc. (C) Typical stress–strain curves, derived from the measurements in (B) for 10 subsequent loading cycles.

2. Methods 2.1. Model An axisymmetric finite element model of cartilaginous tissue (radius rcart ¼ 3:49 mm) was created to simulate dynamical indentation experiments performed on fresh human temporomandibular joint discs (Fig. 1, [14]). This model was loaded by undeformable and im-

M. Beek et al. / Clinical Biomechanics 18 (2003) 69–76

71

ters was assumed to be zero. This allowed for a free fluid flow across these boundaries in an aqueous environment. At the symmetry axis the radial displacement of the solid matrix and the radial velocity of the fluid were assumed to be zero. It was assumed that the disc was perfectly glued to both indenters. The simulations were performed using MARC 7.3.2 (MSC.Software, Los Angeles, USA) on a Cray C916 computer. 2.2. Simulations

Fig. 2. Axisymmetric finite element model for the dynamical indentation simulations. On the symmetric axis, the symmetric boundary conditions are shown. The displacements of the nodes in contact with the bottom indenter were restrained. Indentation was simulated by prescribing a displacement to the nodes in contact with the top indenter.

permeable indenters (radius rind ¼ 1:97 mm). The disc model was larger than the indenters, because the experiments had been performed with intact discs. The mesh consisted of 160 isoparametric quadrilateral elements (Fig. 2). The cartilaginous tissue was supposed to consist of a solid matrix and interstitial fluid, and was modeled using the poroelastic theory. The formulation of poroelasticity used a linear elastic relation between the Green–Lagrange strain e and the second Piola Kirchhoff stress r as the constitutive relation for the solid matrix. Initially, a YoungÕs modulus and a PoissonÕs ratio of the isotropic solid matrix of 0.4667 MPa and 0.1667, respectively, were applied [18]. The results from the experiments [14] showed that the behavior of the temporomandibular joint disc was nonlinearly dependent on the strain. Therefore, hyperelasticity of the solid matrix was assumed, which was modeled by a third-order polynomial relation between the Green–Lagrange strain e and the second Piola Kirchhoff stress r: r ¼ c3e3eq þ c2e2eq þ c1eeq þ c0

ð1Þ

in which eeq ¼ sqrtðe211 þ e222 þ e233 Þ

ð2Þ

This hyperelasticity was implemented in the material model in series with the linear elastic model for the solid phase through a user subroutine. The applied values for c3, c2, c1, and c0 were approximated to 103 , 103 , 0, and 1, respectively, from the experimentally obtained stress– strain curves. The permeability of the solid matrix was 7:5  1015 m4 /N s [18]. Additionally, the porosity of cartilage was assumed to be 0.7, according to a 70% fluid content of cartilaginous structures [21,22]. The pore pressure at the boundary of the disc outside the inden-

In our simulations, the temporomandibular joint discs were loaded by applying a sinusoidal displacement (Dh) in axial direction to the nodes in contact with the top indenter (Fig. 2). The strain value e as measured in the experiments [14] was calculated as the quotient of the axial displacement of the node on the symmetry axis and the unloaded thickness (h0) of the cartilage between the indenters (e ¼ Dh=h0). The nodes in contact with the bottom indenter were fixed. The reaction forces gathered at the node on the symmetry axis by means of kinematic tyings, enabled to instantaneously obtain the total reaction force F . The applied stress (r ¼ F =A) was defined 2 as the reaction force divided by the area (A ¼ prind ) of the indenters. After ten sinusoidal indentation cycles, relaxation was allowed for 200 s to verify the modelÕs return to its original state. The ten indentation cycles were simulated in 1500 increments and the following 200 s of relaxation in 500 increments. At present, there is no agreement on the YoungÕs modulus of the temporomandibular joint disc [8–10,12]. Furthermore, reliable quantitative data concerning the permeability of cartilage and in particular of the temporomandibular joint disc is scarce [18,22,23]. Therefore, the sensitivity of the porohyperelastic model to changes in the YoungÕs modulus and the permeability was evaluated by performing subsequent simulations in which these parameters were varied. For these parameters we applied values that were ten times smaller and ten times larger compared to the values applied by Prendergast et al. [18]. The characteristics of the porohyperelastic material model were further investigated by determining the influence of the indentation amplitude and the indentation frequency, respectively. According to the experimental data [14], we varied the strain amplitude at a frequency of 0.05 Hz. The amplitudes applied were 20%, 30% and 40%. Additionally, we varied the strain frequency at an amplitude of 30%. The frequencies applied were 0.02, 0.05, and 0.1 Hz.

3. Results Duplicating the simulations described in the article by Prendergast et al. [18], produced the same results. Using

72

M. Beek et al. / Clinical Biomechanics 18 (2003) 69–76

between the loading and unloading curve for each subsequent indentation was less after the second loading cycle. Like in the experiments [14], the differences in shape of each subsequent stress–strain curve gradually vanished during the series of indentations. The final shape represents a dynamic equilibrium in which fluid flow during loading and unloading is equal but opposite. 3.1. Sensitivity analyses The results of the sensitivity analyses are shown in Fig. 4. A ten times larger elasticity resulted after ten subsequent cycles in a maximum force, which was ten times larger compared to the reference simulation. This larger maximum force led to a larger area enclosed by the curve and thus in more energy dissipation. A decrease in elasticity of the solid matrix had a reciprocal

Fig. 3. (A) Simulation of cyclic dynamic indentation experiments using a linear poroelastic model. Indentation frequency and amplitude were 0.05 Hz and 30%, respectively. (B) Simulation of dynamical indentation experiments using a poroelastic model with a nonlinear solid matrix, according to r ¼ c3e3eq þ c2e2eq þ c1eeq þ c0, with c3 ¼ 103 , c2 ¼ 103 , c1 ¼ 0, and c0 ¼ 1. Indentation frequency and amplitude were 0.05 Hz and 30%, respectively.

a linear poroelastic model [18], simulation of the cyclic indentation experiments (ten cycles) resulted in an almost linear viscous behavior (Fig. 3A). In contrast to viscoelasticity, poroelasticity predicted the development of an increasing underpressure in the tissue during the unloading phase of each subsequent indentation. The reaction force at 30% indentation was 4 N during the first cycle. The strains in the solid matrix were inhomogeneously distributed in axial direction. After simulating a relaxation of 200 s all mechanical quantities had returned to their original state. Introducing hyperelasticity of the solid matrix according to Eqs. (1) and (2), resulted in superproportional stress–strain curves found experimentally (compare Fig. 1C and Fig. 3B). This means that the observed property of cartilage to become stiffer at higher deformations due to a decrease in permeability of the solid matrix and a subsequent pressurization of interstitial fluid, can be modeled by a porohyperelastic model. The decrease of the difference

Fig. 4. (A) Influence of the YoungÕs modulus E on the stress–strain curves obtained from simulations of ten subsequent indentation cycles. The value of E was varied with respect to Eref ¼ 0:4667 MPa. Indentation frequency and amplitude were 0.05 Hz and 30%, respectively. (B) Influence of the permeability k on the stress–strain curves obtained from simulations of ten subsequent indentation cycles. The value of k was varied with respect to kref ¼ 7:5  1015 m4 /N s. Indentation frequency and amplitude were 0.05 Hz and 30%, respectively. For visibility reasons, only the first cycle is displayed.

M. Beek et al. / Clinical Biomechanics 18 (2003) 69–76

73

effect. An increase in the resistance of the solid phase against deformations will thus increase the load bearing capabilities of the cartilage. A ten-fold increase of the permeability, resulted in a decrease of the resistance to indentation and herewith the maximum reaction force (Fig. 4B). The area enclosed by the stress–strain curve was also decreased. The shape of the stress–strain curve did not change much after ten subsequent cycles. A decrease of the amount of permeability led to an increase in the maximum reaction force in the first compression cycle. In addition, the amount of energy dissipation was larger. While the maximum force at equilibrium was hardly influenced by the value of the permeability, the asymptotical decrease of this maximum force was larger when the permeability was smaller. When fluid flow eases with increasing permeability, the collagen network will be more heavily loaded and the shock absorbing capabilities of the cartilage will decrease. 3.2. Dependency on indentation amplitude and frequency Fig. 5 shows the results of the simulations in which the amplitude of the indentation with a frequency of 0.05 Hz was varied. The reaction force was superproportionally dependent on the indentation amplitude (Fig. 5A) and the shape of the stress–strain curve was hardly influenced by the amplitude of the indentation (Fig. 5B). Unfortunately, the simulations with an indentation amplitude of 40% did not converge after six cycles. Therefore, only the first six cycles are depicted. Fig. 6 shows the results of the simulations in which the frequency of the indentation with an amplitude of about 30% was varied. The maximum reaction force in the first indentation cycle was proportionally dependent on the frequency (Fig. 6A). However, the reaction force after ten cycles (approximation of equilibrium) was larger for an indentation frequency of 0.05 Hz than for a frequency of 0.1 Hz. When the disc was indented at a higher frequency, it behaved stiffer during the loading phase of the indentation (Fig. 6B). The unloading phase was hardly affected by the frequency. Fig. 7 shows the maximum stress obtained at all combinations of indentation amplitude and frequency. This parameter appeared to be more dependent on the amplitude than on the frequency.

4. Discussion In the present study, an axisymmetric finite element model was used to simulate dynamical indentation experiments with human temporomandibular joint discs. Apart from the simulations with an indentation amplitude of 40%, there were no signs of convergence problems or oscillations. Previously, such experiments

Fig. 5. Influence of the indentation amplitude. The values applied for the indentation amplitude were 20%, 30%, and 40%, respectively. The indentation frequency was 0.05 Hz. (A) Time series of the reaction force. (B) Stress–strain curves.

Fig. 6. Influence of the indentation frequency. The values applied for the indentation frequency were 0.02, 0.05, and 0.1 Hz, respectively. The indentation amplitude was 30%. (A) Time series of the reaction force. (B) Stress–strain curves.

showed that the mechanical behavior of the disc was nonlinear and time-dependent [14]. The fluid content in cartilaginous structures has been shown to comprise about 70–85% of the total mass [21,22]. The remaining

74

M. Beek et al. / Clinical Biomechanics 18 (2003) 69–76

Fig. 7. Strain and frequency dependency of the maximum stress in the first cycle, for all combinations of indentation amplitude and frequency.

part mainly consists of a collagen network and proteoglycans. Several studies have shown that this large amount of fluid plays an important role in the complex mechanical behavior of cartilage [17,24]. Therefore, it is crucial to apply a material model that includes both fluid and solid constituents, enabling to distinguish between the mechanical functions of these different constituents. The poroelastic theory has been demonstrated a promising alternative for this purpose. Its quasi-static capabilities have been established earlier [18–20]. Knowledge about its dynamical capabilities is still incomplete. The present study, therefore, was performed to investigate its capability to describe the remarkable characteristics of the behavior of the temporomandibular joint disc as found in dynamical indentation experiments. During the experiments the temporomandibular joint discs were cyclically loaded with frequencies up to 0.1 Hz, as limited by the capabilities of the applied material testing machine [14]. A similar range was applied in the present simulations to be able to compare the results. Nonetheless, the temporomandibular joint disc is loaded at frequencies up to about 1 Hz under physiologic circumstances. It has been shown that the characteristic frequency for cartilage specimen averages between 0.00044 Hz and 0.03 Hz [25,26]. Therefore, it can be assumed that the behavior of the temporomandibular joint disc will not change considerably for loading frequencies higher than 0.1 Hz. Extrapolating the results of our experiments and simulations (Fig. 7), suggests that the disc will become stiffer at higher frequencies. This can be comprehended by a smaller time span available

for fluid flow at higher frequencies and the incompressibility of the fluid [25,27]. In the indentation experiments by Beek et al. [14], intact temporomandibular joint discs were glued between the indenters. Therefore, in the present study, the disc was also modeled beyond the borders of the indenters. In the tissue remote from the volume between the indenters and their near surroundings, the deformations were negligible, ensuring that the results were hardly affected by the limited measures of the model. Free fluid flow across the free boundaries of the disc was allowed. This boundary condition is in agreement with various studies which have shown that fluid flow across the boundaries does indeed occur, for example in confined compression experiments (review: [28]). In the present study, the relatively small amount of temporomandibular joint disc between the indenters was supposed to be homogeneous. However, various studies have shown that the internal organization of cartilaginous structures (e.g., direction and density of collagen fibers) varies with location [29]. Such regional variations in the internal structure have also been reported to exist in the temporomandibular joint disc [30–32]. Furthermore, local increase in elasticity of the disc was found in patients with severe internal derangement [33]. These variations possibly underlie observed regional variations in the mechanical behavior of cartilage [34–37] and the temporomandibular joint disc [8,9,11,14]. The present study enables to model the behavior of a whole disc including heterogeneities by assigning individual material parameters extrapolated from the present sensitivity analyses to relevant regions. To be able to check the poroelastic formulation in MARC, we applied the same values for the material parameters as Prendergast et al. [18]. We are aware that these values may differ from those of human temporomandibular joint disc cartilage during compressive loading. However, at present an unambiguously more reliable complete set of material parameters is still lacking. Furthermore, the stress–strain curves from the indentation experiments [14] were determined for the whole structure, meaning that they did not only involve the stresses in the solid matrix but also the pressurization of the interstitial fluid. Because results from static confined compression experiments (no fluid flow) were not available, the exact mechanical behavior of the solid matrix remained unknown. It must be noted that there is a large biodiversity in biological tissues. Furthermore, the disc is not a homogeneous structure [14]. Therefore, the present results have to be interpreted as qualitative. The performed sensitivity analyses demonstrated that any specific experiment and heterogeneous internal structure can be mimicked by adapting the applied polynomial and its constants. Furthermore, true strains and true stresses need to be incorporated in the poro-

M. Beek et al. / Clinical Biomechanics 18 (2003) 69–76

elastic formulation to increase accuracy when large displacements are involved. In our experiments we detected negative forces during the unloading phase, which have also been described in literature [25]. The present simulations, however, predicted somewhat larger negative forces than found in the experiments. This can be attributed to a better connection between disc and indenters in the simulations than in the experiments. While fluid is expelled out of the cartilage during loading, an underpressure in the cartilage between the indenters develops when the indenters separate during unloading. Therefore, tensile forces are needed to restore the original thickness ðh0Þ of the tissue between the indenters. When the contact between the disc and the indenters is not perfect, which might be the case in the experiments, less tensile forces can be transmitted to the cartilage. The results of the simulations show various similarities with the results obtained from the dynamical indentation experiments. Both the reaction force and the amount of energy dissipation were decreased in subsequent cycles. In the simulations, the maximum stress appeared to be superproportionally dependent on the strain amplitude and only marginally on the frequency (Fig. 7). The same characteristics have also been reported in the experimental study. Consequently, the results indicate that the poroelastic material model adequately describes the strain and frequency dependency of the temporomandibular joint discs as found experimentally.

5. Conclusions The present study enables to implement the specific dynamic characteristics of the temporomandibular joint disc found experimentally [14] into biomechanical models in a relatively easy way. Such enables to perform biomechanical analyses as a tool to develop means for prevention and treatment of temporomandibular disorders. The temporomandibular joint disc has been demonstrated to play an important role in the mechanics of the joint during static tasks [2–7]. As it is assumed that the temporomandibular joints are mostly used dynamically during habitual tasks, dynamic analyses seem to be the most appropriate for this purpose. The results of the present study indicate that the poroelastic material model can be applied to model essential characteristics of the dynamical behavior of the temporomandibular joint disc.

Acknowledgements This work was sponsored by the National Computing Facilities Foundation (NCF) for the use of super-

75

computing facilities, with the financial support from the Netherlands Organization for Scientific Research (NWO). This research was institutionally supported by the Interuniversity Research School of Dentistry, through the Academic Center for Dentistry Amsterdam. We gratefully thank Academic Computing Services Amsterdam (SARA) for technical support. We also thank P. Prendergast and D. Lacroix from Trinity College in Dublin and C. Gelten at MSC. Software for their help. References [1] Carlsson GE. Epidemiology and treatment need for temporomandibular disorders. J Orofac Pain 1999;13:232–7. [2] Chen J, Xu L. A finite element analysis of the human temporomandibular joint. J Biomech Eng 1994;116:401–7. [3] DeVocht JW, Goel VK, Zeitler DL, Lew D. A study of the control of disc movement within the temporomandibular joint using the finite element technique. J Oral Maxillofac Surg 1996;54:1431–7. [4] Chen J, Akyuz U, Xu L, Pidaparti RMV. Stress analysis of the human temporomandibular joint. Med Eng Phys 1998;20:565–72. [5] Nagahara K, Murata S, Nakamura S, Tsuchiya T. Displacement and stress distribution in the temporomandibular joint during clenching. Angle Orthod 1999;69:372–9. [6] Beek M, Koolstra JH, Van Ruijven LJ, Van Eijden TMGJ. Threedimensional finite element analysis of the human temporomandibular joint disc. J Biomech 2000;33:307–16. [7] Beek M, Koolstra JH, Van Ruijven LJ, Van Eijden TMGJ. Three dimensional finite element analysis of the cartilaginous structures in the human temporomandibular joint. J Dent Res 2001;80: 1913–8. [8] Tanne K, Tanaka E, Sakuda M. The elastic modulus of the temporomandibular joint disc from adult dogs. J Dent Res 1991; 70:1545–8. [9] Teng S, Xu Y, Cheng M, Li Y. Biomechanical properties and collagen fiber orientation of temporomandibular joint discs in dogs: 2. Tensile mechanical properties of the discs. J Craniomand Disorders 1991;5:107–14. [10] Chin LPY, Aker FD, Zarrinnia K. The viscoelastic properties of the human temporomandibular joint disc. J Oral Maxillofac Surg 1996;54:315–8. [11] Lai W-FT, Bowley J, Burch JG. Evaluation of shear stress of the human temporomandibular joint disc. J Orofac Pain 1998;12: 153–9. [12] Tanaka E, Tanaka M, Miyawaki Y, Tanne K. Viscoelastic properties of canine temporomandibular joint disc in compressive load-relaxation. Arch Oral Biol 1999;44:1021–6. [13] Eckstein F, Lemberger B, Stammberger T, Englmeier KH, Reiser M. Patellar cartilage deformation in vivo after static versus dynamic loading. J Biomech 2000;33:819–25. [14] Beek M, Aarnts MP, Koolstra JH, Feilzer AJ, Van Eijden TMGJ. Dynamical properties of the human temporomandibular joint disc. J Dent Res 2001;80:876–80. [15] Mow VC, Wang CC-B. Some bioengineering considerations for tissue engineering of articular cartilage. Clin Orthop Rel Res 1999;367S:S204–23. [16] Mow VC, Kuei SC, Lai WM, Armstrong CG. Biphasic creep and stress relaxation of articular cartilage in compression: Theory and experiments. J Biomech Eng 1980;102:73–84. [17] Soltz MA, Ateshian GA. Experimental verification and theoretical prediction of cartilage interstitial fluid pressurization at an impermeable contact interface in confined compression. J Biomech 1998;31:927–34.

76

M. Beek et al. / Clinical Biomechanics 18 (2003) 69–76

[18] Prendergast PJ, Van Driel WD, Kuiper J-H. A comparison of finite element codes for the solution of biphasic poroelastic problems. Proc Inst Mech Eng 1996;210(H):131–6. [19] Simon BR. Multiphasic poroelastic finite element models for soft tissue structures. Appl Mech Rev 1992;45:191–218. [20] Wu JZ, Herzog W, Epstein M. Evaluation of the finite element software ABAQUS for biomechanical modeling of biphasic tissues. J Biomech 1998;31:165–9. [21] Mow VC, Holmes MH, Lai WM. Fluid transport and mechanical properties of articular cartilage: a review. J Biomech 1984;17: 377–94. [22] Fithian DC, Kelly MA, Mow VC. Material properties and structure-function relationships in the menisci. Clin Orthop Rel Res 1990;252:19–31. [23] Mansour JM, Mow VC. The permeability of articular cartilage under compressive strain and at high pressures. J Bone J Surg A 1976;58:509–16. [24] Bursac PM, Obitz TW, Eisenberg SR, Stamenovic D. Confined and unconfined stress relaxation of cartilage: appropriateness of a transversely isotropic analysis. J Biomech 1999;32: 1125–30. [25] Suh J-K, Li Z, Woo SL-W. Dynamical behavior of a biphasic cartilage model under cyclic compressive loading. J Biomech 1995;28:357–64. [26] Soltz MA, Ateshian GA. Interstitial fluid pressurization during confined compression cyclical loading of articular cartilage. Ann Biomed Eng 2000;28:150–9. [27] Kim Y-J, Bonassar LJ, Grodzinsky AJ. The role of cartilage streaming potential, fluid flow and pressure in the stimulation of chondrocyte biosynthesis during dynamic compression. J Biomech 1995;28:1055–66.

[28] Cohen NP, Foster RJ, Mow VC. Composition and dynamics of articular cartilage: structure, function, and maintaining healthy state. J Orthop Sports Phys Ther 1998;28:203–15. [29] Arokoski JPA, Hyttinen MM, Helminen HJ, Jurvelin JS. Biomechanical and structural characteristics of canine femoral and tibial cartilage. J Biomed Mat Res 1999;48:99–107. [30] Kopp S. Topographical distribution of sulphated glycosaminoglycans in human temporomandibular joint disks. A histochemical study of an autopsy material. J Oral Pathol 1976;5:265–76. [31] Piette E. Anatomy of the human temporomandibular joint. An updated comprehensive review. Acta Stomatol Belg 1993;90: 103–27. [32] Minarelli AM, Del Santo Junior M, Liberti EA. The structure of the human temporomandibular joint disc: a scanning electron microscopy study. J Orofac Pain 1997;11:95–100. [33] Tanaka E, Shibaguchi T, Tanaka M, Tanne K. Viscoelastic properties of the human temporomandibular joint disc in patients with internal derangement. J Oral Maxillofac Surg 2000;58:997– 1002. [34] Gore DM, Higginson GR, Minns RJ. Compliance of articular cartilage and its variation through the thickness. Phys Med Biol 1983;28:233–47. [35] Athanasiou KA, Niederauer GG, Schenck Jr RC. Biomechanical topography of human ankle cartilage. Ann Biomed Eng 1995; 23:697–704. [36] Schinagl RM, Gurskis D, Chen AC, Sah RL. Depth-dependent confined compression modulus of full-thickness bovine articular cartilage. J Orthop Res 1997;15:499–506. [37] Jurvelin JS, Arokoski JP, Hunziker EB, Helminen HJ. Topographical variation of the elastic properties of articular cartilage in the canine knee. J Biomech 2000;33:669–75.