Journal of Biomedical Materials Research. Part B, Applied Biomaterials
John Wiley & Sons, Inc.
Computational characterization of the porous‐fibrous behavior of the soft tissues in the temporomandibular joint
Volume: 108, Issue: 5
DOI 10.1002/jbm.b.34558
•
•
•
• Altmetric

### Notes

Ortún‐Terrazas, Cegoñino, and Pérez del Palomar: Computational characterization of the porous‐fibrous behavior of the soft tissues in the temporomandibular joint

## 1. INTRODUCTION

The temporomandibular joint (TMJ) is an articular joint that connects the mandible to the skull and enables jaw movement. From a biomechanical point of view, some of the most important tissues in this joint are the articular disc, the condylar and temporal cartilage, the retrodiscal tissue, and the collateral ligaments (Figure 1). These soft tissues constrain jaw movement and are the tissues most commonly implicated in temporomandibular joint disorders (TMDs). Although TMDs have a multifactorial etiology, approximately 70% of TMDs (Detamore & Athanasiou, 20052005) and/or TMJ malformations are closely linked to biomechanical anomalies (Allen & Athanasiou, 20062006) that end up being resolved through joint replacements. A deeper knowledge of the biomechanical behavior of the soft tissues of the TMJ would facilitate the development of more accurate finite element (FE) models, more effective treatments and designs of implants closer to human joint.

Figure 1
Schematic showing the anatomy of the temporomandibular joint. (a) TMJ disc, (b) condylar cartilage, (c) collateral ligament, and (d) retrodiscal tissue. In each panel, the alignment of collagen networks is shown on the left, and the samples analyzed in this study are indicated on the right

The TMJ disc is a biconcave, elliptical, fibrocartilaginous tissue located in the glenoid fossa of the temporal bone (Figure 1). This tissue distributes load and absorbs shock, preventing stress concentration on the articulating surfaces of the TMJ (condylar and temporal cartilage). The TMJ disc slides between the articular surfaces of the temporal and condylar bone regions which are covered by hyaline cartilage (Singh & Detamore, 20092009) (Figure 1). In normal conditions, the condylar cartilage is mainly subjected to compression stress and the shear stress is often neglected due to extremely low friction (Nickel & McLachlan, 19941994) between the contact surfaces which are lubricated by the synovial fluid. The retrodiscal tissue contains collagen and elastic fibers and connects the posterior band of the disc to the mandibular fossa (Tanaka, Del Pozo, Sugiyama, & Tanne, 20022002). Finally, the medial and lateral collateral ligaments connect the medial and lateral sides of the TMJ disc to the condyle, thereby contributing to joint stabilization during TMJ movement.

Structurally and biomechanically, the TMJ disc (Wright et al., 20162016) and fibrous region of the condylar cartilage (Singh & Detamore, 20092009) can be divided in five regions (Figure 1a,b): the anterior band (ANT), posterior band (POST), central zone (CENT), lateral zone (LAT), and medial zone (MED). Experimental studies (Singh & Detamore, 20082008; Wright et al., 20162016) have shown that collagen fiber density and orientation account for mechanical differences between these regions. In the disc attachments (retrodiscal and collateral ligaments), the fibrous network is particularly dense and is oriented in an anteroposterior direction (Murphy, Arzi, Hu, & Athanasiou, 20132013).

Because experimental manipulation of the TMJ is hampered by the fragility of this joint's soft tissues, the last decade has seen an increase in the number of computational studies of the TMJ (Ingawalé & Goswami, 20092009) (shown in Table S1). However, the mechanical properties of the soft tissues of the TMJ remain to be fully characterized (Ingawalé & Goswami, 20092009), as underscored by the many contrasting findings in the published literature. This variability may be due to viscoelastic behavior is not considered; the dependence of biomechanical behavior on fiber orientation is not considered; and mechanical tests are performed using animal species with a very different masticatory system to that of the human. Therefore, experimental studies that have considered the influence of fibrous network (Murphy et al., 20132013; Singh & Detamore, 20082008; Wright et al., 20162016) and the time‐dependence behavior (Beek, Aarnts, Koolstra, Feilzer, & Van Eijden, 20012001) are particularly important for the present work.

Moreover, few of the computational studies have included all TMJ soft tissues in their FE models (shown in Table S1) and only a few have examined the behavior of the collagen network. Furthermore, no studies have characterized the biomechanical contributions of both, the fibrous network and the interstitial fluid of the TMJ disc, using experimental data obtained from a species with a masticatory system similar to that of the human (Ingawalé & Goswami, 20092009; Wright et al., 20162016).

Given the discrepancies in the literature regarding some of the key biomechanical properties of the soft tissues of the TMJ, we sought to identify and rigorously validate material model parameters for the characterization of the soft tissues of the TMJ. The aims of this study were: (a) to characterize the mechanical parameters that describe the fibrous and porous behavior of the TMJ disc and the role of the fibrous network in other soft tissues of the TMJ (condylar cartilages, retrodiscal tissue, and collateral ligaments); (b) to validate the material properties by comparing our numerical results with experimental data from the literature; (c) to evaluate the simplifications performed during the material model characterization by considering different configurations of the fibrous network. In this study, firstly the characterization of the fibrous network of the TMJ disc, retrodiscal tissue, condylar cartilage, and collateral ligaments is explained. The behavior of the temporal cartilage was not characterized owing to a lack of sufficient data. Next, the contribution of TMJ disc fluid is evaluated by reproducing compression tests described in the literature. Based on the respective contributions of the fibrous network and the TMJ disc fluid, a material model is proposed. Finally, we describe the evaluation of TMJ soft tissue responses for different directions, configurations, and strain rates.

## 2. MATERIALS AND METHODS

### 2.1. Characterization of the fibrous matrix

In soft tissues studied, the collagen network is the main component responsible for tensile stiffness, while the interstitial fluid resists most of the compressive load (Beek et al., 20012001; Wright et al., 20162016). The tensile stiffness of the tissue varies according to the preferential direction of fiber orientation. Thus, collagen fibers contribute significantly to tissue stiffness when stretched, but make little contribution when transverse loads are applied. In this latter situation, the majority of the load must be borne by the hyperelastic isotropic matrix.

Considering x and X, the respective positions of a particle in a current (Ω) and reference configurations 0), the deformation gradient F is defined by F = x/X. Alternatively, the modified deformation gradient F~ can be expressed by the Jacobian deformation gradient (Jel = detF) and the deformation gradient (F) by the formula F~=Jel1/3F. The right Cauchy‐Green tensor C~ can therefore be defined by the modified deformation gradient, F~, as C~=F~TF~. To account for the aforementioned behaviors of the matrix and the fibrous network, we used a transversely isotropic hyperelastic material model (Holzapfel, 20002000) which is expressed by the following strain energy density function (W):

