Influence of the Temporomandibular Joint in the Estimation of Bone Density in the Mandible through a Bone Remodelling Model

The temporomandibular joint (TMJ) plays a key role in the distribution of stresses in the mandible during mastication and consequently in the distribution of bone density, due to the interconnection between both variables through bone remodelling. Two finite element models of the mandible were compared to study the influence of the redistribution of stresses produced by the joint: (1) a model without TMJ, but with simplified boundary conditions to replace the joint, as done in previousmodels; (2) a more realistic model including the articular disc and some ligaments present in the TMJ. The stresses and strains in both models were compared through the strain energy density, used in many bone remodelling models as a measure of the mechanical stimulus. An anisotropic bone remodelling model was used to simulate the behaviour of mandible bone and to estimate its density distribution. The results showed that the TMJ strongly affects the stress distribution, the mechanical stimulus, and eventually the bone density, and not only locally in the condyle, but also in the whole mandible. It is concluded that it is utterly important to include a detailed model of the TMJ to estimate more correctly the stresses in the mandible during mastication and, from them, the bone density and anisotropy distribution.


Introduction
Bone density and stresses (or strains) are intimately related to each other.That is the reason why bone remodelling models (BRMs) have been widely used to predict the density distribution in bones [1][2][3][4].Bone apparent density can be more easily estimated through computer tomography (CT).However, some BRMs can additionally estimate bone anisotropy, not accessible through CT and thus provide a closer estimation of bone elastic properties.BRMs have been traditionally used to study the proximal femur (see Doblaré and García [5], Fernandes et al. [6], among many others).In the case of the jaw, many authors have used models of a mandible section to study bone around dental implants [7,8].Only a few have estimated bone density in the whole mandible from the stresses produced by mastication loads [3,4], but the temporomandibular joint (TMJ) was not included in those models.Instead, equivalent boundary conditions were used to simulate, in a very simplistic way, the constraints imposed by the joint.This has important effects on the stress distribution of the whole mandible, as discussed in the present work.
The main component of the TMJ is the articular disc, which eases the relative movement between the condyle of the mandible and the temporal bone.The articular disc absorbs and distributes the joint reaction force over a larger contact area than a direct contact between bones would achieve and reduces the friction of the bone-on-bone contact, to prevent the damage of the articulating surfaces [9].That redistribution of loads has an effect on the stresses and strains which is not only local but affects the whole mandible.Several studies have proposed different viscoelastic models to characterize the mechanical behaviour of the disc and to analyse its effect on the masticatory tasks [10][11][12][13][14].

Mathematical Problems in Engineering
In one of those works [14], a FE model of the mandible including the TMJ was used to simulate a masticatory cycle, applying the loads exerted by the jaw opening and jaw closing muscles during the cycle.In that study, the activity of the lateral pterygoid, which is hardly accessible with electromyography and has been a subject of discussion in other computational simulations [11,12], was adjusted through inverse analysis to reproduce the movement of the jaw during mastication.That work [14] analysed the evolution of stresses and strains in the mandible during the masticatory cycle, which are used in the present work to study the remodelling response in the mandible bone.
The main objective of the present work is to study how the change in stress distributions affects the remodelling response of the mandible bone.Particularly, the study shows the important effect of the TMJ on the stress distribution during the masticatory cycle.This will be done by comparing the bone density distribution estimated with two models: with and without TMJ.In the last case, the joint constraints were modelled with simplified boundary conditions, restraining certain movements of the condyle, as done in a previous study [3].In contrast, the new model includes the TMJ and simulates the contact between the articular disc and the articulating surfaces, as well as some ligaments that constrain the relative movements between the jaw and the temporal bone.

