Przemysław Stróżyk,Jacek Bałchanowski
Application of numerical simulation studies to determine dynamic loads acting on the human masticatory system during unilateral chewing of selected foods 10.3389/fbioe.2023.993274
dynamic patterns of a muscles,dynamic characteristics of foods,numerical simulation,kinematic-dynamic model of unilateral chewing,experimental tests of food
Introduction: This paper presents its kinematic-dynamic computational model (3D) used for numerical simulations of the unilateral chewing of selected foods. The model consists of two temporomandibular joints, a mandible, and mandibular elevator muscles (the masseter, medial pterygoid, and temporalis muscles). The model load is the food characteristic (i), in the form of the function Fi = f(Δhi )−force (Fi ) vs change in specimen height (Δhi ). Functions were developed based on experimental tests in which five food products were tested (60 specimens per product). Methods: The numerical calculations aimed to determine: dynamic muscle patterns, maximum muscle force, total muscle contraction, muscle contraction corresponding to maximum force, muscle stiffness and intrinsic strength. The values of the parameters above were determined according to the mechanical properties of the food and according to the working and non-working sides. Results and Discussion: Based on the numerical simulations carried out, it can be concluded that: (1) muscle force patterns and maximum muscle forces depend on the food and, in addition, the values of maximum muscle forces on the non-working side are 14% lower than on the working side, irrespective of the muscle and the food; (2) the value of total muscle contraction on the working side is 17% lower than on the non-working side; (3) total muscle contraction depends on the initial height of the food; (4) muscle stiffness and intrinsic strength depend on the texture of the food, the muscle and the side analysed, i.e., the working and non-working sides.
Application of numerical simulation studies to determine dynamic loads acting on the human masticatory system during unilateral chewing of selected foods 10.3389/fbioe.2023.993274

Introduction: This paper presents its kinematic-dynamic computational model (3D) used for numerical simulations of the unilateral chewing of selected foods. The model consists of two temporomandibular joints, a mandible, and mandibular elevator muscles (the masseter, medial pterygoid, and temporalis muscles). The model load is the food characteristic (i), in the form of the function Fi = f(Δhi )−force (Fi ) vs change in specimen height (Δhi ). Functions were developed based on experimental tests in which five food products were tested (60 specimens per product). Methods: The numerical calculations aimed to determine: dynamic muscle patterns, maximum muscle force, total muscle contraction, muscle contraction corresponding to maximum force, muscle stiffness and intrinsic strength. The values of the parameters above were determined according to the mechanical properties of the food and according to the working and non-working sides. Results and Discussion: Based on the numerical simulations carried out, it can be concluded that: (1) muscle force patterns and maximum muscle forces depend on the food and, in addition, the values of maximum muscle forces on the non-working side are 14% lower than on the working side, irrespective of the muscle and the food; (2) the value of total muscle contraction on the working side is 17% lower than on the non-working side; (3) total muscle contraction depends on the initial height of the food; (4) muscle stiffness and intrinsic strength depend on the texture of the food, the muscle and the side analysed, i.e., the working and non-working sides.