$W\left(\stackrel{~}{\mathbit{C}},\stackrel{~}{{I}_{1}},\stackrel{~}{{I}_{4}}\right)={W}_{m}\left(\stackrel{~}{{I}_{1}}\right)+{W}_{f}\left(\stackrel{~}{{I}_{4}}\right)+{W}_{\mathit{vol}}\left({J}_{\mathit{el}}\right)$
where Wm is the component of the hyperelastic matrix substance without fibers, Wf represents the fibrous term, and Wvol is the volumetric term. Wm is defined according to the deviatoric invariant which can be expressed in term of C~ eigenvalues λ^1λ^2λ^3 as follows:
${\stackrel{~}{I}}_{1}=\mathit{tr}\stackrel{~}{\mathbit{C}}={\phantom{\rule{0.25em}{0ex}}\stackrel{^}{\lambda }}_{1}^{2}+{\stackrel{^}{\lambda }}_{2}^{2}+{\stackrel{^}{\lambda }}_{3}^{2}$

On the other hand, Wf is a function of the dilatational invariant which is defined as:

${\stackrel{~}{I}}_{4}={\mathbit{a}}^{0}·\stackrel{~}{\mathbit{C}}·{\mathbit{a}}^{0}$
where a0 is the unitary vector defining the orientation of the collagen fibers in the reference configuration. Therefore, the terms of the energy density function (Equation (1) can be reformulated as follows:
${W}_{m}={C}_{1}·\left({\stackrel{~}{I}}_{1}-3\right);$
$\phantom{\rule{0.25em}{0ex}}{W}_{f}=\frac{{k}_{1}}{2·{k}_{2}}\left\{\mathit{exp}\left[{k}_{2}·{\left({\stackrel{~}{I}}_{4}-1\right)}^{2}\right]-1\right\}$
${W}_{\mathit{vol}}=\frac{1}{D}\left(\frac{{\left({J}_{\mathit{el}}\right)}^{2}-1}{2}-\mathit{ln}{J}_{\mathit{el}}\right)$
where Wm is defined by the Neo‐Hookean hyperelastic material formulation and C1 is the material constant related to the matrix substance. k1 > 0 and k2 > 0 are the parameters which identify the exponential behavior due to the presence of collagen fibers, and D is related to bulk modulus K, by D = 2/K being zero for an incompressible material. For quasi‐incompressible materials (Jel ≈ 1), the volumetric term can be underestimated.

In an uniaxial stretch test, axial stretch ratio λu can be defined by the C~ eigenvalues as λ1 = λu and λ2=λ3=1/λu due to incompressible condition (Jel = λuλ22 = 1). Hence, the nominal stress‐stretch relation (Holzapfel, 20002000) for the transversely isotropic material model can be expressed as follows:

${\mathrm{\sigma }}_{u}=\frac{\partial W}{\partial {\lambda }_{u}}={{\lambda }_{u}}^{2}\left\{{C}_{1}\left(1-\frac{1}{{{\lambda }_{u}}^{3}}\right)+2·{k}_{1}\left({\stackrel{~}{I}}_{4}-1\right)·\mathit{exp}\left[{k}_{2}·{\left({\stackrel{~}{I}}_{4}-1\right)}^{2}\right]\right\}\phantom{\rule{0.25em}{0ex}}$

When the collagen fibers are almost not stretched I~41 the contribution of the fibrous network is practically null and the hyperelastic matrix (first term of Equation (5)) must support most of the load. Conversely, most of the load is supported by the fibers if the tissue is stretched in the direction of the fibrous network (I4~>1). Therefore, it is possible to determine separately the material parameters that define both terms (Wm and Wf) by fitting the stretch experimental response of samples tested longitudinal and transversal to the direction of fibers network. The parameters of Equation (5) were adjusted through a nonlinear least‐squares fitting procedure using MATLAB commercial software (MATLAB 6.0 R12, The MathWorks Inc., Natick, MA, 2000), which minimizes the relative error between the analytical values of Equation (5) and the experimental fitted curves. The toe region (strain <6%) (Singh & Detamore, 20082008) of the fitted curves was defined by the transversal elastic modulus and was used to the characterization of C1. Meanwhile, k1 and k2 parameters were obtained by fitting the near‐linear region (strain >6%) which was defined based on the longitudinal elastic modulus of the experimental studies that are described below.

For characterization of the TMJ disc, we used the experimental data of Wright et al. (20162016) in which 24 human TMJ discs were tested in anteroposterior and mediolateral directions, analyzing the effects of loading direction, region, and sex. Their histological findings revealed that collagen fibers are oriented to form a peripheral ring (Figure 1a). Given that central fibers run in the anteroposterior direction, the matrix parameters can be computed from the elastic modulus of central region tested mediolaterally. Therefore, we defined Wm based on the elastic modulus of the ML‐CENT sample. The parameters of the fibrous term (Wf) of each region were characterized using samples that were stretched along the preferred direction of the fibrous network of each region. The elastic modulus of AP‐MED, AP‐CENT, AP‐LAT, ML‐POST, and ML‐ANT samples was used for the characterization of the fibrous term in medial, central, lateral, posterior, and anterior regions, respectively. Furthermore, the parameters for each region were characterized separately for men and women using experimental data taken from Wright et al. (20162016).

To characterize the behavior of condylar cartilage, we used experimental data from a study by Singh and Detamore (20082008), who tested seven samples of porcine condylar cartilage in anteroposterior and mediolateral directions. Based on the orientation of the cartilage fibers (Figure 1b), Wm parameters were computed from the ML‐CENT elastic modulus. Wf parameters for each of the five regions of the condylar cartilage were computed from the elastic modulus of the samples tested in the direction of fiber orientation.

Finally, the material parameters of collateral ligaments and retrodiscal tissue were computed based on experimental data from Murphy et al. (20132013). That study quantified the elastic modulus of porcine TMJ disc attachments in both anteroposterior and mediolateral directions (Figure 1c,d). In both tissues the fibers are mainly oriented in the anteroposterior direction. Therefore, Wm and Wf were computed based on the results obtained from samples that were tested mediolaterally and anteroposteriorly, respectively. The material parameters obtained are summarized in Table 1.

Table 1
Left to right: parameters for the transversely isotropic hyperelastic material model for the TMJ disc (man and woman), condylar cartilage, medial and lateral collateral ligaments, and retrodiscal tissue; parameters for the hyperfoam material model, which characterizes the compressive response of the solid phase of the TMJ disc; parameters for biphasic characterization of the TMJ disc. Note: (a) Kuo et al. (2010)
Solid phasePorous phase
Trans. Iso. Hype. mat. modelHyperfoam mat. modelPermeability function parameters
C1 (MPa)D (MPa−1)k1 (MPa)k2 (−)ν (−)μ (MPa)α (−)k0 (10−15 m2)M (−)φ (−)γw (N/m3)
TMJ disc (man)
ANT1.4500.430.340.16671.72758.95 ± 1.17 (a)490.80 (a)9800
LAT1.4500.690.430.16673.86783.64 ± 0.65 (a)450.79 (a)9800
CENT1.4500.970.17
MED1.4500.171.68
POST1.4501.250.160.16670.86758.07 ± 1.55 (a)430.83 (a)9800
TMJ disc (woman)
ANT2.400.053.720.16671.72758.95 ± 1.17 (a)490.81 (a)9800
LAT2.400.112.520.16673.86783.64 ± 0.65 (a)450.79 (a)9800
CENT2.400.750.87
MED2.400.082.93
POST2.400.311.440.16670.86758.07 ± 1.55 (a)430.83 (a)9800
Condylar cartilage (porcine)
ANT1.6500.241.95
LAT1.6502.580.43
CENT1.6503.770.21
MED1.6502.520.42
POST1.6500.161.92
Medial collateral ligaments (porcine)
UNIF.0.3800.6810−5
Lateral collateral ligaments (porcine)
UNIF.0.4100.3510−5
Retrodiscal tissue (porcine)
UNIF.1.0500.2710−5

### 2.2. Characterization of the porous parameters of the TMJ disc

To account for the biphasic nature of the TMJ disc, the poroelastic contribution was added in Abaqus software (Abaqus 6.14, Simulia, Rhode Island) to the transversely isotropic material model defined in Equation (1). In a biphasic material (Holmes & Mow, 19901990), as here studied, media is composed of solid (ϕs) and fluid phases (ϕf). The sum of both phases in fully saturated media, ζ, must satisfy ζ=ϕs + ϕf = 1. Abaqus provides capabilities for modelling a biphasic material when the media is completely saturated with an incompressible fluid by the conventional poroelastic theory which considers the medium as a multiphase media with the interaction between the solid and the fluid phases. The balanced equations for mass, linear momentum and energy in soil mechanics (Wu, 19671967) are additionally provided in the supplementary material. This approach uses the principle of effective stress which states that the total stress acting in a mixture can be decomposed into two components. The neutral stress, also called pore pressure, p, which acts in every direction and the effective stress, σ¯, whose effect is exclusively on the solid matrix. Therefore, the total stress acting at a point, σ, is defined from the Clausius‐Duhem entropy inequality as follows:

$\mathbit{\sigma }\mathbf{=}\left(1-e\right)\overline{\mathbit{\sigma }}-e·p·\mathbit{I}-\zeta ·p·\mathbit{I}$
where I is the identity tensor, p is the pore pressure, and σ¯ is the effective stress tensor which is determined from the strain energy density function for a hyperelastic materials by (Systèmes, 20142014):
$\overline{\mathbit{\sigma }}=\frac{2}{{J}_{\mathit{el}}}\stackrel{~}{\mathbit{F}}\frac{\partial W}{\partial \stackrel{~}{\mathbit{C}}}{\stackrel{~}{\mathbit{F}}}^{\mathrm{T}}$

To compute p, Abaqus requires to specify the permeability, k¯, as a function of the void ratio which is the ratio of fluid volume to total volume of the element (e = dVf/dVt). In our simulations, the exponential permeability function for biphasic materials (Argoubi & Shirazi‐Adl, 19961996) was used:

$\overline{k}={k}_{0}{\left[\frac{e\left(1+{e}_{0}\right)}{{e}_{0}\left(1+e\right)}\right]}^{2}\mathit{exp}\left[M\left(\frac{1+e}{1+{\mathrm{e}}_{0}}-1\right)\right]$
where e is the strain‐dependent void ratio, M is a dimensionless material parameter and, k0 and e0 are the permeability and the void ratio at zero strain (shown in Table 1). The initial void ratio in each region of the disc,e0=ϕ0f/1ϕ0f, was computed from the porosity values at zero strain proposed in the literature (Kuo, Zhang, Bacro, & Yao, 20102010). The initial permeabilities of each region were also determined from this experimental study (Kuo et al., 20102010). The void ratio at each instant of time, e, is dependent on the deformation (Systèmes, 20142014) and was obtained from the computational results of our simulations according to the following equation:
$e=\frac{1-{{J}_{\mathit{el}}}^{-1}\left(1-{\varphi }_{0}^{f}\right)}{{{J}_{\mathit{el}}}^{-1}\left(1-{\varphi }_{0}^{f}\right)}$

On the other hand, M parameter was iteratively determined through a least‐squares fitting procedure, minimizing the relative error between the nominal stress obtained from our compressive simulations and the experimental compressive response. The computational nominal stress values in each fitting‐iteration were obtained through an iterative compression test of the cylindrical 3D model which is subsequently explained. This iterative process was executed in Python (Python 3.5.2, Python Software Foundation). Meanwhile, the experimental data was obtained from Beek et al. (20012001) study in which human samples were tested considering both strain‐rate and region dependencies.

As it will show in Section 3 (Figure 6a), the porous transversely isotropic material model did not fit well the compression data. This incompatibility has been previously reported for other soft tissues of the masticatory system, such as the periodontal ligament (Ortún‐Terrazas, Cegoñino, & Pérez del Palomar, 20202020; Ortún‐Terrazas, Cegoñino, Santana‐Penín, Santana‐Mora, & Pérez del Palomar, 20182018; Ortún‐Terrazas, Cegoñino, Santana‐Penín, Santana‐Mora, & Pérez del Palomar, 20192019). This poor fit is largely due to the fact that a Neo‐Hookean hyperelastic material formulation is not valid for high compressive stretch values (Ogden, 19721972; Yeoh, 19931993). Consequently, we developed a new material model defined by a porous transversely isotropic hyperelastic material description for tensile efforts I~4>1, and a porous hyperfoam material description for compressive efforts I~41. This constitutive material model was implemented in a UMAT user subroutine, as previously described (Ortún‐Terrazas et al., 20182018).

The hyperfoam strain energy density function (Storakers, 19861986) is defined as follows:

$W=\frac{2\mu }{{\alpha }^{2}}\phantom{\rule{0.25em}{0ex}}\left[{\stackrel{^}{\lambda }}_{1}^{\alpha }+{\stackrel{^}{\lambda }}_{2}^{\alpha }+{\stackrel{^}{\lambda }}_{3}^{\alpha }-3+\frac{1}{\beta }\left({J}_{\mathit{el}}^{-\mathit{\alpha \beta }}-1\right)\right]$

Where μ and α are material parameters and the coefficient β determines the degree of compressibility being related to Poisson's ratio, ν, by β = ν/(1 − 2ν). In an uniaxial compression test of a compressible isotropic material (Jel ≠ 1), the Jacobian deformation gradient can be then defined as leading the following nominal stress‐stretch relation (Systèmes, 20142014):

${\mathrm{\sigma }}_{u}=\frac{\partial W}{\partial {\lambda }_{u}}=\frac{2\mu }{{\lambda }_{u}\alpha }\left({\stackrel{^}{\lambda }}_{u}^{\alpha }-{J}_{\mathit{el}}^{-\mathit{\alpha \beta }}\right)$

In our model, the coefficient β was computed through a Poisson ratio, ν = 0.1667, suggested in Beek, Koolstra, and van Eijden (20032003) for the whole disc. Meanwhile, μ and α parameters were defined for each region through a nonlinear least‐squares fitting procedure in MATLAB, which minimizes the relative error between the analytical values of Equation (11) and the experimental fitted curves. The input experimental data for fitting procedure were obtained from the unloading curves of anterior, intermediate and posterior tests at 0.05 Hz performed by Beek et al. (20012001). Additionally, the upper limit of μ variation range was limited to the shear modulus of maximum stretch value because of the relation between μ and the shear modulus (Systèmes, 20142014).

### 2.3. FE models for validation tests

To validate the material parameters computed previously, the experimental data taken from the literature and used for the characterization were mimicked numerically in Abaqus. In this section, we first summarize the geometry of the test samples and then describe the FE models used by us to simulate the experimental tests. We simulated the tensile tests (Murphy et al., 20132013; Singh & Detamore, 20082008; Wright et al., 20162016) using a rectangular 2D FE model (Figure 2a), and the compressive test (Beek et al., 20012001) using a cylindrical 3D model (Figure 2b).

Figure 2
(a) Top: geometry and FE model mesh of the rectangular specimen used to mimic tensile tests. Bottom: schematic showing fiber orientation and the regions of the TMJ disc, condylar cartilage, collateral ligaments, and retrodiscal tissue. (b) 3D FE model of the cylindrical test sample used for compression analysis and schematic showing the boundary conditions applied

In the stress‐relaxation tests described in Wright et al. (20162016), rectangular specimens from six female and six male human TMJ discs were stretched (Figure 1a). The mean thickness of the samples was 2.2 mm in the anteroposterior (AP) direction and 2.5 mm in the mediolateral (ML) direction. In our simulations, sample width (w ) (Figure 2a) was set at 2 mm and mean lengths (l0) at 7.1 mm and 10.6 mm for samples tested in the anteroposterior and mediolateral directions, respectively.

In the study by Beek et al., several cylindrical samples of human TMJ discs were compressed (Beek et al., 20012001) by two impermeable indenters with a radius of 1.97 mm (rind). In a subsequent computational study (Beek et al., 20032003), the authors calculated that the cylindrical samples had an external radius of 3.49 mm (rdisc) and a thickness of 2.14 mm (h0 ) (Figure 2b).

In experiments performed by Singh and Detamore (20082008), the mean dimensions of the condylar cartilage samples tested in mediolateral and anteroposterior directions, respectively, were as follows: thickness, 0.47 and 0.49 mm; length, 9.6 and 9.1 mm; width, 2.0 and 1.9 mm.

Regarding the dimensions of the samples stretched by Murphy et al. (20132013), they reported a mean thickness of 1.7 mm for both retrodiscal tissue and collateral ligaments. Width (1.9 mm) and length (12 mm) values were the same in the anteroposterior and mediolateral directions.

To simulate tensile experimental tests, rectangular samples were 2D modeled because they were less thick than they were long/wide. To define different fiber orientations, each FE model was divided into 10 sections by making five cuts along the length of the sample and one along the midline. The preferential direction of the fiber bundles, a0, was defined using a rectangular reference system (Figure 2a). Each sample was meshed by applying a mesh convergence process. This process evaluates changes in the mesh in samples with the largest and smallest geometric dimensions (i.e., geometry of TMJ disc attachments and AP samples of the TMJ disc, respectively). After the mesh convergence analysis, the FE model consisted of 550 second‐order quadrilateral elements (S8‐type element in Abaqus). In all cases, the longitudinal movement of one end of the sample was restricted in the direction of the load. The opposite side was displaced by 20, 30, 100, and 110% of the sample length, in accordance with the maximum stretch values described in the experimental test mimicked.

To evaluate the sensibility of considering different fibers properties and orientations in the same sample, two different scenarios were simulated for the samples of TMJ disc and the condylar cartilage. The first modeled a single preferential direction of fiber orientation in whole sample (direction of the central sections of the samples shown in Figure 2a), while the second assumed that each sample consisted of three different regions with distinct fiber properties and orientations (Figure 2a). The sections with an oblique orientation shown in Figure 2a were defined with a mean angle φ of 30° for vector a0. The mean angle was computed based on polarized light microscopy images from studies by Wright et al. (20162016) and Singh and Detamore (20082008). Owing to a lack of published experimental data, we made no distinction between regions for collateral ligaments and retrodiscal tissue. Nonetheless, the dependence of their mechanical behavior on fibrous network direction was evaluated by considering different fiber orientations (0, 20, 40, 70, and 90°). The material parameters shown in Table 1 were assigned to each section according to the color code shown in Figure 2a.

Regarding the cylindrical 3D model used to mimic the compressive experimental test, it was meshed with 53,940 porous second‐order elements (C3D20P‐type element in Abaqus). The mesh size was also determined based on a mesh convergence analysis, as described above, and was finer in the intermediate region (Figure 2b). The response of the test sample was recorded for both the porous transversely isotropic hyperelastic material model and the porous hyperfoam material model. In the case of the porous transversely isotropic hyperelastic material model, the orientation of the fiber network was considered perpendicular to the load (φ = 00).

To reproduce Beek et al. (20012001) experimental compression test, of interstitial fluid through the free sides of the disc was allowed and the nodes of the upper surface of the computational sample were thus subjected to a sinusoidal (ε(t) = ε · sin(2πft)) compressive displacement equivalent to 35, 30, and 25% strain (ε). The compressive displacement cycle was evaluated at 0.02, 0.05, and 0.1 Hz. Two full cycles were simulated to allow time for model stabilization (Bergomi et al., 20112011). The second cycle was used for comparison with the experimental data. The compression test was performed for the anterior, posterior, and intermediate regions of the TMJ disc.

## 3. RESULTS

The results of the numerical mimicking of the selected experimental data are presented here. As aforementioned, in our study, the properties of the transversely isotropic hyperelastic material model were characterized assuming the orientation of all fibers in the same direction throughout the whole sample. This assumption was computationally evaluated by comparing results obtained considering one or several concurrent fibrous network orientations and region properties. Figures 3, 4, 5 show the different outcomes obtained for models defined with a single fiber orientation versus with multiple regions with distinct fiber orientations. Numerical tensile responses and Normalized Root Mean Squared Error (NRMSE) from data of Wright et al. (20162016), Singh and Detamore (20082008) and Murphy et al. (20132013) are respectively shown in Figures 3, 4, 5. On the other hand, the results of the evaluation of the compressive response and its associated NRMSEs are shown in Figure 6. Figure 6a shows the fitted curve for the experimental data of Beek et al. (20012001), at 0.05 Hz, while Figure 6bd) shown the different mechanical responses obtained using the porous transversely isotropic hyperelastic and porous hyperfoam material models. Additionally, Figure S1 shows the minimum principal stress values of our simulations and Beek et al. (20012001) experimental results for different compressive strain values (35, 30, and 25%) and frequencies (0.1, 0.05, and 0.02 Hz).