Finite Element Model.
A finite element (FE) model of the mandible, called M1 and built in a previous work [3], was compared with a more recent model, named M2 and described in detail in Commisso et al. [14].Both were built using Abaqus FEA5.The difference between them is that M2 includes a detailed model of the TMJ: ligaments, articular disc, and temporal bone.By contrast, in M1 certain displacements were constrained on the condylar surface to simulate the conditions imposed by the TMJ [3].The inclusion of the TMJ entailed some improvements.First, the disc-condyle and disc-temporal bone contacts were modelled, resulting in a better approximation of the local boundary conditions.Second, the implementation of dynamic loads with M1 is not easy, for example, to simulate a masticatory cycle in which the boundary conditions vary with time.On the contrary, M2 allows the simulation of a masticatory cycle in a straightforward manner.
Both M1 and M2 model comprise 77,490 eight-node linear brick elements (named C3D8 in Abaqus FEA element library) for the mandible.In addition, M2 model comprises 965 elements for each articular disc and 1,788 elements for the ligaments of each side, all of them of type C3D8H (eightnode linear brick, hybrid with constant pressure).In M2, the temporal bones were modelled as rigid surfaces, using 1,324 type M3D4 elements for each side.

Material Models.
The articular disc was modelled with a quasi-linear viscoelastic (QLV) model, where the uniaxial stress response to a step stretch is factorized as thus, separating the dependence of time and deformation [15]. is the reduced relaxation function defined with a 4term Prony series: The elastic response function,   , is the instantaneous stress produced by a uniaxial step stretch.The strain energy function, Ψ, proposed by Humphrey and Yin [16] was used to define   .
where , , and  are material constants and  1 is the first invariant of the left Cauchy-Green tensor.The first term was added to model quasi-incompressibility, so that  must be chosen small enough.This function was previously applied to the articular disc [17] and provides the following elastic response function in the incompressible case: The QLV model was easily implemented in Abaqus FEA by defining a hyperelastic behaviour given by (3) and using a viscoelastic behaviour with a time domain definition given by the Prony series (2).The following constants were taken from a previous experimental work [17]:  1 = 0.28,  2 = 0.37,  3 = 0.27,  4 = 0.08;  1 = 0.01 ,  2 = 0.1 ,  3 = 1.0   4 = 10 ;  = 0.16MPa,  = 4.18, and  = 0.01 MPa −1 .The ligaments of the TMJ were also modelled as in [13], following the approach of Gardiner and Weiss [18].A friction coefficient of  = 0.015 was assumed for the contact between the disc and the bony parts [19].

Bone Remodelling Models.
The BRM used here to simulate the behaviour of the mandible bone (termed later as ABRM) was proposed by Doblaré and García [5] and is an extension of the isotropic model developed by Beaupré et al. [1] to include anisotropy.It relates density, anisotropy, and mechanical properties of bone with the loads the tissue is daily subjected to.The basic aspects of the models are explained next: first the isotropic BRM [1] and next the anisotropic model [5], focusing on the modifications made to the isotropic BRM.However, consulting the original papers is advised for further comprehension of both.

Isotropic Bone Remodelling Model (IBRM).
In the IBRM Beaupré et al. [1] defined a daily mechanical stimulus based on the strain energy density (SED) accumulated in each point of the tissue during the daily activity.This stimulus considers the contribution of each load  of the normal activity.At the continuum level, it is given by where  is an empirical constant, adjusted to 4 by Whalen and Carter [20],   is the daily number of cycles of load , and   is an effective stress at the continuum level, defined as   = √2  , where  is Young's modulus and   is the SED corresponding to load  at a given point.The stimulus that controls the bone remodelling response is the daily stress stimulus measured at the tissue level,   , which is related to the stimulus at the continuum level, , through porosity  using the relation: where  is the bone apparent density and ρ is the density of bone matrix, which is assumed to be fully mineralized in this model and, thus, ρ = 2 / 3 constant.The remodelling response is measured in terms of the bone resorption/apposition rate, ṙ . This rate provides the volume of bone matrix resorbed/formed per day and per unit surface available for bone resorption/formation.It is given by the mechanostat theory [21], as a function of   and a reference stimulus,  *  , close to which no net remodelling response is observed (see Figure 1).Bone density change rate is where the specific surface,  V , is the free bone surface (where bone resorption/formation occurs) per unit volume and was correlated with porosity by Martin [22].Finally, Young's modulus and Poisson's ratio are related to the apparent density by the following correlations [23]:

Anisotropic Bone Remodelling Model (ABRM).
This ABRM, proposed by Doblaré and García [5], is an extension of the IBRM to the anisotropic case.The anisotropy is measured with the fabric tensor Ĥ, normalized such that det ( Ĥ) = 1.Then, a tensor H is defined to consider jointly the porosity and the orientation of the pores, which is measured by Ĥ: where () and () represent, respectively, the constant and exponent in the correlation (8a), which depend on  (e.g.( ρ) = 3.2, ( ρ) = 1763).
The ABRM establishes the strain tensor, , as the mechanical variable that drives bone remodelling.The tensorial mechanical stimulus is given by with Ĝ and λ the Lamé constants corresponding to bone with density ρ and obtained from (8a) and (8b).To weigh the relative influence of the spherical and deviatoric parts of the stimulus, a new stimulus tensor, J, is defined as where the parameter  ∈ [0, 1] must be chosen a priori.If  = 0, the model is purely isotropic and if  = 1, J = V(Y) and the spherical part of the stimulus has no influence on the bone response.Doblaré and García [5] defined two functions,   and   , to establish the resorption and formation criteria, respectively.These functions were dependent on the stimulus tensor, J, the reference stimulus,  *  , and the width of the dead zone, 2w (see Figure 1), and define the domains of the stimulus J for which formation, resorption, or no net remodelling response (dead zone) take place, analogously to the domains defined in Figure 1 for   : The remodelling response ṙ is given by the active remodelling criterion, that is, the condition in ( 12) which is currently accomplished: in formation 0 in the dead zone −     2−()/2 in resorption (13) where   and   are, respectively, the slopes of the resorption and formation ramps (Figure 1).Finally, this remodelling response, ṙ , is used to calculate Ḣ, which leads to the variation of density and anisotropy.

Modified Anisotropic Bone Remodelling Model (MABRM).
ABRM was modified by Ojeda [24] to consider cyclic loads more appropriately.In cyclic loads like those shown in Figure 2 the remodelling response would depend, according to Carter et al. [25], on the typical peak of stimulus and the number of cycles performed during the daily activity.
As far as the authors know, all the BRMs proposed in the literature treat cyclic loads in a simplified way: by solely considering one instant of the cycle, usually when the applied forces are maximum.Thus, a pseudo-static elastic problem is solved at this instant to calculate the stresses, the mechanical stimulus, and the remodelling response in every point of the bone.However, in a chewing cycle, the peaks of stimulus can be out of phase; that is, each point of the mandible can reach its peak of stimulus at a different instant of the chewing cycle.That way, the simulation of a single instant would overlook the peaks of stimulus in some points of the bone.For instance, in the previous study implementing M1 [3] it was assumed that the whole mandible was subjected to the higher stresses (and stimulus) at the instant of centric occlusion (CO); but that might not be the instant of maximum stimulus for the area near the mandibular foramen, where the depressor muscles are attached.The stresses in this area are relatively low at CO, but they are higher when the mandible is being pulled by the depressor muscles to open the mouth.Cases like that must be analysed with MABRM that considers appropriately the variation of stimulus through the cycle.
Let us consider three individuals (I1, I2, and I3) whose daily activities produce evolutions of the stimulus like those shown in Figure 2. I1 (in blue) carries out an intense exercise inducing net bone formation.I3 (in green) performs an activity of low intensity leading to net bone resorption.Finally, I2 (in red) performs an activity of moderate intensity which does not produce a net change of bone mass (dead zone).In cyclic loads peaks of stimulus coincide in time with maximums of the function   , named here  Finally, the bone resorption/apposition rate is obtained as in (13): and Ḣ is evaluated using (14).