Food consumption is one of the most important activities (Lund, 1991), necessary to sustain life processes at the same time as a complex kinematic-dynamic process (Daumas et al., 2005; Hedjazi et al., 2013; Peck et al., 2000; Piancino et al., 2008; Slager et al., 1997; Stróżyk and Balchanowski, 2018), which is controlled by the central nervous system (Dellow and Lund, 1971; Lund, 1991). Mastication is also a multi-parametric issue encompassing many complex and synchronised processes for preparing a bolus for swallowing. The most important processes may include: 1) the dynamic processing of the food (Stokes et al., 2013) dependent on its position on the dental arch (Manns and Díaz, 1988), 2) the continuous change in the geometric dimensions and mechanical properties of the food (Lenfant et al., 2009; Lillford, 2000; Lucas et al., 2004; Mioche et al., 2002; Stokes et al., 2013) and 3) the performance of complex movements by the mandible during chewing (Hylander, 2006; Koolstra, 2002; Piancino et al., 2012; Posselt, 1952; Quintero et al., 2013; Slavicek, 2010; Weinberg, 1963). There are many papers in the literature where the authors (Çakir et al., 2011; Hutchings et al., 2012; Kohyama et al., 2007; Kohyama et al., 2008; Santana-Mora et al., 2014; Wang et al., 2010) focus on various parameters related to the act of chewing, e.g., 1) change in food texture, 2) chewing effort, 3) muscle activity, 4) some cycles needed to prepare a bolus of food, 5) degree of wetting of food with saliva, 6) age, 7) pathogenic changes, but most importantly 8) loads that occur during chewing of food. On the other hand, based on the data reported by (Agrawal et al., 1998; Hiiemae et al., 1996; Koolstra, 2002; Mathevon et al., 1995; Shimada et al., 2012), it appears that the mechanical properties of food are important for the function of the masticatory system during symmetric incisal biting and chewing. Furthermore, the work results (Mioche et al., 1999; Stróżyk and Bałchanowski, 2018) indicate that food imposes individual patterns of muscular force that must adapt to different functional requirements (Fitts et al., 1991). Their values are controlled by the central nervous system (Dellow and Lund, 1971; Lund, 1991), based on remembered external stimuli (patterns) perceived by different receptors (Quintero et al., 2013) so that it is adapted to the food consumed (Hiiemae et al., 1996; Knudson, 2007; Murray, 2004). One of the main ways to determine muscle forces is through electromyography (EMG)-based method and mathematical equations (Ferrario and Sforza, 1996; Itoh et al., 1997; Mioche et al., 1999; Murray et al., 1999; Pruim et al., 1980; Weijs and Hillen, 1985). The method is mainly applied to the masseter muscle and the temporalis muscle. In the case of the medial and lateral pterygoid muscles, measurements are more difficult because they cannot be externally accessed (Koole et al., 1990; Murray et al., 1999; Widmalm et al., 1987; Wood et al., 1986). Consequently, determining muscle forces requires intra-oral access, which makes measurements during chewing difficult due to the possibility of interference with the natural act of chewing. Since using EMG during chewing makes it impossible to determine the forces in all the mandibular elevator muscles, it is a good option to use a numerical computational model of the masticatory system. The model developed for numerical calculations provides reproducibility but requires the preparation of basic data, e.g., the geometry of the mandible, the muscle attachment sites, the muscle model and the mode of support. As chewing is a dynamic issue, a model should be developed based on solid mechanics, i.e., a kinematic-dynamic model (Stróżyk and Bałchanowski, 2016; Stróżyk and Bałchanowski, 2018; Stróżyk and Bałchanowski, 2020). The load of the model should be the characteristics of the food in the form of the function F i = f(Δh i )—force (F i ) vs. change in specimen height (Δh i ), determined during experimental tests, while the movement of the mandible is set in the path of mastication determined for the tested food. Preliminary analysis of the function of the masticatory system in terms of: 1) the position of the food (bite force) on the dental arch, 2) the number of cycles and 3) the muscle activity on the left and right sides of the mandible indicate that, in mechanical terms, unilateral chewing will be an interesting case of physiological load on the masticatory system. In the literature, it is possible to find publications in which the authors model unilateral chewing, but these are static models in which the position of the mandible corresponds to the closure of the mouth (Korioth et al., 1992; Pinheiro and Alves, 2015; Pachnicz and Stróżyk, 2021; Reina et al., 2007). Interesting information on unilateral chewing can be found in (Harrison et al., 2014). However, the authors focus on dynamic changes in the mechanical properties of the food rather than on changes in the values of the muscular forces acting on the masticatory system. Analysing other publications on chewing (Choi et al., 2005; Jahadakbar et al., 2016; Lee et al., 2017; Luo et al., 2017; Marková and Gallo, 2016; Reina et al., 2007), we find that none of them present changes in muscle force values (muscle functioning patterns) to dynamic changes occurring in the food during mastication. The primary aim of this study was to determine computational model of unilateral chewing, loaded with the dynamic patterns of food (i) determined in experimental tests. The model prepared in this way made it possible to determine, firstly, selected dynamic parameters of the mandibular elevator muscles, i.e., total muscle contraction (q i ) and dynamic muscle forces (F i ) as functions, q i = qi(t) and F i = F i (t), respectively. The parameters obtained allowed the determination of muscle dynamic patterns in the function F i = f(q i ). From the results obtained, parameters such as muscle stiffness (K i ) and intrinsic strength (k i ) were also determined as a function of food. As mastication is a complex mechanical process (Stokes et al., 2013), the study was limited to first cycle. The results demonstrate the feasibility of using a hybrid model (based on experimental tests of food and numerical simulation of unilateral chewing) to show how the dynamic parameters of the mastication system change as a function of the mechanical properties of the food and its geometric dimensions. The proposed computational model, the adopted boundary conditions and the calculation method can be one possibility for determining muscular forces. On the other hand, the results obtained can be used to model complex issues concerning the biomechanics of the masticatory system requiring knowledge of loads.

The use of numerical simulation to determine force patterns for the mandibular elevator muscles (the masseter muscle-M, the medial pterygoid muscle-MP and the temporalis muscle-T) during unilateral chewing required, first of all, the construction of a numerical model of the human masticatory system consisting of the mandible, TMJ and muscles. In addition, this necessitated: 1) determining how to support the model and 2) defining and determining the parameters (input data) responsible for the load and initial position (for time t = 0s) of the model.