Figure 3
Comparison of the tensile response computed from the planar FE model sample by using a single (continuous line) or multiple (discontinuous line) fibers orientations regarding the experimental data (Wright et al., 20162016) for the TMJ disc of man (green) and woman (purple). Regions: ANT, anterior; CENT, central; POST, posterior; LAT, lateral; MED, medial; Directions: AP, anteroposterior; ML, mediolateral
Figure 4
Comparison of the tensile response computed from the planar FE model sample by using a single (continuous line) or multiple (discontinuous line) fibers orientations regarding the experimental data (Singh & Detamore, 20082008) for the condylar cartilage. Regions: ANT, anterior; CENT, central; POST, posterior; LAT, lateral; MED, medial; Directions: AP, anteroposterior; ML, mediolateral
Figure 5
Comparison of the tensile response computed from the planar FE model sample by using a single (continuous line) or multiple (discontinuous line) fibers orientations regarding the experimental data (Murphy et al., 20132013) for the (a) medial collateral ligament; (b) lateral collateral ligament, and (c) retrodiscal tissue. Directions: AP, anteroposterior; ML, mediolateral
Figure 6
(a) Compressive response of the anterior, intermediate, and posterior regions of the TMJ (0.05 Hz) evaluated using a porous transversely isotropic hyperelastic model (narrow discontinuous line) and a porous hyperfoam material model (dotted line); distribution of (b) minimum principal stress; (c) porous pressure, and (d) void ratio in the cylindrical FE model using the porous transversely isotropic hyperelastic material model (left) and the porous hyperfoam material model (right)