Loads.
A complete mastication cycle with the right molars (RM) was simulated.In M2, the viscoelastic nature of the articular disc could cause differences between the first and subsequent cycles in a mastication sequence.However, the fast stress relaxation occurring in the disc allowed analysing only the first cycle as representative of the whole sequence, as concluded in [14].The model is symmetric and therefore the stresses produced by mastication with the left molars were obtained from RM by applying a symmetry operation.The closing and opening phases of the cycle were simulated by applying the muscle activation patterns given by Hylander [26], except for the lateral pterygoid.The activity of this muscle is difficult to measure and has been widely debated in the literature, with no agreement in the activation pattern the muscle has during the mastication cycle.In a recent work, we have adjusted that activity with an inverse procedure [14].This activation pattern was used here and can be seen in Figure 5, along with the activities of the rest of masticatory muscles.These activities were multiplied by the maximum forces of each muscle, F  0 , the peak of activity,   , and the direction cosines shown in Table 1, to get the components of the muscular forces.The muscle forces were imposed as external loads, distributed over the insertion area of each muscle as can be seen in Figure 3 for the closing muscles and Figure 4 for the opening muscles.In the figures, the areas highlighted in different colours correspond to various groups of nodes where the different muscles are inserted.Also, an arrow indicates the approximate orientation of each muscle.
Those direction cosines give the orientation of muscle fibres (in average as some muscles have a fan shape) relative to the mandible and corresponding to the instant of centric occlusion (CO).However, that relative orientation may vary with mouth opening and needed to be updated throughout the cycle, in the following way.For each muscle the origin and insertion were distinguished, the latter being its attachment to the mandible.First, the origin of each muscle was estimated using the position of the insertion point in the mandible, the direction cosines at CO, and the length shown in Table 1.The origins are located in the skull and were thus assumed fixed during the simulation, so that the updated direction cosines were easily calculated from the position of those fixed origins and the position of the mandible (more details can be found in [14]).
The instant of CO is indicated in Figure 5 with a red vertical arrow (around  = 0.19 ).This is approximately the instant of maximum mastication force and was the only instant analysed with M1 in the previous work of Reina et al. [3].As stated above, M1 was not adequate to simulate a mastication cycle, given that the displacement constraints used in that model were only valid for the instant of CO.Therefore, M1 simulated the masticatory activity as a static load with the forces exerted by the closing muscles applied in the highlighted areas of Figure 3 and their corresponding magnitudes given in Figure 5 at CO.In contrast, M2 was able to simulate a mastication cycle with the complete record of muscular activity.Inertial forces were very small and thus neglected, leading to a pseudostatic simulation of the cycle.

Displacement Boundary Conditions.
The boundary conditions applied to simulate CO in each model are depicted in Figure 6.At CO, the mouth is closed and the condyles are at their back position, in contact with the articular eminence of the temporal bone.In M1 this contact was not considered as such, but it was simplified by constraining the displacements of the articular surface of the condyle.The pressure in this contact is different for each side, being usually smaller in the ipsilateral (working) side [30,31], where food thickness interposes between both dental arches, thus separating the Table 1: Forces exerted by masticatory muscles.Only the orientation of the forces of the right side are given (symmetry with respect to the sagittal plane was assumed).The magnitude of the muscle force peak in CO for the ipsilateral (I) and contralateral (C) sides are also provided.The +x-axis is directed anteriorly, the +y-axis rightward, and the +z-axis downward.CP: closing phase, OP: opening phase.1: length with the mouth closed [27].2: Korioth et al. [28] for the jaw closing muscles.3: van Eijden et al. [27] for the jaw opening muscles.4: Nelson [29] for the jaw closing muscles.5: proportional to the values adopted by Hannam et al. [12] for the jaw opening muscles.Muscle  articular surfaces.So, for instance, in the simulation of RM with M1 the articular surface of the left condyle was fixed and the right condyle was assumed to move freely (Figure 6).It must be noted that this is a strong simplification of the real situation in which the ipsilateral joint bears a lower reaction force, though not null [30,31].
The new models aimed to avoid this simplification by modelling the contact at the joints.In M2, temporal bones were fixed, condyles were constrained by the contact interaction with the articular disc, and, in turn, articular discs were constrained by the contact interaction with temporal bones.This way reaction forces arose naturally from the contact at the joints, with a lower value in the ipsilateral condyle.Additionally, the ligaments of the joint limit the movement of the mandible by preventing condyles from being pulled apart from the articular eminence.Figure 7 shows the finite element model of the collateral and temporomandibular ligaments and the posterior part of the articular capsule.Additionally, the direction of the collagen fibers are indicated with arrows (more details on the joint model can be found in Commisso et al. [14]).
On the other hand, the tooth-food contact was not modelled.Instead, the vertical displacements were constrained in the occlusal face of ipsilateral molars (e.g., right molars in RM) in order to simulate grinding forces as reaction forces [3].2, were performed.1 was the same performed in a previous study [3], using model M1 (without TMJ), in a static problem and applying the muscle forces at the instant of CO (red vertical arrow in Figure 5).Simulation 2 −  uses the model M2, including the TMJ, to perform a pseudostatic analysis of the instant of CO.Since the articular disc has a viscoelastic behaviour, the time variable was important, but the analysis 2 −  was simplified to ignore the temporal evolution of muscle forces.Instead, these muscle forces were varied from 0 to the values at CO in a ramp-like manner of length  = 0.2, small enough to focus just on the instant of CO, but long enough to be considered a pseudostatic analysis and not to introduce spurious stresses due to the application of instantaneous loads in viscoelastic models.