During the numerical simulations, forces in the muscles were determined for loads corresponding to different foods; therefore the basic parameter responsible for the mandibular force is the food characteristic (i), in the form of the function F i =f(Δh i )—force (F i ) vs. change in high of food specimen (Δh i ), determined in experimental studies. In the numerical model, the above function is decomposed into two functions in which force and change in high of food specimen are time-dependent, i.e., F i (t) and Δh i (t). The general algorithm for determining food characteristics was similar to that presented in (Stróżyk and Bałchanowski 2018). Five foods were prepared for the study (i = c, a, d, b, s): 1) vegetables—carrot (c); 2) fruits—apple (a); 3) sweets—dark chocolate (d) and a chocolate bar (b) and 4) meat and cold cuts—sausage (s), different in terms of structure, mechanical properties and method of production (natural and artificial). The products were purchased from a single selected grocery shop and stored at 7°C until the specimens were prepared. The specimens were then normalised at 21°C for approximately 60 min until the tests began. The aim of the experimental study was not to determine typical mechanical parameters but only the function F i =f(Δh i ). The specimens in terms of dimensions [height (h i ), width (w i ). and length (l i )] were similar to a typical bite of food, while the shape depended on the product (Table 1). To determine the characteristics, 60 specimens were prepared for each food. In the experimental tests, we used our own test stand design (Figure 1) [developed on the basis of a patent application—(Stróżyk, 2021)], consisting of two measuring heads, i.e., upper (maxilla) and lower (mandible), imitating fragments of the dental arches. The lower head is mounted on a cantilever and thus can perform a rotational. In addition, acrylic dental prostheses, used in dental prosthetics, were attached to the heads to make the act of chewing similar to natural chewing. During testing, the test stand was mounted on an Instron 5944 test machine. Based on the information provided in (Fitts et al., 1991; Dick and Wakeling, 2017; Koolstra, 2002), muscle force is dependent on the speed of mandibular movement, while speed is dependent on the mechanical parameters of the food, i.e., high force and low speed (hard food) or low force and high speed (soft food). Furthermore, based on the analysis of (Anderson et al., 2002; Foegeding and Drake, 2007; Foster et al., 2006; Meullenet et al., 2002; Williams et al., 2005) it appeared that the range of reported mandibular movement speeds (in this case, chewing speed) for different foods is extensive. Therefore, it was decided that all foods would be masticated at a constant velocity of v t = 0.02 m/s (Stróżyk and Bałchanowski, 2016; Stróżyk and Bałchanowski, 2018; Stróżyk and Bałchanowski, 2020). Maintaining a constant v t required performing a kinematic analysis of the simulator and, based on this, developing a control program for the Instron 5944 machine (Bluehill ver. 3). In addition, the chewing time (t i ) was determined for each food product based on the nominal food height h i ) and the set velocity v t (Table 1). Figure 1 shows the characteristics of selected products (i), in the form of the function F i = f(Δh i ), developed on the basis of measurement results and elementary statistical calculations.

In order to carry out a kinetostatic analysis of the chewing process of selected foods, a computational model of the human masticatory system was developed, consisting of two elements: a fixed skull (Synbone 8500) and a movable mandible (Synbone 8596). The temporomandibular joints (TMJ) were modelled as shaped joints, with constraints in the form of contact forces between the articular surface of the condyle of the mandible and the articular tuber of the temporal bone (Strozyk and Balchanowski, 2018). The forces in the muscles were modelled using non-linear force vectors applied at the anatomical points of muscle attachment. The numerical model had 4 degrees of freedom, which means that 4 active muscle force unknowns can be determined from the equilibrium equations. Hence, the model assumed the unknowns were the two temporalis muscle forces (FT Wi , FT Ni ). In contrast, the forces in the masseter muscle (FM Wi , FM Ni ) and the medial pterygoid muscles (FMP Wi , FMP Ni ) were replaced by the resultant forces (FV Wi , FV Ni ) (Prium et al., 1980), on the left and right sides (Figure 2). The muscles were modelled using linear kinematic excursions to simulate the elongation and contraction of the muscles that allow the mandible to move in relation to the maxilla during mastication. In the numerical model of the skull−mandible, four linear kinematic excursions were used for the temporalis muscle and the resultant. The active forces FT Wi , FT Ni , in the enforcing qT Wi , qT Ni represent the temporalis muscle, while the active forces FV Wi , FV Ni in the enforcing qV Wi and qV Ni represent the results of the masseter and medial pterygoid muscles (Figure 2). The external force on the model [modelling the occlusal force (Stróżyk et al., 2018; Stróżyk and Bałchanowski, 2016; Stróżyk and Bałchanowski, 2018)] was the force F i , as a function of F i = f(Δh i ), applied at point IW on the occlusal surface of the first molar (46) (Figure 2).