As shown in Figure 3, in both cases the numerical responses fall within the experimental deviation with NRMSEs between 0.88–3.46% and 0.95–9.61% for simulations with singular and multiple fiber orientations, respectively. For both sexes, the greatest deviation between both approaches is observed for the AP‐CENT sample. No differences were observed for the ML‐CENT sample, the fibers were not stretched and the C1 parameter in all regions was characterized using data from the same sample (ML‐CENT).

In the case of the condylar cartilage sample under tensile load, marked differences were observed between the two sample types (Figure 4). These differences were particularly significant for the three samples that were stretched in the anteroposterior direction (AP‐LAT, AP‐CENT, and AP‐MED) with NRMSE values higher than 10%. As a result of the orientation of the fibrous network in these three samples, the anterior and posterior regions are stretched more than the central region. The fibers of the anterior and posterior regions pose almost no resistance because they are not oriented in the direction in which the sample is stretched. Consequently, the degree of stretching in these two bands is greater, leading to a reduction in the stiffness of the sample.

The effect of fiber network orientation on the tensile mechanical response of disc attachments is illustrated in Figure 5. As its analytical statement describes, tissue stiffness increases when the fibers are oriented in the direction of the load. Conversely, when the fibers are oriented transversally relative to the direction of load, the resulting behavior resembles that observed for the hyperelastic material model.