Description of the Performed Simulations. Three simulations, summarized in Table
Simulation 2 −  uses the model M2 as well, but now to perform a pseudostatic analysis of the whole masticatory cycle.Muscle forces were varied following the pattern depicted in Figure 5.
The effect of the simplification of the loads could be assessed through comparison of the stresses obtained in 2 −  and those obtained in 2 −  at CO.The differences were negligible and, thus, the simplification was validated.Nonetheless, the interest of 2 −  mostly lies in simulating the whole masticatory cycle and specifically the opening phase during which some points of the mandible could reach its peak of stimulus.
As stated before, the vertical displacements of the occlusal face of ipsilateral molars were constrained to simulate grinding forces as reaction forces.This was made in simulations 1 and 2 − , while in 2 −  vertical displacements were constrained only during the closing phase and released during the opening phase.A summary of the performed simulations is given in Table 2.

Simulation of the Bone Remodelling
Response.In the previous work [3] (simulation 1) and in other similar works, the simulations started from a bone with uniform density ( 0 = 0.5 / 3 in that particular case) and initially isotropic.The final density and anisotropy distribution were estimated by simulating the daily masticatory activity until a remodelling equilibrium was achieved, with no further changes in bone density.(It could also be checked that the convergence of density implied the convergence of anisotropy.) In simulations with model M2 it was not advisable to start from a low uniform density ( 0 = 0.5 / 3 ).In such case, numerical problems arise from the disc-condyle contact, due to the low stiffness of both materials.Moreover, the computing cost of simulating the large number of days needed to achieve the remodelling equilibrium is excessive in the more complex model M2.For these reasons, simulations 2 −  and 2 −  started from the bone density and anisotropy distribution obtained in simulation 1 at the remodelling equilibrium.The simulation of the masticatory cycle with M2 produced a redistribution of stresses, caused by the TMJ, thus leading to a change in the density distribution until a new remodelling equilibrium was achieved.The density distribution obtained at the new equilibrium was compared with that obtained in 1 to analyse the effect of the TMJ.
The masticatory pattern assumed here was an alternating unilateral molar chewing followed by 75% of the population [32]: chewing with the right molars (RM) followed by chewing with the left molars (LM), thus, defining the sequence: RM, LM, RM, LM,. .., until the remodelling equilibrium was achieved.A total of   = 500 (see ( 5)) daily masticatory cycles (250 RM + 250 LM) were assumed, like in Reina et al. [3].The convergence of the density distribution was checked by computing the following parameter: where   represents the density of a given point at day .The remodelling equilibrium was assumed to occur when  < 1%.

Comparison of Simulations
1 and 2 − .The strain energy density (SED) distribution obtained in 1 and 2 −  at the remodelling equilibrium is compared in Figure 8.A detail of this distribution in the ramus and the body of the mandible (cross-sections 1-1' and 2-2') can be seen in Figures