Based on an analysis based on solid mechanics and analytical geometry, it appears that the parameters that will have a significant effect on the alignment of the mandible with respect to the maxilla will be: 1) product height (h i ) and 2) unilateral contraction of the lateral pterygoid muscle. The height of h i is responsible for establishing the distance between the upper and lower incisors, whereas contraction of the muscle on the non-working side (N) is responsible for lateral movement of the mandible, resulting in the appearance of a path of mastication and asymmetrical displacement of the processes in the TMJ during unilateral chewing. Since there is a lack of data in the available literature on the preload (muscle force) and corresponding stiffness of the lateral pterygoid muscle, it was decided to use the effect of its action as a parameter defining mandibular preposition and movement. The appropriate association of the height of the h i and the mastication path (effect) and the appropriate alignment of the mandibular processes at the TMJ enable the initial position of the model to be determined during unilateral chewing. Based on the analysis of data reported in (Bhatka et al., 2004; Buschang et al., 2007; Nishigawa et al., 1997; Piancino et al., 2012; Slavicek, 2010), a hypothetical path of mastication was prepared for each product, in the frontal plane (Figure 3A), in shape close to an ellipse (Hedjazi et al., 2013), in such a way that its geometric parameters (height and width) were synchronised with the food height (h i ), i.e., the initial distance between a pair of corresponding first molars. The path of mastication prepared in this way are the parameter by which it is possible to reproduce the trajectory of the incisal point during a single or one cycle of unilateral chewing. The results of preliminary calculations showed that the introduction of the mastication path significantly affects the values of the length of the condylar trajectories (Δd i ). In addition, analysis of the results showed that the condyles move according to the mandibular kinematics reported in (Weinberg, 1963), i.e., the working condylar rotates and little lateral Bennett’s movement. In contrast, the non-working condylar movement is sloped downward, forward, and medially. For the non-working side, the length of the trajectory of the condylar is: chocolate Δd d = 0.0031 m, chocolate bar Δd b = 0.0068 m, apple Δd a = 0.0079 m, carrot Δd c = 0.0065 m, and sausage Δd s = 0.0102 m, and is on average as much as seven times greater than on the working side. Therefore, in the numerical simulation, it was assumed that there would be a fixed centre of rotation of the mandible in the condylar on the working side, through which the instantaneous axes of rotation would pass. On the non-working side, on the other hand, the condylar will have the possibility of rotation and translation (Figure 3B). The mandibular model is fixed in the TMJ and the origin attachment sites of the masseter, medial pterygoid, and temporalis muscles. Based on the developed chewing paths, there are kinematic contact pairs with 5 degrees of freedom (3 rotations and 2 displacements) in the TMJ on the working and non-working sides. The displacement values depend on the type of food and the chewing side (Table 1). The proposed way of supporting the model, especially in the TMJ, made it possible to simulate the complex movement of the lower incisors during chewing (Figures 2, 3). A detailed description of the model is given in (Stróżyk and Bałchanowski, 2018).

Simulation studies consisted in moving the incisal point from the lower position (point IP Figure 2A) to the endpoint IS (the origin of the coordinate system Figure 2). During the simulation of food chewing, the molar to which the force F i has applied moves with a constant velocity v t = 0.02 m/s. The developed numerical model is universal and can model variable (nonlinear) chewing velocity. In the simulation studies chewing was modelled at a constant speed because the food patterns were determined experimentally on a testing machine for such speed. Numerical simulations were carried out in three stages. In the first, the inverse task of the kinematics of motion of the mandibular-cranial system was solved. This consisted in applying a forced mandibular movement along the longitudinal mastication path and determining the mandibular movement at the TMJ on the W and N sides. During this simulation, the changes in muscle length (contractions) necessary to force the mandibular closure movement were also determined (Figure 4). In the second stage, proper calculations were performed modelling a simple dynamics task simulating mandibular closure resulting from muscle contractions. Forcing qT Wi , qT Ni , qV Wi and qV Ni were used as changes in muscle length (contractions) determined in the simple task. As a result of calculations, the unknown muscle forces were determined: FM Wi , FM Ni , FMP Wi , FMP Ni and FT Wi , FT Ni . Simulations were performed separately for each food. In the third step, based on elementary trigonometric calculations, the principal vector of the temporalis muscle was decomposed (FT Ni ) into three components, respectively for: the anterior temporalis (FAT Ni ), the middle temporalis (FMT Ni ) and the posterior temporalis (FPT Ni ). The distribution was performed based on the assumption that muscle force is proportional to physiological cross-sectional area (PCSA) (Koolstra et al., 1988); Figure 5A shows the schematic position of the temporalis muscle (outline—ABC) in relation to the sagittal plane. On the other hand, Figure 5B shows the position of the temporalis muscle (outline—A'B'C) in the sagittal plane and gives the angles determining the position of the B'C edge (α = 11°) and the line of action of FT Ni (ϕ = 58°) and its components, i.e., FAT Ni (θ = 76°), FMT Ni (γ = 42°) and FPT Ni (β = 22°), in relation to the sagittal axes.

