Derivation of a finite-element model of lingual deformation during swallowing from the mechanics of mesoscale myofiber tracts obtained by MRI

Abstract

To demonstrate the relationship between lingual myoarchitecture and mechanics during swallowing, we performed a finite-element (FE) simulation of lingual deformation employing mesh aligned with the vector coordinates of myofiber tracts obtained by diffusion tensor imaging with tractography in humans. Material properties of individual elements were depicted in terms of Hill's three-component phenomenological model, assuming that the FE mesh was composed of anisotropic muscle and isotropic connective tissue. Moreover, the mechanical model accounted for elastic constraints by passive and active elements from the superior and inferior directions and the effect of out-of-plane muscles and connective tissue. Passive bolus effects were negligible. Myofiber tract activation was simulated over 500 ms in 1-ms steps following lingual tip association with the hard palate and incorporated specifically the accommodative and propulsive phases of the swallow. Examining the displacement field, active and passive muscle stress, elemental stretch, and strain rate relative to changes of global shape, we demonstrate that lingual reconfiguration during these swallow phases is characterized by (in sequence) the following: 1) lingual tip elevation and shortening in the anterior-posterior direction; 2) inferior displacement related to hyoglossus contraction at its inferior-most position; and 3) dominant clockwise rotation related to regional contraction of the genioglossus and contraction of the hyoglossus following anterior displacement. These simulations demonstrate that lingual deformation during the indicated phases of swallowing requires temporally patterned activation of intrinsic and extrinsic muscles and delineate a method to ascertain the mechanics of normal and pathological swallowing.

Fig. 1.

Development of finite element (FE) mesh from images of human lingual myofiber tracts obtained by diffusion tensor MRI with tractography. A: anatomy of the muscles present in the midsagittal plane of the tongue, shown in schematic view displaying several principal muscle groups and their bony attachments. Distinguished muscle groups include genioglossus (gg), verticalis (v), geniohyoid (gh), superior longitudinalis (sl), and inferior longitudinalis (il), with connections shown to the symphysis and the hyoid bones. B: diffusion tensor imaging (DTI) tractography images of myofiber tracts noted in A. For these imaging experiments, diffusion-weighted gradients were applied in 90 directions, employing single-shot echo-planar spatial encoding with repetition time = 3,000 ms, echo time = 80 ms, field of view 192 mm × 192 mm, slice thickness 3 mm, and b-value of 500 s/mm², followed by the streamline construction of multivoxel myofiber tracts along the maximum diffusion vector per voxel. C: DTI tractography myofiber tracts displayed along with the points of insertion into midsagittal and non-midsagittal structures: geniohyoid (gh), hyoglossus (hg), styloglossus (sg), and inferior longitudinalis (il). D: FE mesh whose elemental alignment is derived from the principal diffusion direction per voxel obtained through DTI, including midsagittal and out-of-plane muscles, as well as boundary conditions such as attachments to the symphysis and hyoid bones. Note specifically that the hyoid bone is attached to fixed structures via elastic (connective tissue) structures.

Fig. 2.

Multiscale model of tongue contraction, displayed in the manner in which lingual mechanical function is characterized across relevant spatial scales. A: DTI tractography image of the midsagittal plane of the tongue discretized into FEs. Actual tongue dimensions are denoted along horizontal and vertical axis in millimeters. B: diagram depicting the muscle fibers contained within a characteristic three-dimensional (3D) FE, including denoted integration points and the principal direction muscle fibers, ξ. C: elongation of an individual muscle fiber at the indicated spatial scale under stress is denoted by σ_ξξ. L, length; ΔL, change in length; ^tL, length at current time t; L_o, slack (relaxed) length. D: Hill's three-element model, displayed as the contractile element (CE) connected to series elastic elements (SEE) and parallel elastic elements (PEE).

Fig. 3.

Hill's functional model of skeletal muscle. Muscle mechanical function will be represented in terms of the Hill's phenomenological model. By this model, active muscle is composed of three elements: CE, SEE, and PEE. CE represents contractile (sarcomeric) part of muscle, which generates active force and increases overall muscle stiffness. In the relaxed state, this element generates zero tension and zero stiffness. The SEE includes the elasticity of actin and myosin filaments and tendons, and the PEE includes the elasticity of the surrounding connective tissue and the elasticity of the noncontractile part of the muscle cell's cytoskeleton. The lengths of CE, SEE, and PEE in the relaxed state are denoted as L_mo, L_so, and L_po, respectively. The relaxed muscle length is equal to L_po = L_mo + L_so. The lengths of the Hill's model components during contraction are ^tL_m = L_mo + ^tU_m, ^tL_s = L_so + ^tU_s and ^tL_p = L_po + ^tU_p, where ^tU_m, ^tU_s, and ^tU_p are displacements, i.e., change of length from undeformed configuration of CE, SEE, and PEE, respectively.