8(c) and 8(d)
. SED is shown because it is closely related to bone density distribution, given that the mechanical stimulus is defined in terms of SED.
Starting from a uniform density distribution ( 0 = 0.5 / 3 ) and an initially isotropic bone, around 300 days of mastication were needed to achieve the remodelling equilibrium state in 1.Starting from that state, SED changed noticeably in 2 − , inducing changes in density distribution.For that reason, 10 further days of mastication were needed to achieve a new remodelling equilibrium state ( < 1%).The new density distribution is compared with that of 1 in Figure 9.A detail of both density distributions at a premolar section is compared with a CT scan of the real mandible in Figure 10.
Mastication forces at the instant of CO could be estimated from reaction forces at the nodes of the occlusal faces of the first and second right molars, where the displacements were constrained.The maximum mastication force was 501  in 1 and a little lower, 458 , in 2− .

Comparison of Simulations 𝑀2 − 𝐶𝑂 and 𝑀2 − 𝑐𝑦𝑐𝑙𝑒.
Comparison of simulations 2 −  and 2 −  allows analysing the effect of MABRM and to check whether the variation of the mechanical stimulus through the masticatory cycle affects the remodelling response of the mandible.The difference of bone density between both simulations was hardly noticeable and, therefore, instead of showing  its distribution for 2 − , the following variable was represented in Figure 11: where  2− and  2− represent the density of a certain point at the remodelling equilibrium in simulations 2 −  and 2 − , respectively.