Based on the determined food characteristics (Figure 1), the developed model of the masticatory system (Figure 2), the chewing path (Figure 3), the determined changes in muscle length (Figure 4) and numerical calculations, muscle force patterns were determined, for the masseter muscle, the medial pterygoid muscle and the temporalis muscle and its components (the anterior temporalis, the middle temporalis, the posterior temporalis) (Figure 6). In addition, Tables 2, 3 give the values for maximum bite force and the corresponding values for muscle contraction and total muscle contraction. In contrast, Tables 4, 5 show muscle stiffness and intrinsic strength values. In the general case, the total contraction of the muscle (Δq C ), the contraction corresponding to the maximum force (Δq H ) and the muscle contraction (Δq E ) were determined from Eqs 1–3. During the calculations, the equations were modified in such a way that the above parameters could be determined on selected muscles depending on the food, separately for the working and non-working sides. where: q(0)—initial muscle length for time t = 0s—open mouth, q(t)—muscle end length for time t—closed mouth, q(t H )—muscle length for time t H , corresponding to the position of the mandible in which the muscle force reaches its maximum value (F H ), q(t E )—muscle length determined for time t E from the interval 0÷t,

Preliminary analysis indicates significant differences between patterns, force values and contractions. The variation in muscle behaviour is not surprising, as the diversity of characteristics (Figure 1) generates individual conditions for the masticatory system (Mioche et al., 1999), adapted to functional requirements, depending on the input data (Fitts et al., 1991). For efficient interpretation of the results, muscle patterns were prepared separately for each muscle (Figure 7) and for the components of the temporalis muscle (Figure 5), and separately for the working and non-working sides. Thanks to this, there are 5 muscle patterns on each graph, depending on the food, while the data are given in Tables 2, 3 enable the determination of additional mechanical parameters characterising the muscle functioning.

In the general case, the stiffness (K) is determined from the elementary Eq. 4, in which the force (F) and displacement (x) are directly proportional to each other. Eq. 4 is linear, so its use is limited to problems with a linear relationship between F and x. This means that the stiffness determined for a range of load (0÷F) or displacement (0÷x) has a constant value. The above equation can be used for issues of a non-linear nature, but then the stiffness is only determined for a specific value of F or x. Therefore, the value of stiffness has the characteristics of a comparative parameter based on which it is possible to assess the influence of various factors on the functioning of the muscular system, e.g., 1) different ways of biting, 2) the geometry of the masticatory system resulting from individual characteristics and pathological conditions, 3) foods, 4) the state of dentition and 5) changes resulting from surgical procedures (Pachnicz and Stróżyk, 2021). Based on the considerations above, this paper is limited only to the determination of muscle stiffness (Eq. 4) for the maximum muscle forces (Table 2, 3) for the working (K Wi ) and non-working (K Ni ) sides, respectively. Of course, Eq. 4 was adjusted so that stiffness could be determined for selected muscles, i.e., F→F H and x→ΔqH.

Based on publications (Koolstra et al., 1988; Korioth and Hannam, 1990; May et al., 2001; Pruim et al., 1980) it can be stated that intrinsic strength (k), is used to determine maximum muscle force (F H ) based on Eq. 5 and the relationship between muscle force (F) and EMG measurements according to Eq. 6. where: (EMG F )—scale factor To calculate intrinsic strength, Eq. 5 was used, which was adjusted for the selected muscles depending on the selected foods, for the working (k Wi ) and non-working (k Ni ) sides, respectively—Table 5. To determine intrinsic strength, the data given in Table 1–3 were used and it was assumed that the PCSA values on the working and non-working sides are the same. The values of the temporalis muscle components are not included in Table 5 because their values are identical to those of the principal vector.

In the work presented here, we used our model of unilateral chewing (Figure 2) developed based on: 1) the anatomical structure of the masticatory system (Synbone skull 8500 and mandible 8596), 2) food patterns Figures 1, 3 chewing paths Figure 3A. Based on the assumptions made, dynamic patterns were determined for the mandibular elevator muscles in the form of muscle forces and corresponding contractions. Furthermore, the results obtained allowed the calculation of stiffness of the muscle and intrinsic strength. All values of the above parameters were determined in relation to the food and the working and non-working sides for a single chewing cycle.

The model proposed and the assumptions made are a compromise between a correctly anatomical structure of the masticatory system and the physiological mechanism of mastication (act of chewing) and the possibility of reproducing the conditions (parameters) necessary to reproduce physiologically correct unilateral chewing. Despite the use of simplifications (e.g., same chewing velocity, the muscle is represented by a single principal vector, no friction coefficients, no consideration of the muscle damping coefficient, use of a single geometrical model of the mandible, mechanically correct but basic in terms of the anatomical structure temporomandibular joint), the model meets the basic requirements to carry out calculations in accordance with the principles of solid mechanics and function of the masticatory system. Since the model does not take into account all the parameters responsible for the correct course of the chewing process, it must be assumed that the results of the calculations show how the mandibular elevator muscles may function during the first unilateral chewing cycle.