Despite the differences between both approaches, the NRMSEs are lower than 7.2% for all simulations that considered singular orientation of the collagen network and all material responses fall within the experimental deviation.

Figure 6 shows the compressive response of the TMJ disc obtained using the cylindrical FE model. An unfitted compressive response was obtained using the porous transversely isotropic hyperelastic material description to describe the experimental data of Beek et al. (20012001). In this scenario, energy dissipation during the loading–unloading cycle (known as hysteresis) was almost inappreciable, possibly owing to the quasi‐incompressible behavior of the tissue and the minimal contribution of the fluid phase. Meanwhile, the porous hyperfoam model produced a better curve fit, and more marked hysteresis, as illustrated in the color maps in Figure 6bd. Comparison of the two models reveals that the minimum principal stress of the solid matrix (Figure 6b) is greater for the porous transversely isotropic material model, while the interstitial fluid pressure is greater for the porous hyperfoam material model (Figure 6c). This outcome demonstrates the marked contribution of the fluid phase in the porous hyperfoam model. This contribution is also reflected in the capacity of the fluid to fill surrounding areas, resulting in a higher void ratio gradient in the porous hyperfoam model (Figure 6d). To verify the porous effect using the porous hyperfoam model, we also compared results obtained at different frequencies (Figure S1) with those of Beek et al. (20012001). Although, the computed response of the intermediate region at a high strain rate (0.05 Hz) and strain (35%) is the most deviated from the experimental results of Beek and coworkers, all values fell within the experimental deviation.

## 4. DISCUSSION

The first part of this study describes the characterization of the fibrous properties of the soft tissues of the TMJ using a transversely isotropic hyperelastic material model (Holzapfel, 20002000). The material parameters were computed for the collateral ligaments, the retrodiscal tissue, the five regions of the TMJ disc, and the condylar cartilage. While many studies have sought to characterize these tissues (see Table S2), the existing literature contains many conflicting findings, owing to the complex behavior of these tissues.

Several experimental studies have investigated the biomechanical behavior of the TMJ disc. While most of these studies take into account the biomechanical differences between different regions of the disc (shown in Table S2), few consider the influence of fiber orientation on the outcomes obtained (shown in Table S2). Thus, in the present study, we examined the biomechanical responses of the TMJ disc in the context of the different potential orientations of the collagen fiber network. The anterior and posterior bands of the TMJ disc are stiffer in the mediolateral direction because of the ring‐like distribution of the collagen fibers in these regions. By contrast, the remaining regions are stiffer in the anteroposterior direction owing to the anteroposterior orientation of the fibers. These properties can be related to the physiological function of the disc. When the disc is compressed against the mandibular fossa by the chewing muscles, the ring‐shaped fibers help maintain the shape of the disc (Angelo et al., 20162016; Detamore & Athanasiou, 20032003). On the other hand, the fibers of the central region (mainly anteroposteriorly oriented) support protrusion and retrusion movements (Palla, Gallo, & Gossi, 20032003) of the mandible. The fact that both the non‐fibrous hyperelastic matrix and the collagen network are stiffer in women than in men may have etiological implications for the onset and progression of TMDs in women (Wright et al., 20162016).

Despite inconsistencies in the literature regarding the compressive properties of the TMJ disc, it is clear that the disc response is softer in compression than in tension (Allen & Athanasiou, 20062006). The study by Beek et al. (20012001) is one of the most important experimental studies of the compressive response of the TMJ published to date. In that study, the authors demonstrated that the compressive response is dependent on strain‐rate, a relationship they propose may help explain some of the discrepancies in the literature. They examined the effects of compression in different regions (intermediate, anterior, and posterior bands) of seven human TMJ discs. When the tissue is compressed, the fibers are not stretched, and the tissue softens, resulting in the absorption of most of the stress by the interstitial fluid. In the present study, we included the porous contribution in a transversely isotropic hyperelastic material model in order to quantify the contribution of the interstitial fluid. However, the porous formulation of the fiberless matrix defined by a Neo‐Hookean hyperelastic material formulation does not fit well with the compressive experimental data of Beek et al. (20012001). We therefore assessed the compressive response of the tissue using a porous hyperfoam material model, which allowed us to characterize the hysteresis behavior generated by the viscoelastic effects of the interstitial fluid. Using this material model, we also correctly mimicked the compressive responses of the TMJ disc at different frequencies and strain rates.