Discussion
The comparison of simulations 1 and 2 −  revealed significant differences in all aspects: SED, bone density distribution and mastication forces.SED is more uniformly distributed in 2− and not so concentrated in the ipsilateral molar region, as it was in 1 (Figure 8).The cross-sections in Figures 8(c) and 8(d) also revealed great differences.In the ramus (plane 1-1'), 2 −  estimates a continuous layer of high SED surrounding the  central part of the section, not seen in 1.In the molar region (plane 2-2'), the inner area of low SED is larger in 2 −  and, again, there is a continuous layer of high SED enclosing the section.This is not seen in 1, due to a small portion of the inferior medial surface which has a very low SED.Given that SED is directly related to the mechanical stimulus, the differences of SED can explain the different estimated density.A significant change is seen at the condyles, where bone density is more uniformly distributed in 2 − , due to the TMJ and the redistribution of loads it produces.Constraining displacements, as done in 1, always produce a stress concentration in the FE solution, which leads to a discontinuous bone density across elements (see detail of the condyle in Figure 9(a)).In the rest of the mandible, the TMJ induces smaller changes, but they are still worth a mention, for example, in the chin.Here, the area of density close to 1 / 3 (green in Figure 9) spreads in 2−, reducing the area of maximum density (2 / 3 ), but, more importantly, reducing the area of minimum density (0.3 / 3 ), in light blue.In conclusion, the simulations that use M2 tend to close the so-called cortical shell, a continuous layer of intermediate to dense bone that surrounds the real mandible.
If numerical results are compared with CT scan at the premolar region (Figure 10), bone density distribution obtained with 2 −  is, again, closer to the actual one.It is true that 2 −  produces less cortical bone of high density (in red) than 1, but, it corrects an important problem seen in the old model: the cortical shell was missing at the inferior side of the canine region, where a portion of very light bone (∼ 0.5 / 3 ) was obtained.This tubular structure of dense bone surrounding an inner central area of bone of low density is usual in the diaphysis of long bones and is also seen in the mandible, to resist bending and torsion loads that mastication produces [3].
Mastication force was obtained as the resultant reaction force in the occlusal face of molars, in both 1 and 2−.Despite the applied muscle forces being the same, mastication force was almost 10% higher in 1 (501 ) than in 2 −  (458 ).Both values were within the range measured in experimental studies, 430 − 650  [33,34].
The difference found between 1 and 2 −  shows that the TMJ, modelled in M2, does not only redistribute the stresses in the condyle, but also reduce the stresses in the alveolar region around molars.The same conclusion can be drawn by comparing the distribution of SED (Figure 8).With TMJ (M2) SED is more uniformly distributed, and not only in the condyles, as could be expected, but also in the ramus and the body of the mandible.In the model without TMJ (M1) SED was mainly concentrated around the molar region of the ipsilateral side.
In this regard, the flexibility of the articular disc might plays an important role to redistribute the stresses, as well as the stress relaxation that the disc experiences, which increases that flexibility.For that reason, in 2− it was important to apply loads at a rate close to the real one.In this sense, a ramp of length  = 0.2s seemed adequate as it produced very similar results compared to the more realistic simulation 2 − .
2 −  also improves the estimation of density in the alveolar region right below teeth, which is clearly noticeable in the premolar section (Figure 10).1 produces a bone of intermediate density in that area, unlike 2 − , which predicts the actual light bone that can be seen in the CT scan.
The influence of analysing only the instant of CO or the whole masticatory cycle was studied through the comparison of the bone density distributions obtained in 2 −  and 2 − .This comparison revealed very little differences (Figure 11), which were explained by the low mechanical stimulus the mandible was bearing during the opening phase.This stimulus was much lower than that produced during the closing phase for two reasons: (1) jaw opening muscles exert forces of lower magnitude than jaw closing muscles and (2) during the closing phase the displacements were constrained in the molars, giving rise to the mastication force.The mechanical stimulus was particularly high at CO, when the activity of closing muscles is maximum.Consequently, most of the mandible reaches its peak of stimulus at CO and bone formation is driven in simulation 2 −  by the stress state at that instant, which is the one captured in 2 − .There was a very small area in the lingual side of the chin, close to the insertion of opening muscles, where the peaks of stimulus were reached during the opening phase.Nonetheless, these peaks were similar to those obtained at CO and they did not lead to significative differences in the estimated bone density.In conclusion, the simulation of the whole mastication cycle did not yield relevant differences with respect to previous analysis in which only the instant of CO was studied.In other words, the simplification of the analysis to focus just on CO was justified from the perspective of the remodelling response of the mandible.This fact does not mean that the modification introduced in MABRM to study cyclic loads is irrelevant on a general basis.It was in the case of mastication, but not in the case of the femur during gait, for example, where the peaks of stimulus were reached in different phases of the cycle for different locations in the femur [24].In this particular case, it was shown that MABRM led to important differences from the bone remodelling perspective.
Some limitations of the model are discussed next.First, the articular disc was assumed as isotropic.Actually, it is an anisotropic material, composed of an extracellular matrix reinforced with collagen fibres.These fibres run in anteroposterior direction in the central portion of the disc and in mediolateral direction, in both the anterior and posterior regions [35].In any case, fibres run mostly perpendicular to the thickness.In the load cases simulated in this work the disc was mainly subjected to compression in vertical direction (across its thickness).With such loads, fibres would be stretched, so contributing to the stiffness of the tissue.For that reason, they should be taken into account in the disc constitutive model, but they were not.Nonetheless, compression along thickness was the type of load applied in the experiments where the constitutive model was adjusted [13].Therefore, the stiffness of the fibres was indirectly included in the material model.In any case, the anisotropy should be explicitly considered (by using a fibre reinforced material model) if the type of load varied during the mastication cycle.
Another limitation of the study was the omission of the articular cartilage that covers the surface of articulating bones.The articular cartilage helps to reduce the stresses in the condyle [36].The mechanical properties of this layer has not been established appropriately, but it is generally accepted that it is more flexible than the articular disc [10].For this reason, it is thought to have an effect on the redistribution of stresses over the condyle that needs to be addressed in future works.Other limitations are the simplified definition of the teeth-food contact and not considering the pennation angle of muscles, though this was simplified by using the average direction of muscle fibres and taking into account that pennation angle is close to zero in all the muscles involved.
Overcoming these limitations is the aim of the short-term future research, but another interesting topic to be addressed is to analyse the influence of pathologic conditions of the articular disc in the reaction forces and stress distribution of the mandible.It is known that bruxism may cause damage on the disc [37], likely altering its viscoelastic properties.Since the influence of the TMJ has been shown to be important on the stress distribution, this alteration of viscoelastic properties could, in turn, affect the bone density distribution within the mandible, specially if damage is not symmetrical.As a matter of fact, another interesting topic, not yet addressed in the literature, is to study the effect of other types of mastication patterns, such as unilateral, on bone density distribution.