Similar to the authors’ previous work (Stróżyk and Bałchanowski 2016; Stróżyk and Bałchanowski, 2018; Stróżyk and Bałchanowski, 2020) and other publications (Agrawal et al., 1998; Hiiemae et al., 1996; Koolstra, 2002; Mathevon et al., 1995; Shimada et al., 2012) and the present study, it has been shown that mechanical [F i = f(Δh i )] and geometric [height (h i )] food parameters have a significant influence on the function of the masticatory system, especially the muscular system. Based on the elementary comparative analysis of the food characteristics (Figure 1) with the muscle force patterns (Figures 6, 7), it can be shown that there is a high similarity between them related to the course of the graphs. Such a relationship occurs in the so-called simple mechanisms (Stróżyk and Bałchanowski, 2016; Stróżyk and Bałchanowski, 2018; Stróżyk and Bałchanowski, 2020) characterised by the fact that the response of the system is in the form of displacements and enforcing forces (muscle forces) is correlated with the force (resistance force posed by food). The correlation between the above-mentioned parameters is manifested, among other things, by their values increasing at the same rate and reaching maximum values at the same time. Of course, the above statement is fulfilled for the kinematically correct movement of the masticatory system. Furthermore, it can also be seen that the patterns of the same muscle on the working and non-working sides have similar characteristics but different values (Tables 2, 3). This is due to the asymmetrical alignment of the mandible concerning the sagittal plane, which is a consequence of including the chewing path. In addition, significant differences can be observed, irrespective of the muscle, between the values of maximal muscle forces and maximal contractions, depending on the food. However, the comparison of the values of the contractions corresponding to the maximum forces shows that the differences are not as explicit as for the muscle forces (Tables 2, 3).

The results indicate that not only muscle force patterns but also maximum muscle forces are dependent on food (mechanical parameters). Moreover, analysing their values (Tables 2, 3) it was observed that on the working side, they are higher than on the non-working side. This result proves the correctness of the obtained results and compliance with the elementary static principle, from which it follows that the resultant of external forces (occlusion force) is located closer to the components (muscular forces) with higher values. A detailed analysis taking into account the values of the maximum muscle forces shows that the forces on the non-working side are on average 14% lower than those on the working side, irrespective of the muscle (the masseter, medial pterygoid and temporalis muscles), the division of the temporalis muscle (Figure 5) and the food (Table 1). The difference is not large, but it balances the non-working and the working side, allowing a stable mandible elevation during unilateral chewing for the given boundary conditions. The knowledge of the maximum muscle forces (F HWi and F HNi ) and the maximum occlusion forces (F i , max ) made it possible to show that there is a correlation between them, which can be presented in the form of their ratio (Table 6), corresponding for the working (R FW ) and non-working (R FN ) sides. Based on the data (Figure 1; Tables 2, 3) R FW and R FN can be determined considering: (1) the principal vectors of the masseter, medial pterygoid and temporalis muscles (case I) and/or (2) the principal vectors of the masseter muscle and the medial pterygoid muscle with simultaneous consideration of the division of the temporalis muscle into three parts (case II)—(Figure 5). The results given in Table 6 indicate that the values of R FW and R FN are constant for the muscles and components of the primary vector of the temporalis muscle, of food, but different for the working and non-working sides. The relationship between maximum muscle force and maximum occlusion force can also be represented by a linear relationship (Eqs 7, 8) in which the proportionality coefficient is R FW and R FN . In the literature, it is possible to find publications in which authors model bilateral and unilateral chewing, among others (Choi et al., 2005; Pinheiro and Alves, 2015; Jahadakbar et al., 2016; Lee et al., 2017; Luo et al., 2017; Marková and Gallo, 2016; Orassi et al., 2021; Pachnicz and Stróżyk, 2021; Reina et al., 2007); however, these are models corresponding to static occlusion and the muscle force values used are based, primarily, on the data reported in (Korioth et al., 1992; Korioth and Hannam, 1994) and in (Ferrario et al., 2004). Despite the use of static models in the above-mentioned publications, it was decided to relate the obtained maximum values of muscle forces to those works in which the authors take into account the unilateral loading of the mandible—Table 7. Taking into account the data provided by (Pinheiro and Alves, 2015), it follows that the values of forces in the muscles differ significantly from those reported in Tables 2, 3. According to the authors of the above-mentioned publication, the main role in generating occlusal force is played by the temporalis muscle, which muscular force is greater than that of the masseter muscle and the medial pterygoid muscle, both on the working and non-working sides. On the basis of the results obtained in Tables 2, 3, it follows that during unilateral chewing the muscle forces generated by the masseter muscle and the medial pterygoid muscle are significantly greater than the values determined for the temporalis muscle. A similar regularity regarding the values of muscle forces can be found in the work of (Korioth et al., 1992). Based on the muscle force values reported by (Reina et al., 2007), it appears that the medial pterygoid muscle is the dominant muscle. The analysis performed indicates that on the working side, the medial pterygoid muscle has the highest values, while the superficial the masseter muscle and the temporalis muscle have similar values. On the non-working side, the medial pterygoid muscle and the posterior temporalis muscle have similar values, while the masseter muscle has a value ∼2 times lower than them. In contrast, the results are shown in Tables 2, 3 indicate that regardless of the product and side, the masseter muscle has the highest value, and the temporalis muscle has the lowest value. The muscle force values of the medial pterygoid muscle are only 10% less than those of the masseter muscle. Furthermore, in the works (Korioth et al., 1992; Reina et al., 2007), the temporalis muscle is decomposed into 3 components (Table 7), between which there is the same relationship as shown in Table 3, i.e., the anterior temporalis muscle has the highest value of muscle force, the middle temporalis muscle intermediate and the posterior temporalis muscle the lowest, regardless of the side.