We propose therefore a complex model for TMJ soft tissues that considers the contribution of collagen fibers and the compression of the interstitial fluid. This model is based on a model generated by our group for another soft tissue, the periodontal ligament (Ortún‐Terrazas et al., 20182018, 20192019). When the tissue is compressed I~41 the strain energy function is defined by the porous hyperfoam material model (Equation (10)). When the tissue is stretched I~4>1, its behavior is defined by the porous transversely isotropic hyperelastic material model (Equation (1)). Therefore, the behavior of the disc can be described as quasi‐incompressible for tensile loading (Beek et al., 20012001; Tanaka et al., 20042004) but highly compressible under slow compressive loads (Hagandora, Chase, & Almarza, 20112011; Spilker, Nickel, & Iwasaki, 20092009).

Few experimental studies have analyzed the biomechanical response of condylar cartilage, and almost none have evaluated the temporal cartilage (shown in Table S2). This dearth of information is largely due to the complex structure of this hyaline articular cartilage, which differs to that of other forms of hyaline cartilages in the appendicular skeleton (Milam, 20032003). The distribution of collagen fibers suggests that the fibrous zone is mainly subjected to tensile load, while the mature and hypertrophic zones support compressive loads (Singh & Detamore, 20092009). The zones proximal to the subchondral bone are rich in aggrecan and poor in collagen fibers. Therefore, the mature and hypertrophic zones of condylar cartilage primarily serve to resist compressive load. On the other hand, the collagen fibers of the fibrous region are predominantly oriented in the anteroposterior direction (Singh & Detamore, 20082008) following the main sliding direction of the condyle. The condylar cartilage is likely subjected to anteroposterior shear stress (Singh & Detamore, 20092009), which results in a preferred anteroposterior collagen fiber alignment in the fibrous zone (Singh & Detamore, 20082008). This study was based on the aforementioned relationship between shear/tensile stresses and the fibrous region of the cartilage. Even though the specimens stretched by Singh and Detamore (20082008) were cartilage samples with fibrous and hyaline‐like cartilage regions, the authors assumed that the fibrous region was the only part that supported tensile stress. Consistent with this simplification, this tissue shows greater stiffness when stretched anteroposteriorly as a result of the significant contribution of the fibrous network. Furthermore, the deviation between the results obtained for single versus multiple fiber orientations was greater for the condylar cartilage than any other soft tissues of the TMJ. This can be explained by the greater differences in stiffness between the fiberless and fibrous matrices of this tissue than the other TMJ soft tissues.

Regarding the tissues that attach the TMJ disc to the surrounding tissue, although they are extremely important for coordinated movements of the joint (Willard, Arzi, & Athanasiou, 20122012), little is known about their mechanical properties. Most studies performed to date have focused on identifying the functional requirements of the lateral and posterior attachments (Liu & Herring, 20002000; Sun, Liu, & Herring, 20022002; Tanaka et al., 20022002, 20032003). Although these studies consider the composition of the fibrous tissue, they do not account for fiber orientation. However, in their study Murphy and coworkers tested the response of discal attachments to tension in two different directions (Murphy et al., 20132013). By using those data, we were able to separately characterize the contributions of the hyperelastic matrix and the fibrous network in collateral ligaments and retrodiscal tissue.

The collateral ligaments join the disc to the head of the condyle, where they combine with the joint capsule. These lateral attachments keep the disc fixed to the condyle during retrusion and protrusion movements. The fiber bundles in the collateral ligaments are mainly oriented in an anteroposterior direction to withstand anteroposterior sliding movement of the condyle, maintaining the TMJ disc positioning (Murphy et al., 20132013). To date, the few computational studies that have modelled the collateral ligament (Commisso, Martínez‐Reina, & Mayo, 20142014; Pérez del Palomar & Doblaré, 20062006) have simulated its behavior based on the known properties of the knee ligament. However, Murphy et al. (20132013) did conduct some experimental studies of porcine collateral ligament. By using those data, which were generated from samples that closely resemble the human collateral ligament, we can separately consider the properties of the hyperelastic fiberless matrix and the fibrous network.

Our numerical data corroborate this function–orientation relationship, demonstrating greater stiffness of the collateral ligament in the anteroposterior direction. Remarkably, our data indicate that the medial collateral ligament is softer than the lateral ligament. This may be due to the lower collagen fiber content of the medial region (Willard et al., 20122012). According to our results, the response of the collateral ligaments is highly dependent on the orientation of the fibrous network. Thus, the connection formed by the collateral ligaments between the TMJ disc and the condyle is less stiff when the fibers in these ligaments are oriented mediolaterally. This phenomenon may be implicated in the progression from normal to abnormal joint motion observed in patients with TMDs.

In their study of the mechanical properties of the posterior attachment, Tanaka and coworkers examined the viscoelastic response to compressive efforts of bovine retrodiscal tissue, but did not consider fibrous reinforcement (Tanaka et al., 20022002). However, a 2013 analysis (Murphy et al., 20132013) of the response of porcine retrodiscal tissue to tensile stress did consider the influence of fiber orientation. Thus, by using those data we were able to characterize the parameters of the transversely isotropic hyperelastic material model. Our results show that greater stiffness is associated with fiber orientation in the anteroposterior direction. Moreover, we found that the stiffness of the fiberless matrix was similar to that of the fibrous network, possibly because fibers are more randomly dispersed in this tissue, resulting in stiffer isotropic behavior.

### 4.1. Study limitations

The findings presented here constitute an important step toward characterizing the fibrous properties of the soft tissues of the TMJ and the porous‐fibrous properties of the TMJ disc. Nonetheless, several limitations of the study should be considered. First, for the material characterization of each sample, we assumed a preferential direction of fiber orientation. Although, the consideration of different regions and fibers orientations did not cause several differences, it is better to characterize the properties of the samples based on specific histological studies for each sample and biaxial mechanical test. Another limitation is that the experimental data derived from different species were used to test the computational models described here. Furthermore, the study from which we obtained the data used to characterize the compressive behavior of the TMJ disc (Beek et al., 20012001) did not test the medial and lateral regions of the TMJ disc. Therefore, the compression parameters characterized in this study are also applied to the anterior, intermediate, and posterior bands of the disc in part 2 of this study (without distinguishing between medial and lateral regions). Future experimental studies should test samples from the same species and separately analyze all five regions of the TMJ disc. It should be noted that we did not consider the hyaline‐like cartilage region in our characterization of the properties of condylar cartilage owing to a lack of sufficient reliable experimental data to describe permeability functions. Finally, the specific test protocols used in this study most likely influenced our findings. The establishment of standard protocols for future experimental tests would considerably facilitate the comparison of results across studies.