Fig. 4.

The FEs mesh and assigned mechanical characteristics in the Gaussian points for each FE. A: the mesh overlaid on the original DTI image slice. The thick yellow lines indicate muscle fiber groups for which mechanical characteristic are assigned. The specific muscle groups and subgroups employed in the current simulations are identified. B: Gordon's tension-length relationship indicates the fraction of the active muscle tension at arbitrary length, ^tL, in relation to the tension at a slack (relaxed) length, L_o. Muscle length can be prescribed as a sarcomere length at a micrometer scale that, after normalization by the slack length, can be represented as stretch, λ_mξ, at the macroscopic scale. C: activation function denotes fraction of the isometric force at time t to the maximum active force that can be generated at the same length, ^tL. The tabular values of the activation functions associated with the FE groups denoted in A are shown in Table 1. f_Lm, maximum isometric force and muscle length; f_a, activation function.

Fig. 5.

Simulated displacement field relative to tissue deformation during swallowing. For Figs. 5–9, the deformed mesh (black lines) is superimposed on the original mesh (gray lines in background) for perspective to the initial configuration of the tongue. Elemental displacement is defined as a vector of the difference between current position of a material position (at time t) and the position of the same point in initial configuration (at t = 0). The color code represents magnitude of the displacement vector in millimeters. Shown is prototypical lingual deformation during swallowing, initiating with lingual tip contact with the hard palate and extending through the accommodative and propulsive phases. These results demonstrate the significance of the upward deformation of the lingual tip and anterior-posterior shortening, followed by inferior displacement (principally related to hyoglossus contraction) and clockwise rotation (principally related to superior-inferior sequential contractions of the genioglossus regions).

Fig. 6.

Simulated active muscle stress during swallowing. Active stress is defined in Eq. 8 and denoted as σ_m. The color code represents the magnitude of active muscle tension (kPa) in the direction of muscle fibers ξ. The dominant patters of active stress (largest magnitude) is associated with lingual rotation about the genioid process and is associated with sequential regional contraction of the regions of the genioglossus.

Fig. 7.

Simulated passive stress field during swallowing. Passive stress is defined by Eq. 12 and represents the stress in the parallel elastic element, σ_ξξ^E, in direction of muscle fiber ξ. The color code represents the magnitude of σ_ξξ^E in kilopascals. σ^E is the passive stress tensor, and σ_ξξ^E is the most relevant component of this tensor. These results demonstrate that, within the current two-dimensional framework, passive stress tends to align with active stress in the apparent distribution of the superior regions of the genioglossus. The passive stress in both the longitudinal and the transverse directions considers only the passive elasticity and does not account for the increase of stiffness during contraction caused by bound cross bridges.

Fig. 8.

Simulated local stretch during swallowing. Local stretch is defined as λ_p = ^tL_p/L_op = (1 + ^tU_p)/L_op. Because, as the length of parallel elastic component is equal to the muscle length (Fig. 3), i.e., ^tL_p = ^tL, the local parallel elastic stretch λ_p is equal to macroscopic stretch in Gaussian point of a FE, λ (Eq. A1). The color code represents the magnitude of the muscle stretch in direction of muscle fibers ξ. Stretch in the direction of indicated muscle fibers may be displayed as regions of compression (gray) or extension (red). The pattern of stretch displayed is similar to that exhibited by passive stress (Fig. 7), although relative magnitudes are skewed due to material and geometric nonlinearities.

Fig. 9.

Simulated local strain rate during swallowing. Local strain rate is computed to serve as a comparison with prior experimental work deriving aligned strain rate from lingual pressure-gated phase-contrast MRI and DTI tractography (4, 5). In the current simulation, local strain rate is calculated as a difference between consecutive images of the local stretches divided by time intervals between times when these images were obtained. The images of local stretches were collected at 1-ms intervals, and then the strain rates were collected at specified times denoted in the left corners of the presented images. The color code represents the magnitude of the strain rate in direction of muscle fibers ξ. Strain rate is displayed as regions of contraction rates (red) or elongation rates (gray).