On the basis of the data given in Tables 2, 3, it can be shown that for the working side, the total contraction of the mandibular elevator muscles is 17% lower on average than on the non-working side. Analysis of the results also showed that total contraction was dependent on the initial food height (Table 1), i.e., the distance between a pair of corresponding molars. Detailed analysis showed that the height of the food influences: 1) the size of the chewing path (mainly height, but also width (Figure 3) and 2) muscle contractions (Tables 2, 3). The relationship between food height and chewing path height in the sagittal plane can be described by Tales’ theorem. On the other hand, the relationship between the height of the food (h i ) and the total muscle contraction (Δq Ci ) can be determined: 1) by solving tasks in which the movable part (mandible) is suspended on a tendon (muscle) and is simultaneously supported at a point corresponding to the TMJ or 2) using a parameter that is the ratio (ratio) of their values (Δq Ci /h i ). The ratio was determined separately for the working side (R qW ) and the non-working side (R qN ). The results of the calculations (Table 8) showed that the values of R qW and R qN are: 1) constant for each muscle and 2) dependent on the considered side, i.e., working or non-working side. Considering the temporalis muscle (Figure 5), the values of R qW and R qN (Table 8) can be presented for 2 cases: 1) without considering the division of the temporalis muscle (case I) and 2) with considering the division of the temporalis muscle into three parts (case II) (Figure 5). The relationship between (Δq Ci ) and (h i ) can be represented by the linear relationship (Eqs 9, 10) in which the proportionality coefficient is R qW and R qN . Based on the results presented, it can be assumed that the muscle contraction during unilateral chewing, with the dimensions of the masticatory system fixed, depends on the height of the food bite. Since the height (h i ) changes during biting, therefore, the contraction can be written in the form of a function Δq E = f(Δh i ) - the muscle contraction (Δq E ) vs. the height of food specimen (Δh i ). Moreover, the above parameters are time-dependent, i.e., Δq E (t) and Δh i (t). The above conclusion can also be applied to symmetric incisal biting of foodstuffs. There was no clear relationship between maximum force and corresponding contraction (Tables 2, 3).

As with the other parameters, the stiffness of the muscle (K) is also strongly dependent on: 1) the mechanical properties of the food (texture), 2) the muscle analysed and 3) the side analysed, i.e., working and/or non-working. Preliminary analysis of the results indicates that similar values of stiffness of the muscle (Table 4) can be obtained for products different in texture, characteristics (Figure 1) and method of manufacture. The results given in Table 4 indicate that with respect to K the selected products can be divided into two groups, i.e., 1) chocolate and carrot and 2) bar, apple and sausage. It can also be observed that the K in the first group is, on average 2.7 times higher than that of the second group, irrespective of the type of muscle and the working and non-working side. Furthermore, analysing the K values in relation to the side, it can be seen that the stiffness of the muscles is lower under the non-working side (K Ni ), compared to the working side (K Wi ), by an average of 29%. Analysis of the results given in Tables 2–4 indicates that the greater the stiffness of the muscle, the greater the maximum muscle force. It should be added that the determined patterns of muscle forces (Figures 6, 7) enable the determination of the non-linear stiffness (K i ) according to the general Eq. 11. However, first, with the help of an appropriate tool (mathematical program), an appropriate mathematical equation should be fitted to the data (patterns), i.e., F Ei = f(ΔqEi) or F Ei = F Ei (ΔqEi) (muscle force vs. muscle contraction).where: dF Ei —muscle force gains in relation to food and side of the mandible, d(Δqi)—increases in muscle contraction according to food and side of the mandible The resulting K Wi and K Ni values can be used in numerical calculations where the muscle is modelled as a spring only (Antic et al., 2015; Pachnicz and Stróżyk, 2021). If the calculations involve static loading, it is sufficient to provide a single value corresponding to a given bite force or muscle force, depending on the target. In dynamic calculations, the non-linear function of the muscle must be taken into account (Figure 6), so the stiffness must be entered as a function according to the requirements of the programme. On the other hand, an accurate representation of the performance characteristics of the muscle requires the use of complex models, e.g., 1) Maxwell model, 2) Voigth model, 3) Kelvin model, and 4) Hill model, in which additional parameters regarding elements such as spring elements, damping elements, contractile element and parallel elastic element (elements describing the passive elastic properties of the muscle fibres) must be taken into account (Romero and Alonso, 2016). It is also possible to use the model proposed by Röhrle and Pullan (2007), in which the muscle is a viscoelastic material.