## 5. CONCLUSION

This study describes the characterization of the porous fibrous properties of the TMJ disc and the fibrous properties of the condylar cartilage, collateral ligaments, and retrodiscal tissue. The transversely isotropic hyperelastic material model showed limited validity for the characterization of the compressive behavior of the TMJ disc. Consequently, we evaluated the validity of a porous hyperfoam material model for this purpose. These two material models were combined, building upon an existing material model proposed by Ortún‐Terrazas et al. (20182018) and later validated (Ortún‐Terrazas et al., 20192019, 20202020). The present study further develops the approach proposed by Pérez del Palomar and Doblaré (20062006) by considering the fibrous‐porous properties of the TMJ disc, and extends this approach to the other soft tissues of the TMJ. Based on our findings, the following conclusions can be drawn:

• Transversely isotropic hyperelastic material models can be used to explain the differences in the behavior of the TMJ disc, condylar cartilage, collateral ligaments and retrodiscal tissue associated with loading direction.
• The porous component of the TMJ disc mediates reduced disc stiffness at low rates of deformation and increased stiffness at higher rates.
• Complex material models are required for characterization of the biomechanical properties of the TMJ disc. The behavior of the collagen fibers and the contribution of the interstitial fluid cannot be defined using a single material model. This study demonstrates that the biphasic nature of the TMJ disc can be defined by combining porous transversely isotropic and porous hyperfoam material models.
• Compared with other TMJ soft tissues, the assumption of a preferential orientation of fiber direction resulted in poorer fitting of experimental data for the condylar cartilage.
• Fibers play a fundamental role in connecting the TMJ attachments to the movements of the condyle‐disc system, modifying their stiffness according to load direction. Collateral ligaments are less stiff when the condyle compresses the disc against the fossa, but are stiffer, facilitating guidance of the disc, during protrusion and retrusion movements.

### Conflict of interests

The authors have no conflicts of interest related to the present work.

## ACKNOWLEDGEMENTS

This work was supported by the Spanish Ministry of Economy and Competitiveness (project DPI 2016‐79302‐R), the European Social Funds and Regional Government of Aragon (grant 2016/20) and Ibercaja‐Cai Fundation (grant IT 4/18).

## REFERENCES

1

Allen, K. D., & Athanasiou, K. A. (2006). . Tissue engineering of the TMJ disc: A review. Tissue Engineering, 12, 060509063358001 Available from: http://www.liebertonline.com/doi/abs/10.1089/ten.2006.12.ft-62

2

Angelo, D. F., Morouço, P., Alves, N., Viana, T., Santos, F., González, R., … Pinho, M. (2016). . Choosing sheep (Ovis aries) as animal model for temporomandibular joint research: Morphological, histological and biomechanical characterization of the joint disc. Morphologie, 100, , pp.223–233 Available from: http://linkinghub.elsevier.com/retrieve/pii/S1286011516300261

3

Argoubi, M., & Shirazi‐Adl, A. (1996). . Poroelastic creep response analysis of a lumbar motion segment in compression. Journal of Biomechanics, 29, , pp.1331–1339 Available from: http://linkinghub.elsevier.com/retrieve/pii/0021929096000358

4

Beek, M., Aarnts, M. P., Koolstra, J. H., Feilzer, A. J., & Van Eijden, T. M. G. J. (2001). . Dynamic properties of the human Temporomandibular joint disc. Journal of Dental Research, 80, , pp.876–880 Available from: http://jdr.sagepub.com/content/80/3/876.abstract

5

Beek, M., Koolstra, J. H., & van Eijden, T. M. G. J. (2003). . Human temporomandibular joint disc cartilage as a poroelastic material. Clinical Biomechanics, 18, , pp.69–76 Available from:http://linkinghub.elsevier.com/retrieve/pii/S0268003302001353

6

Bergomi, M., Cugnoni, J., Galli, M., Botsis, J., Belser, U. C., & Wiskott, H. W. A. (2011). . Hydro‐mechanical coupling in the periodontal ligament: A porohyperelastic finite element model. Journal of Biomechanics, 44, , pp.34–38.

7

Commisso, M. S., Martínez‐Reina, J., & Mayo, J. (2014). . A study of the temporomandibular joint during bruxism. International Journal of Oral Science, 6, , pp.116–123 Available from: http://www.nature.com/ijos/journal/v6/n2/full/ijos20144a.html/nhttp://www.nature.com/ijos/journal/v6/n2/pdf/ijos20144a.pdf

8

Detamore, M. S., & Athanasiou, K. A. (2003). . Motivation, characterization, and strategy for tissue engineering the temporomandibular joint disc. Tissue Engineering, 9, , pp.1065–1087 Available from: http://www.liebertonline.com/doi/abs/10.1089/10763270360727991

9

Detamore, M. S., & Athanasiou, K. A. (2005). . Evaluation of three growth factors for TMJ disc tissue engineering. Annals of Biomedical Engineering, 33, , pp.383–390 Available from: http://link.springer.com/10.1007/s10439-005-1741-y

10

Hagandora, C. K., Chase, T. W., & Almarza, A. J. (2011). . A comparison of the mechanical properties of the goat temporomandibular joint disc to the mandibular condylar cartilage in unconfined compression. Journal of Dental Biomechanics, 2011, 212385 Available from: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=3134092&tool=pmcentrez&rendertype=abstract

11

Holmes, M. H., & Mow, V. C. (1990). . The nonlinear characteristics of soft gels and hydrated connective tissues in ultrafiltration. Journal of Biomechanics, 23, , pp.1145–1156 Available from: https://linkinghub.elsevier.com/retrieve/pii/002192909090007P

12

Holzapfel, G. (2000). Nonlinear solid mechanics: A continuum approach for engineering. New York: Wiley.

13

Ingawalé, S., & Goswami, T. (2009). . Temporomandibular joint: Disorders, treatments, and biomechanics. Annals of Biomedical Engineering, 37, , pp.976–996 Available from: http://journals.sagepub.com/doi/10.1177/154405910808701101

14

Kuo, J., Zhang, L., Bacro, T., & Yao, H. (2010). . The region‐dependent biphasic viscoelastic properties of human temporomandibular joint discs under confined compression. Journal of Biomechanics, 43, , pp.1316–1321 Available from: , doi: 10.1016/j.jbiomech.2010.01.020

15

Liu, Z.‐J., & Herring, S. W. (2000). . Masticatory strains on osseous and ligamentous components of the temporomandibular joint in miniature pigs. Journal of Orofacial Pain, 14, , pp.265–278.

16

Milam, S. B. (2003). . Pathophysiology and epidemiology of TMJ. Journal of Musculoskeletal & Neuronal Interactions, 3, , pp.382–390.