Conclusions
A FE model of the mandible including the TMJ was used to simulate the mastication process.The most relevant features of the joint were modelled.A viscohyperelastic model fitted from experimental tests was used to simulate the behaviour of the TMJ articular disc.In contrast to previous similar studies, mandible bone was considered a deformable material, which, additionally, is able to remodel and change its mechanical properties to adapt itself to the mechanical environment.The model presented here allows analysing the influence of the joint on the stresses produced in bone during the masticatory cycle and the effect of those stresses on the bone remodelling response.
By comparing the FE models with and without TMJ it was concluded that the joint plays an important role by redistributing the stresses and strains and not only over the condyles but also throughout the whole mandible.
Regarding bone remodelling, it was observed that the joint has also an important effect on the bone density distribution predicted by the model.If the joint was modelled, the estimated density was closer to the actual one, specially below the teeth and in the chin.Here, a cortical shell confining a central area of trabecular bone was seen in the whole body of the mandible, unlike the model without TMJ, which failed to close the cortical shell in the premolar and incisive region.
Furthermore, considering the complete mastication cycle instead of only the instant of centric occlusion had a small effect in the estimation of density and could be disregarded.This is an important conclusion, as it simplifies this type of studies of bone remodelling in the mandible.In other words, the analysis of the whole cycle does not seem necessary in this case from the remodelling perspective, as it was in other cyclic loads, such as gait cycle, where it had a strong influence on femoral bone remodelling [24].

Figure 1 :
Figure 1: Remodelling response as a function of the daily stress stimulus at the tissue level.

Figure 2 :
Figure 2: The evolution of stimulus in three individuals are represented in solid line.The peaks correspond to the maximums of   and minimums of   .Dashed lines: alternative activities with the same peaks but different valleys.
, and with minimums of   , named    .(They only coincide in time, because    and    have different values, in general.)Similarly, valleys of stimulus coincide in time with minimums of the formation function,    , and with maximums of the resorption function,   .According to Carter et al.[25] the remodelling response depends only on the peaks of stimulus.Thus, the activities plotted in dashed lines for I1 and I2 would lead to the same remodelling response as those plotted in solid lines.In MABRM the whole cycle is simulated to calculate the evolution of   and   .With those evolutions,    and    are calculated in every point of the model to apply the remodelling criterion (15) that places the peaks in one of the three regions: formation, resorption, or dead zone.

Figure 3 :Figure 4 :
Figure 3: Insertion and schematic orientation of different portions of the closing muscles used in models M1 and M2.

Figure 5 :
Figure 5: Activation pattern of the jaw closing and opening muscles during unilateral mastication.The ipsilateral (I) muscles are above, in the upper y-axis, while the contralateral (C) ones are represented in the inferior y-axis.The amplitudes are normalized to the corresponding peak of muscle force, given in Commisso et al. [14].The closing and opening phases are named after those muscles which are active during the corresponding phase (jaw opening or jaw closing muscles) and not after the direction of the movement, which is shown down in the figure.The instant of CO is indicated with a red vertical arrow.

Figure 6 :
Figure 6: Boundary conditions applied at the instant of CO in a mastication with the right molars: (a) in M1 and (b) in M2.The triangle represents a rigid fixation of temporal bone (contact pairs temporal bone-disc and disc-condyle are defined in the TMJ), while the arrows indicate the direction in which the corresponding displacements are constrained.

Figure 7 :
Figure 7: Finite element model of the articular disc and ligaments in the TMJ.The direction of the collagen fibers in the ligaments is indicated with arrows.

5 Figure 8 :
Figure 8: Distribution of SED at the instant of CO in 1 and 2 − , for the whole mandible (a, b) and at the cross-sections 1-1' and 2-2'(c, d).Plane 1-1 cuts the ramus and plane 2-2 is a cross-section located approximately at the second molar.