Analysis of the data given in Table 5 has shown that the values of intrinsic strength (k): 1) depend on the food, above all on its mechanical properties, 2) depend on the side analysed, i.e., working and/or non-working, 3) regardless of the side and food considered, are the highest for the temporalis muscle and the lowest for the medial pterygoid muscle. However, for the masseter muscle they are smaller than the medial pterygoid muscle by an average of 28%, and 4) they are identical both for the main vector of the temporalis muscle and its components (this results from an assumption related to the decomposition of the main vector into its components (Figure 5). In the available literature (in recent years), there are few publications concerning the determination of the k-parameter, especially in relation to physiological forces. An interesting publication on the k-parameter, despite the passage of more than 40 years, is the paper by Pruim et al., 1980. In the above-mentioned paper, the authors provide k values determined during experimental studies (in vivo) but during bilateral static occlusion at different levels of occlusion force. Moreover, the authors assume the same k value for all muscles, which makes it difficult to compare the results. However, taking the values obtained for chocolate as a reference (lowest height and maximum occlusion force reached in the mandibular position close to the close-up of the teeth) and considering the temporalis muscle (lowest k), it was decided to compare the results given by Pruim et al., 1980 with the values presented in Table 5. Analysis of the value of the k parameter showed that during unilateral chewing, the temporalis muscle has significantly lower average values (more than 20 times lower than bilateral chewing). Considering the masseter muscle and the medial pterygoid muscle, it turns out that the mean values are also smaller, respectively: 4-fold and 3-fold. Taking into account some of the oldest data concerning the parameter k given in the works by Fick et al. (1904)-1 × 106 N/m2, Morris (1948)-0.9 × 106 N/m2, Ikai and Fukunaga (1968)-0.7 × 106 N/m2 it follows that a value close to the data given in Table 5 is 0.4 × 106 N/m2 determined by Hettinger. (2017). A similar value in the 1980s was given by Weijs and Hillen (1985) i.e., k = 0.37 × 106 N/m2, which is widely used by many authors (Hattori et al., 2003; Koolstra et al., 1988; Koolstra and Eijden, 1992; Korioth and Hannam, 1990; Langenbach and Hannam, 1999; Peck et al., 2000; Zheng et al., 2019). Analysis of the aforementioned papers indicates that the authors use only one value of k = 0.37 × 106 N/m2 or 0.4 × 106 N/m2, although the paper (Slager et al., 1997) used k = 0.35 × 106 N/m2. Furthermore, the authors used one value of k for all muscles. However, according to the data given in Table 5, it follows that the values of k are different and depend on the food product. This indicates that the mechanical parameters of the product also have a significant influence on the values of the k parameter. The same observations but concerning chewing, can be found in (Stróżyk and Bałchanowski, 2018). It should also be added that the equations (Eqs 5, 6) from which the value of k was determined are linear and restricted to a point only, i.e., F H . It is interesting to note that k is very strongly dependent on PCSA, i.e., for the same muscle force but different values of PCSA, different values of k can be obtained, and each will bear the hallmarks (or give the impression) of a correct value. It is also worth noting that there is an analogy between intrinsic strength (k) and the classical definition of strain (σ), i.e., the same formula and the same physical unit in the SI system. Considering the general function σ = f(ε) (tension vs. strain), which is used for the mathematical notation of the characteristics obtained in the classical compression test, it is possible, on its basis, to propose an elementary function for the parameter k, i.e., k = f(ε), where:ε = Δq E /q(0). On the basis of the above function, it can be shown that k is not a constant value during muscle contraction and changes with the change in load similarly to σ. The equation k = f(ε), does not change the physical meaning of intrinsic strength but extends its applications to the full range of changes in muscle force values, not only the maximum values.

The present work shows that the use of a hybrid model based on numerical simulation and experimental studies can provide an advanced (promising) tool for modelling complex processes such as the act of chewing. This paper also shows that muscle forces and mechanical parameters of food can be represented as functions, F i = f(q i ) and F i = f(Δh i ), respectively. In addition, the advantage of the presented model is that it allows rapid modification of geometric parameters and loading functions tailored to individual patient characteristics. The results of such calculations will allow for better fabrication and fitting of implants or prosthetic components to the patient. The determined muscle or food force patterns can also be used in numerical simulations based on deformable body mechanics. Such simulations will make it possible to determine stress, strain and displacement at selected points in the mandible, maxilla, teeth and temporomandibular joint. The conclusions of the present study confirm and complement theories about the mechanisms of the masticatory system during chewing. On the other hand, the computational technique used, using inverse kinematics and dynamics analysis and the application of food patterns, made it possible to understand the influence of mechanical food parameters on muscle force patterns.