17

Murphy, M. K., Arzi, B., Hu, J. C., & Athanasiou, K. A. (2013). . Tensile characterization of porcine temporomandibular joint disc attachments. Journal of Dental Research, 92, , pp.753–758 Available from: http://www.ncbi.nlm.nih.gov/pubmed/23783320

18

Nickel, J. C., & McLachlan, K. R. (1994). . In vitro measurement of the stress‐distribution properties of the pig temporomandibular joint disc. Archives of Oral Biology [Internet, 39, , pp.439–448 Available from: https://linkinghub.elsevier.com/retrieve/pii/0003996994901759

19

Ogden, R. W. (1972). . Large deformation isotropic elasticity: On the correlation of theory and experiment for compressible rubberlike solids. Proceedings: Mathematical, Physical and Engineering Sciences, 328, , pp.567–583 Available from: http://rspa.royalsocietypublishing.org/cgi/doi/10.1098/rspa.1972.0096

20

Ortún‐Terrazas, J., Cegoñino, J., & Pérez del Palomar, A. (2020). . In silico study of cuspid' periodontal ligament damage under parafunctional and traumatic conditions of whole‐mouth occlusions. A patient‐specific evaluation. Computer Methods and Programs in Biomedicine, 184, 105107 Available from: https://linkinghub.elsevier.com/retrieve/pii/S0169260719311599

21

Ortún‐Terrazas, J., Cegoñino, J., Santana‐Penín, U., Santana‐Mora, U., & Pérez del Palomar, A. (2018). . Approach towards the porous fibrous structure of the periodontal ligament using micro‐computerized tomography and finite element analysis. Journal of the Mechanical Behavior of Biomedical Materials, 79, , pp.135–149 Available from: , doi: 10.1016/j.jmbbm.2017.12.022

22

Ortún‐Terrazas, J., Cegoñino, J., Santana‐Penín, U., Santana‐Mora, U., & Pérez del Palomar, A. (2019). . A porous fibrous hyperelastic damage model for human periodontal ligament: Application of a micro‐computerized tomography finite element model. International Journal of Numerical Methods in Biomedical Engineering, 35, e3176 Available from: http://doi.wiley.com/10.1002/cnm.3176

23

Palla, S., Gallo, L., & Gossi, D. (2003). . Dynamic stereometry of the temporomandibular joint. Orthodontics and Craniofacial Research, 6, , pp.37–47.

24

Pérez del Palomar, A., & Doblaré, M. (2006). . 3D finite element simulation of the opening movement of the mandible in healthy and pathologic situations. Journal of Biomechanical Engineering, 128, , pp.242 Available from: http://biomechanical.asmedigitalcollection.asme.org/article.aspx?articleid=1415491

25

Singh, M., & Detamore, M. S. (2008). . Tensile properties of the mandibular condylar cartilage. Journal of Biomechanical Engineering, 130, 011009 Available from: http://biomechanical.asmedigitalcollection.asme.org/article.aspx?doi=10.1115/1.283e8062

26

Singh, M., & Detamore, M. S. (2009). . Biomechanical properties of the mandibular condylar cartilage and their relevance to the TMJ disc. Journal of Biomechanics, 42, , pp.405–417 Available from: http://linkinghub.elsevier.com/retrieve/pii/S0021929008006623

27

Spilker, R. L., Nickel, J. C., & Iwasaki, L. R. (2009). . A biphasic finite element model of in vitro plowing tests of the temporomandibular joint disc. Annals of Biomedical Engineering, 37, , pp.1152–1164 Available from: http://link.springer.com/10.1007/s10439-009-9685-2

28

Storakers, B. (1986). . On material representation and constitutive branching in finite compressible elasticity. Journal of the Mechanics and Physics of Solids, 34, , pp.125–145 Available from: http://linkinghub.elsevier.com/retrieve/pii/0022509686900335

29

Sun, Z., Liu, Z.‐J., & Herring, S. W. (2002). . Movement of temporomandibular joint tissues during mastication and passive manipulation in miniature pigs. Archives of Oral Biology, 47, , pp.293–305.

30

Systèmes, D. (2014). Abaqus analysis user's guide, version 6.14 In Dassault Systèmes. Simulia: Rhode Island, United States.

31

Tanaka, E., Del Pozo, R., Sugiyama, M., & Tanne, K. (2002). . Biomechanical response of retrodiscal tissue in the temporomandibular joint under compression. Journal of Oral and Maxillofacial Surgery [Internet, 60, , pp.546–551 Available from: http://linkinghub.elsevier.com/retrieve/pii/S0278239102844792

32

Tanaka, E., Del Pozo, R., Tanaka, M., Asai, D., Hirose, M., Iwabe, T., & Tanne, K. (2004). . Three‐dimensional finite element analysis of human temporomandibular joint with and without disc displacement during jaw opening. Medical Engineering & Physics, 26, , pp.503–511 Available from: http://linkinghub.elsevier.com/retrieve/pii/S1350453304000414

33

Tanaka, E., Hanaoka, K., Tanaka, M., Van Eijden, T., Iwabe, T., Ishino, Y., … Tanne, K. (2003). . Viscoelastic properties of bovine retrodiscal tissue under tensile stress‐relaxation. European Journal of Oral Sciences, 111, , pp.518–522 Available from: http://doi.wiley.com/10.1111/j.0909-8836.2003.00081.x

34

Willard, V. P., Arzi, B., & Athanasiou, K. A. (2012). . The attachments of the temporomandibular joint disc: A biochemical and histological investigation. Archives of Oral Biology, 57, , pp.599–606 Available from: , doi: 10.1016/j.archoralbio.2011.10.004

35

Wright, G. J., Coombs, M. C., Hepfer, R. G., Damon, B. J., Bacro, T. H., Lecholop, M. K., … Yao, H. (2016). . Tensile biomechanical properties of human temporomandibular joint disc: Effects of direction, region and sex. Journal of Biomechanics, 49, , pp.3762–3769 Available from: , doi: 10.1016/j.jbiomech.2016.09.033

36

WuTH. Soil mechanics. 1967. Allyn and Bacon, Inc: Boston, MA, United States.

37

Yeoh, O. H. (1993). . Some forms of the strain energy function for rubber. Rubber Chemistry and Technology, 66, , pp.754–771 Available from: http://rubberchemtechnol.org/doi/abs/10.5254/1.3538343

Citing articles via
https://www.researchpad.co/tools/openurl?pubtype=article&doi=10.1002/jbm.b.34558&title=Computational characterization of the porous‐fibrous behavior of the soft tissues in the temporomandibular joint&author=Javier Ortún‐Terrazas,José Cegoñino,Amaya Pérez del Palomar,&keyword=hyperfoam constitutive model,porohyperelastic behavior,soft tissues,temporomandibular joint,transversely isotropic hyperelastic constitutive model,&subject=Original Research Report,Original Research Reports,