In this paper, we describe an open-source physical simulation system that combines both FEM and multibody capabilities, with an emphasis on computational efficiency. We then show the application of this system to a novel integrated model of the tongue, jaw, and hyoid complex that requires FEM/multibody capability to physically model the dynamic interactions when soft tissue and hard bodies are connected and interact. The physical simulation system is described in Section 3, along with specific discussions of attachments between bodies, contact handling, solution techniques, and a comparison of solution results with ANSYS. Our demonstration of the simulation system with a model of coupled jaw-tongue-hyoid dynamics is described in Section 4, along with a number of simulations of combined tongue-jaw actions to illustrate the types of complexities that arise in this region and how our system manages them.

) simulations demonstrating effects of coupled jaw-tongue-hyoid dynamics, and (7) preliminary comparisons of simulations results with published human data (Section 4.3).

control strategies are not completely understood, especially for complex motor actions, such as speech production, that involve the coordination of multiple structures in very fast movements up to 80 cm/s [START_REF] Perkell | Physiology of speech production: results and implications of a quantitative cineradiographic study[END_REF]. Comprehensive dynamic models of this region will enhance our understanding of both its normal physiological function and its dysfunctions, such as obstructive sleep apnea [START_REF] Huang | Use of computational modeling to predict responses to upper airway surgery in obstructive sleep apnea[END_REF], swallowing disorders [START_REF] Lundy | Aspiration: cause and implications[END_REF], and speech pathologies [START_REF] Buchaillard | Simulations of the consequences of tongue surgery on tongue mobility: implications for speech production in post-surgery conditions[END_REF]. These models will also help to evaluate the consequences of maxillofacial surgical interventions [START_REF] Chabanas | Patient specific finite element model of the face soft tissues for computer-assisted maxillofacial surgery[END_REF].

Modeling complex structures for biomedical applications is enhanced by an interactive environment capable of simulating the dynamics of coupled rigid and deformable structures, together with contact and constraints, in a reasonable amount of time (seconds as opposed to minutes or hours). Traditionally, deformable structures are modeled using finite element method (FEM) approaches [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF][START_REF] Bonet | Nonlinear Continuum Mechanics for Finite Element Analysis[END_REF] that emphasize boundary conditions and nodal loading. Rigid-body structures, on the other hand, are modeled using multibody techniques [START_REF] Shabana | Dynamics of Multibody Systems[END_REF] that emphasize constraints and contact. Existing software packages tend to be oriented to either the FEM (e.g. ANSYS, SIMULIA) or multibody (e.g. SolidWorks, ADAMS) approach and creating hybrid models within either is generally not easy. Moreover, solution times can be very slow: a 1-s simulation of the FEM tongue model described in [START_REF] Buchaillard | A biomechanical model of cardinal vowel production: muscle activations and the impact of gravity on tongue positioning[END_REF] can require many hours of computing time.

RELATED WORK

Biomechanical models of individual structures in human upper airway have been developed and used since the 1970s. Model complexity has increased due to both the acquisition of new knowledge about anatomical, neurophysiological, and physical characteristics of the articulators, and the vast growth in computational capacity. Early models were based on a 2D mid-sagittal representation of the airway [START_REF] Perkell | A physiologically-oriented model of tongue activity in speech production[END_REF][START_REF] Laboissiere | The control of multi-muscle systems: human jaw and hyoid movements[END_REF][START_REF] Payan | Synthesis of VV sequences with a 2D biomechanical tongue model controlled by the Equilibrium Point Hypothesis[END_REF][START_REF] Sanguineti | A dynamic biomechanical model for neural control of speech production[END_REF][START_REF] Payan | Biomechanical models to simulate consequences of maxillofacial surgery[END_REF][START_REF] Perrier | Influences of tongue biomechanics on speech movements during the production of velar stop consonants: a modeling study[END_REF], whereas recent models attend to the 3D structure of the articulators. 3D biomechanical models have been proposed for the face [START_REF] Chabanas | Patient specific finite element model of the face soft tissues for computer-assisted maxillofacial surgery[END_REF][START_REF] Sifakis | Automatic determination of facial muscle activations from sparse motion capture marker data[END_REF][START_REF] Nazari | Simulation of dynamic orofacial movements using a constitutive law varying with muscle activation[END_REF], jaws [START_REF] Hannam | A dynamic model of jaw and hyoid biomechanics during chewing[END_REF][START_REF] Tuijta | Differences in loading of the temporomandibular joint during opening and closing of the jaw[END_REF], tongue [START_REF] Buchaillard | A biomechanical model of cardinal vowel production: muscle activations and the impact of gravity on tongue positioning[END_REF][START_REF] Dang | Construction and control of a physiological articulatory model[END_REF][START_REF] Gerard | 3D Biomechanical Tongue Modelling to Study Speech Production[END_REF], and larynx [START_REF] Hunter | A three-dimensional model of vocal fold abduction/adduction[END_REF]. Most authors have used the FEM to discretize the partial differential equations governing the continuum mechanics of soft tissues. However, because of the significant technical challenges with properly modeling the dynamical coupling between soft bodies (tongue, lips, soft palate) and hard structures (jaw, hyoid bone, hard palate), most 3D models have proposed to focus on one single articulator. Hence, to our knowledge, our approach is unique for providing a 3D modeling framework simulating full dynamical interaction between soft and rigid articulators. Such dynamic coupling is important, especially in speech production where it has been clearly shown that consideration of the dynamic interaction of vocal tract bony and soft structures is needed to correctly account for orofacial dynamics (see [START_REF] Sanguineti | A dynamic biomechanical model for neural control of speech production[END_REF] with the use of a quite simple 2D FEM model).

Deformable bodies often form stiff systems that require implicit integrators for computationally efficient solutions, and hence rigid-deformable coupling requires the ability to apply implicit integration techniques to constrained multibody problems. Theoretical work on this problem has been reported in [START_REF] Anitescu | A time-stepping method for stiff multibody dynamics with contact and friction[END_REF][START_REF] Potra | A linearly implicit trapezoidal method for integrating stiff multibody dynamics with contact, joints, and friction[END_REF][START_REF] Lunk | Solving constrained mechanical systems by the family of Newmark and -methods[END_REF]. Rigid-deformable body coupling has also been studied by the computer graphics community [START_REF] Shinar | Two-way coupling of rigid and deformable bodies[END_REF]. For surgical simulation, fast rigid-deformable contact is required for instrument interactions; methods to achieve this include penalty forces [START_REF] Meseure | A physically-based virtual environment dedicated to surgical simulation[END_REF] and a linear complementarity approach optimized for linear materials [START_REF] Saupin | Contact Model for Haptic Medical Simulations[END_REF].

Several open simulation systems and architectures have been presented to the biomedical community in recent years. OpenSim [START_REF] Delp | OpenSim: opensource software to create and analyze dynamic simulations of movement[END_REF] is a multibody simulator with inverse modeling capabilities designed for musculoskeletal work, whereas FEBio [START_REF] Weiss | FEBio-a finite element tool for biomechanics[END_REF] is a finite element toolkit with special support for tissue modeling and some support for rigid-bodies, contact and constraints. Systems geared toward surgical training include Gipsi [START_REF] Cavusoglu | Gipsi: a framework for open source/open architecture software development for organ-level surgical simulation[END_REF] and Spring [START_REF] Montgomery | Spring: a general framework for collaborative, real-time surgical simulation[END_REF]. Sofa [START_REF] Allard | Sofa an open source framework for medical simulation. In Medicine Meets Virtual Reality 15-In Vivo[END_REF] provides a general software architecture in which models can be partitioned into different submodels for simulating appearance, behavior, and/or haptic response. To the best of our knowledge, none of these open simulation systems yet provide an interactive environment with fully coupled FEM/multibody capabilities.

COUPLED HARD-SOFT TISSUE SIMULATION WITH CONTACT AND CONSTRAINTS

The ArtiSynth system

Our physical simulation system is embedded within ArtiSynth (www.artisynth.org), an opensource Java-based biomechanical modeling toolkit developed at the University of British Columbia. Originally designed for speech applications [START_REF] Fels | Artisynth: towards realizing an extensible, portable 3d articulatory speech synthesizer[END_REF], ArtiSynth has evolved into a tool for physiological research (particularly neuromotor control) and clinical treatment planning. Artisynth models are created in Java code, using the packages and classes of the ArtiSynth API. Graphical editing and model creation is also supported. Applications to date have focused on the jaw and oral region [START_REF] Hannam | A dynamic model of jaw and hyoid biomechanics during chewing[END_REF][START_REF] Stavness | Predicting muscle patterns for hemimandibulectomy models[END_REF], but it is broadly applicable to biomechanical modeling in general. Key system features include [START_REF] Perkell | Physiology of speech production: results and implications of a quantitative cineradiographic study[END_REF] an architecture that supports extensive interactivity, including graphically based model editing and simulation control (Figure 1(b)), and (2) a physics engine that combines FEM and multibody capabilities, with constraints and contact handling, as described below.

Physical simulation framework

Artisynth models consist of a hierarchy of components, which include dynamic components such as particles, FEM nodes, or rigid-bodies, force effectors such as point-to-point muscles (including Hill and other types), linear or nonlinear finite elements, and constraints such as joints or collision specifiers. FEM capabilities include support for tetrahedral, hexahedral, and some higher-order elements, along with both linear and large-deformation behaviors, including corotated linear [START_REF] Müller | Interactive virtual materials[END_REF] and hyperelastic materials. We now describe the mathematical framework for the dynamic simulation of these models.

Let q and u be the generalized positions and velocities of all the dynamical components, with q related to u by q = Qu (Q generally equals the identity, except for components such as rigid-bodies, where it maps angular velocity onto the derivative of a unit quaternion). Let f(q, u, t) be the force produced by all the force effector components (including the finite elements), and let M be the (block-diagonal) composite mass matrix. By representing rigid-body velocity and acceleration in body coordinates, we can ensure that M is constant. Newton's second law then gives M u = f(q, u, t).

(1)

In addition, bilateral and unilateral constraints give rise to locally linear constraints on u of the form

G(q)u = 0, N(q)u 0. ( 2
)
Bilateral constraints include rigid-body joints, FEM incompressibility associated with the mixed u-P formulation [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF], and point-surface constraints, whereas unilateral constraints include contact and joint limits. Constraints give rise to constraint forces (in the directions G(q) T and N(q) T ) that supplement the forces of (1) in order to enforce the constraint conditions. In addition, for unilateral constraints, we have a complementarity condition in which Nu>0 implies no constraint force, and a constraint force implies Nu = 0. Any given constraint usually involves only a few dynamic components and hence G and N are generally sparse. Solving the equations of motion requires integrating (1) together with [START_REF] Huang | Use of computational modeling to predict responses to upper airway surgery in obstructive sleep apnea[END_REF]. The presence of deformable bodies generally makes this system stiff, implying the need for an implicit integrator to obtain efficient performance ‡ . For the work described in this paper, we use a semi-implicit second-order Newmark integrator [START_REF] Lunk | Solving constrained mechanical systems by the family of Newmark and -methods[END_REF], with = 1 2 and = 1 4 (also known as the trapezoidal rule). Letting k index values at a particular time step, and h denote the time step size, this leads to the update rules

u k+1 = u k + h 2 ( uk + uk+1 ), q k+1 = q k + h 2 (Q k u k +Q k+1 u k+1 ), (3)
subject to

G k+1 u k+1 = 0, N k+1 u k+1 0. (4)
Since G and N tend to vary slowly between time steps, we can approximate (4) using

G k u k+1 = g k , N k u k+1 n k , (5)
where

g k ≡ -h Ġk u k and n k ≡ -h Ṅk u k . Likewise, we use the approximation Q k+1 ≈ Q k +h Qk .
For uk+1 , recalling that M is constant, an estimate of the (unconstrained) value of uk+1 can be obtained from uk+1 ≈ M -1 f k+1 , with f k+1 approximated by the first-order Taylor series Placing this into the expression for u k+1 in (3), multiplying by M, noting that q = h/2(Q k u k + Q k+1 u k+1 ) and u = u k+1 -u k , and incorporating the constraints [START_REF] Chabanas | Patient specific finite element model of the face soft tissues for computer-assisted maxillofacial surgery[END_REF], we obtain the system

    Mk -G kT -N kT G k 0 0 N k 0 0         u k+1 k z     +     -Mu k -h fk -g k -n k     =    0 0 w    , 0 z ⊥ w 0. ( 6
)
where w is a slack variable, k and z give the average constraint impulses over the time step, and

Mk ≡ M- h 2 *f k *u - h 2 4 *f k *q Q k+1 and fk ≡ f k - 1 2 *f k *u u k + h 4 *f k *q Q k u k .
The complementarity condition for unilateral constraints is enforced by 0 z ⊥ w 0. A more detailed explanation of this formulation can be found in [START_REF] Potra | A linearly implicit trapezoidal method for integrating stiff multibody dynamics with contact, joints, and friction[END_REF]. System ( 6) is a mixed linear complementarity problem (LCP), a single solve of which is required to determine u k+1 for each semi-implicit integration step. Other types of integrators give rise to similar systems. A fully implicit integrator (not currently implemented in ArtiSynth) would require (6) to be applied iteratively at each time step.

For finite element models, the localized stiffness and damping matrices are embedded within *f k /*q and *f k /*u, which means that for models dominated by FEM components M will have an FEM sparsity structure. ArtiSynth FEMs also use a lumped mass model, which ensures that M is block-diagonal and makes it easier to interconnect FEMs with mass-spring and rigid-body components.

Friction, damping, and stabilization

Coulomb (dry) friction can be included by extending [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF] to include either a linearized friction cone [START_REF] Anitescu | A time-stepping method for stiff multibody dynamics with contact and friction[END_REF][START_REF] Potra | A linearly implicit trapezoidal method for integrating stiff multibody dynamics with contact, joints, and friction[END_REF] or a (more approximate but easier to solve) box friction [START_REF] Lacoursière | Ghosts and machines: regularized variational methods for interactive simulations of multibodies with dry frictional contacts[END_REF]. ArtiSynth currently implements box friction, and since the friction in our system tends to be quite small, we apply this as a post hoc correction to u k+1 (in the manner of [START_REF] Shinar | Two-way coupling of rigid and deformable bodies[END_REF]), using a simplified version of ( 6), with M instead of M and extra constraints added in the tangential directions at contact points.

Different forms of viscous damping are available, including translational and rotary damping applied directly to particles and rigid-bodies, and damping terms embedded in point-to-point springs and muscle actuators. For FEM models, Rayleigh damping is available, which takes the form

D F = M F + K F ,
where M F is the portion of the (lumped) mass matrix associated with the FEM nodes and K F is the (instantaneous) FEM stiffness matrix. D F is then embedded within the overall system matrix *f/*u.

In addition to solving for velocities, it is also necessary to correct positions to account for drift from the constraints, including interpenetrations arising from contact. This can be done at each time step using a modified form of (6), which computes an impulse q that corrects the positions while honoring the constraints:

    Mk -G kT -N kT G k 0 0 N k 0 0        q k z    +    0 d g d n    =    0 0 w    , 0 z ⊥ w 0, (7)
where d g and d n are the constraint displacements that must be corrected. If the corrections are sufficiently small, it is often permissible to use M in place of Mk , which improves solution efficiency since M is constant and block-diagonal. While such stabilization can sometimes be incorporated directly into (6) [START_REF] Anitescu | A constraint-stabilized time-stepping approach for rigid multibody dynamics with joints, contact and friction[END_REF], we prefer to perform the position correction separately as this (a) allows for the possibility of an iterative Figure 2. Problems with decoupling constraints from the velocity solve: in (a), a uniform 3×6×3 FEM grid of linear material with a Poisson ratio of 0 is about to be compressed by a block. The decoupled solve causes the top contacting nodes to bunch up on the surface (b), completely squashing the top two element layers, while the lower nodes hardly move at all; the coupled solve (c) causes the correct uniform displacement for all nodes. In (d), a decoupled solve causes a tongue model attached to a jaw to exhibit large vertical errors when the jaw clenches upwards against the bite plane.

correction in the case of larger errors, and (b) explicitly separates the computed velocities from the impulses used to correct errors.

System solution and complexity

For notational convenience, in this section we will drop the k superscripts from M, G, N, g, n, and f in [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF] and assume that these quantities are all evaluated at time step k. System ( 6) is a large, sparse mixed LCP [START_REF] Richard | The Linear Complementarity Problem[END_REF] that is not particularly easy to solve, given the unilateral constraints and the fact that M is not block-diagonal. If M is symmetric positive definite (SPD), it is equivalent to a convex quadratic program. If there are no unilateral constraints (N = ∅), then it reduces to a linear Karush-Kuhn-Tucker (KKT) system.

Generally, M is symmetric (unsymmetric terms sometimes arise from rotational effects, but these are usually small enough to ignore) and hence will also be SPD for small enough h (since M is SPD). However, the resulting system is still harder to solve than non-stiff multibody systems where M = M. This is because M, while still sparse, is not block-diagonal. Multi-body systems are often solved using the projected Gauss-Seidel method [START_REF] Lacoursière | Ghosts and machines: regularized variational methods for interactive simulations of multibodies with dry frictional contacts[END_REF]. However, this involves a sequence of iterations, each requiring the computation of

G i M-1 G T i or N i M-1 N T i ,
which is easy to do for a block-diagonal M but much more costly for M.

It is tempting to follow the approach we use for friction and decouple the velocity and constraint solves, by first computing u * = M-1 (Mu k +h f) and then applying constraints to u * in a post hoc fashion, using a version of ( 6) in which M is replaced with M. This can be done by various methods, including the Gauss-Seidel iteration, and is equivalent to projecting u * onto the space of legal velocities. Unfortunately, this does not propagate constraint effects properly throughout the system, and can result in very large errors when the constraint forces are large, as illustrated in Figure 2.

At present, ArtiSynth solves (6) by using a Schur complement to turn it into a dense regular LCP

NA -1 NT z+ NA -1 b-n = w, 0 z ⊥ w 0, (8)
where which is solved using Keller's algorithm [START_REF] Lacoursière | Ghosts and machines: regularized variational methods for interactive simulations of multibodies with dry frictional contacts[END_REF]. u k+1 and k can then be obtained using backsubstitution:

A ≡ M -G T G 0 , N ≡ (N 0), b ≡ Mu k +h f g Copyright ᭧
u k+1 k = A -1 (b+ NT z). (9)
Keller's algorithm is a pivoting method with an expected complexity of O(m 3 ), where m is the number of unilateral constraints. In addition, forming [START_REF] Shabana | Dynamics of Multibody Systems[END_REF] and the back solve of (9) requires m +1 solves of a system involving A. This is done using the Pardiso sparse direct solver [START_REF] Schenk | Solving unsymmetric sparse systems of linear equations with PARDISO[END_REF], and entails a once-per-step factoring of A, plus m +1 solve operations. Experimentally, we have determined that the complexity of factoring A (using Pardiso) for 3D FEM type problems is roughly O(n 1.7 ), where n is the size of A (Figure 3). Similarly, we have also determined that the complexity of solving a factored A is roughly O(n 1.3 ). Hence, we can expect the overall complexity for solving [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF] to be O(m 3 )+m O(n 1.3 )+ O(n 1.7 ).

This works well provided that the number of unilateral constraints m is small. To help achieve this, we can sometimes treat the unilateral constraints arising from contact as bilateral constraints (i.e. entries in G) on a per-step basis, as described further in Section 3.6. Specific solution times for our jaw-tongue-hyoid model, along with a discussion of possibly more efficient solution methods, are presented in Section 5.

Attachments between bodies

In creating comprehensive anatomical models, it is often necessary to attach various bodies together. Most typically, this is done by connecting points of one body to specific locations on another body. For example, FEM nodes may be attached to particular spots on a rigid body, or to other nodes of a different FEM model.

To facilitate this, ArtiSynth provides the ability to attach a dynamic component to one or more master components. Let the set of attached components be denoted by , and the remaining set of unattached active components be denoted by . In general, the velocity u j of an attached component is related to the velocities u of the active components by a locally linear velocity constraint of the form u j +G j u = 0. G j will be sparse except for entries corresponding to the master components to which j is attached. Letting G denote the composite matrix formed from G j for all attached components, we have

Iu +G u = 0
for the constraints that enforce all attachments.

We could simply add these constraints to [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF] and solve the resulting system, but this would increase both the system size and solution time. Instead, we use the attachments to actually reduce the size of [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF]. Consider first the subsystem involving only bilateral constraints. As in Section 3.4, we drop the k superscripts from M, G, N, g, n, and f in ( 6) and assume that these quantities are all evaluated at time step k. Letting b ≡ Mu k +h f and partitioning the system into active and attached components yields

       M M G T G T M M G T I G G 0 0 G I 0 0               u k+1 u k+1 k k        =       b b g 0      
.

The identity submatrices make it easy to solve for u k+1 and k :

u k+1 = -G u k+1 , k = b -M u k+1 + M G u k+1 -G T k
and hence reduce the system to

M G T G 0 u k+1 k = b g , (10)
where

M ≡ P MP T , G ≡ GP T , b ≡ Pb with P ≡ (I -G T ).
Similarly, unilateral constraints can be reduced via N = NP T . The reduction operation can be performed in O(n) time and results in a system that is less sparse but generally faster to solve than the original.

Most attachments in ArtiSynth are point-based, with the most common kind being the attachment of an FEM node to a rigid body. It is also possible to attach FEM nodes to the faces and edges of an FEM element, allowing us to handle the so-called 'T-junction' problem and create FEM models with non-conforming element faces. This is quite useful for creating localized subdivisions of particular elements, particularly hexahedrons.

Contact handling

Collision detection can be enabled between any combination of rigid or deformable bodies. It is assumed that the bodies in question contain a triangular surface mesh. A bounding-box hierarchy is used to determine if any two surfaces' meshes intersect. If they do, then a tracing algorithm (similar to [START_REF] Aftosmis | Robust and efficient cartesian mesh generation for component-based geometry[END_REF]) is used to compute all the intersection contours between the two meshes as shown in Figure 4. Such contour tracing can be done relatively quickly but does require the use of robust geometry predicates similar to those in [START_REF] Edelsbrunner | Simulation of simplicity: a technique to cope with degenerate cases in geometric algorithms[END_REF]; this is particularly true because collision conditions tend to drive the contacting surfaces into degenerate mesh configurations.

Determining the intersection contour allows us to easily create a set of constraints for correcting the interpenetration and preventing interpenetrating velocities. It also provides a good estimate of the contact area, which can be used for determining contact pressure.

For collisions involving a deformable body, we locate all mesh vertices that are interior to the contour. Each such vertex corresponds to surface FEM node, which is interpenetrating the other body. For each interpenetrating node, we then find the nearest point and face on the opposite mesh, and use the face's normal n as a contact normal. A linear one-dimensional constraint is then created that prohibits relative motion in the negative normal direction between the node and nearest point. If the opposite face is located on a deformable body, this results in a constraint between the velocity of the node v n and the velocities v 0 , v 1 , v 2 of three nodes associated with the nearest face:

n T v n -w 0 n T v 0 -w 1 n T v 1 -w 2 n T v 2 0,
where w 0 , w 1 , and w 2 are the barycentric coordinates of the nearest point with respect to the face.

If the opposite face is located on a rigid body, then the constraint is between v n and the body's translational and angular velocities, v b and x b (expressed in body coordinates):

n T v n -n T Rv b +n T R(p×x b ) 0,
where R is the rotational transformation from body to world coordinates and p is the location of the nearest point in body coordinates. Each of these constraints can be expressed in the general form N i u 0, where N i is a row of N and u is the vector of all velocities. These constraints serve both to prevent interpenetrating velocities in the contact direction, and to remove interpenetrations during the position correction step [START_REF] Bonet | Nonlinear Continuum Mechanics for Finite Element Analysis[END_REF], with the correction distance taken to be the distance d between the interpenetrating node and its nearest point on the other body.

Intersection contours are also used to determine contact constraints for rigid-body/rigid-body contact, although we omit the details here for brevity.

As mentioned in Section 3.4, the solution time of (6) can be greatly improved if some contact constraints can be temporarily treated as bilateral constraints within a particular time step. By default, ArtiSynth does this for contact involving deformable bodies, since such bodies have many degrees-of-freedom and their contact constraints tend to be somewhat decoupled. To prevent sticking, each contact's vertex-face pair is stored between time steps, and if it reappears in the next step, it is used as a contact constraint only if its corresponding k value computed in ( 6) is greater than or equal to 0, implying that there is no force trying to make it separate. This is in effect an active set method, with the active set used to solve [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF] being updated between steps.

It should be noted that we do not claim that the collision handling scheme described here is optimal for all applications. In particular, we do not currently implement edge-edge type contacts, and so there can be some interpenetration that depends on the coarseness of the surface meshes. However, the collision handling is properly isolated from the rest of the simulation, and other collision handling schemes can be easily used as long as they provide a set of constraints for enforcing the contact and resolving interpenetrations.

Simulation engine summary

The complete ArtiSynth simulation engine is summarized below. It uses the concept of [START_REF] Shinar | Two-way coupling of rigid and deformable bodies[END_REF][START_REF] Guendelman | Nonconvex rigid bodies with stacking[END_REF] whereby velocities are computed in advance of positions, subject to constraints, to help prevent constraint violations during the subsequent position computation.

1. Compute contacts (as per Section 3.6) and the bilateral and unilateral constraint matrices G k and N k . 2. Correct positions q k to remove interpenetration and drift errors, using [START_REF] Bonet | Nonlinear Continuum Mechanics for Finite Element Analysis[END_REF]. 3. If necessary, adjust G k and N k to reflect changes in q. 4. Solve for u k+1 using (6). 5. Adjust velocities u k+1 for dry friction, as described in Section 3.3. 6. Compute new positions:

q k+1 = q k +h/2(Q k+1 u k+1 +Q k u k ).
This algorithm is generally applicable to any rigid-deformable body dynamics. In Section 4, we demonstrate that it can be effectively used for simulating the combined dynamics of the jaw, tongue, and hyoid bone.

In the absence of constraints, the above system turns into a trapezoidal rule solution of a regular ODE, for which global errors can be expected to be proportional to O(h 2 ). The inclusion of constraints, particularly non-smooth unilateral ones, makes formal convergence and error analysis more difficult. However, the main velocity update ( 6) is the same as that described in [START_REF] Potra | A linearly implicit trapezoidal method for integrating stiff multibody dynamics with contact, joints, and friction[END_REF], which is shown to be stable and has second-order convergence under certain assumptions.

Generally speaking, our method is a time-stepping scheme that uses fixed (or adapatively varying) time steps, as opposed to an event-driven scheme in which the integration time intervals are precisely aligned with contact events but which becomes impractical in the presence of large numbers of contacts. To the extent to which results exist, time-stepping schemes are typically shown to have less accuracy but better convergence properties than event-driven ones [START_REF] Brogliato | Numerical Methods for Nonsmooth Dynamical Systems[END_REF].

Validation using ANSYS

To help assess the performance of our integration scheme, we compared it against the commercial finite element package ANSYS for two test examples: a beam, fixed at one end and allowed to fall under gravity, and a cube, resting on a flat surface and subjected to a downward load applied to several top nodes. It should be noted that ArtiSynth uses several simplifications compared with ANSYS, notably the use of semi-implicit integration and a lumped mass model.

The beam example (Figure 5(a)) consisted of a beam with dimensions 0.1×0.02×0.02 m divided uniformly into 8×4×4 hexahedral elements, with a density of 1040 kg/m 3 , Rayleigh damping coefficients of = 20/s and = 0.015 s, and a five-parameter Mooney Rivlin material with c 10 = 10 370, c 20 = 486 and c 01 = c 11 = c 02 = 0 Pa. Incompressibility in both system was enforced using a mixed u-P formulation [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF], and time integration was performed for 0.4 s using a 1 ms time step. To assess dynamic performance, the resulting z displacement and velocity of a reference node located in the middle of the free end were compared over time between ArtiSynth and ANSYS (Figure 6 (left)). The dynamic behavior was essentially identical: the resulting displacement error (relative to the maximum displacement) had maximum and average values of 0.3 and 0.08%. Likewise, the resulting velocity error (relative to the maximum velocity) had maximum and average values of 1.3 and 0.4%. We also determined the errors in total displacement and Von Mises stress for all the nodes in the final position: the maximum and average displacement errors (relative to the maximum displacement) were 0.06 and 0.04%, whereas the maximum and average Von Mises stress errors (relative to the maximum stress value) were 0.9 and 0.13%.

The cube example (Figure 5 (b)) used a cube with a width of 0.1 m in all directions and divided uniformly into 6×6×6 hexahedral elements, with a density of 1040 kg/m 3 , Rayleigh damping coefficients of = 20/s and = 0.015 s, and a five-parameter Mooney Rivlin material with c 10 = 1037, c 20 = 486, and c 01 = c 11 = c 02 = 0 Pa (identical to the material used for our tongue model described below). In addition to gravity, an immediate external load of -0.8 N was applied in the vertical direction to the nine nodes located in the middle of the top surface, resulting in the deformation shown in Figure 5. Incompressibility in both system was enforced using the B-bar method [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF], and the example was integrated for 0.2 s with a 1 ms time step. To assess dynamic performance, a reference node was selected in the middle of the top surface and its z displacements and velocities were compared (Figure 6 values of 22 and 1.4%. The large value for the velocity error occurred at the beginning, where ANSYS computed an unexpected initial upward velocity for the node. Compared with ANSYS, the ArtiSynth behavior was slightly more damped. For all nodes in the final position, the maximum and average displacement errors were 3.5 and 0.5%, whereas the maximum and average Von Mises stress errors were 5.6 and 0.5%. Much of this error was due to differences in the way ArtiSynth and ANSYS compute pressure for the B-bar method, resulting in different dilational displacements: in ArtiSynth the model compressed slightly, whereas in ANSYS it inflated slightly. These results help to demonstrate that our simulation approach is competitive with commercially available codes, in addition to be considerably more efficient: ArtiSynth was 20 and 10 times faster for the beam and cube examples, respectively.

APPLICATION TO JAW-TONGUE DYNAMICS

Building on our experience with the 3D jaw-hyoid [START_REF] Hannam | A dynamic model of jaw and hyoid biomechanics during chewing[END_REF] and the 3D tongue [START_REF] Buchaillard | A biomechanical model of cardinal vowel production: muscle activations and the impact of gravity on tongue positioning[END_REF] models, we introduce here the first 3D jaw-tongue-hyoid dynamical model taking into account full coupling between the FEM tongue model and the jaw-hyoid bony structures. As described in Section 3, we depend on the various components of Artisynth to provide dynamic simulations of interactions due to muscle forces of the jaw-tongue-hyoid complex and contact phenomena such as tongue-palate collisions.

A desire to better understand structure-function relationships and motor control strategies in upper airway movement has led a number of scientists to develop biomechanical models of individual structures in the human oral, pharyngeal, and laryngeal complex. For example, in the context of upper airway pathophysiology, biomechanical modeling enables an analysis of the structural causes of observed dysfunction. Biomedical applications typically require patient-specific models of the corresponding structures (jaw, lips, tongue, soft palate, larynx, and pharynx). Upper airway modeling has been applied to obstructive sleep apnea [START_REF] Huang | Use of computational modeling to predict responses to upper airway surgery in obstructive sleep apnea[END_REF][START_REF] Chouly | Simulation of the Retroglossal Fluid-Structure Interaction During Obstructive Sleep Apnea[END_REF] and dysphagia (disorders of swallowing) [START_REF] Lundy | Aspiration: cause and implications[END_REF], as well as planning and evaluating the consequences of maxillofacial surgery [START_REF] Chabanas | Patient specific finite element model of the face soft tissues for computer-assisted maxillofacial surgery[END_REF][START_REF] Luboz | Orbital and maxillofacial computer aided surgery: patient-specific finite element models to predict surgical outcomes[END_REF], segmental jaw surgery [START_REF] Stavness | Predicting muscle patterns for hemimandibulectomy models[END_REF][START_REF] Hannam | A comparison of simulated jaw dynamics in models of segmental mandibular resection versus resection with alloplastic reconstruction[END_REF], and glossectomy [START_REF] Buchaillard | Simulations of the consequences of tongue surgery on tongue mobility: implications for speech production in post-surgery conditions[END_REF]. Further, in the context of speech production, modeling the physics of speech articulators (jaw, lips, tongue, soft palate, and larynx) enables an evaluation of the motor control strategies underlying the production of speech gestures. Physical speech signals are indeed the consequence of the interaction between the mechanical structure of the articulatory apparatus and the forces that act on it. Part of the forces involved are centrally controlled, while other forces depend on interactions between articulators and on peripheral feedback. Thus, to study the control of speech production from the articulatory or acoustic signals, researchers have proposed to model separately the biomechanics of the articulators and the neurophysiology of the control [START_REF] Laboissiere | The control of multi-muscle systems: human jaw and hyoid movements[END_REF][START_REF] Perrier | Influences of tongue biomechanics on speech movements during the production of velar stop consonants: a modeling study[END_REF][START_REF] Perrier | The equilibrium point hypothesis and its application to speech motor control[END_REF].

Section 4.1 describes the jaw-tongue-hyoid model implemented in Artisynth, whereas Sections 4.2 and 4.3 report simulations of simple lingual and mandibular motor tasks, focusing on the dynamical interaction between the tongue soft structure and the jaw rigid body.

Jaw-tongue-hyoid model

For the dynamic tongue model, we implemented the published model by Buchaillard et al. [START_REF] Buchaillard | A biomechanical model of cardinal vowel production: muscle activations and the impact of gravity on tongue positioning[END_REF] in ArtiSynth, as pictured in Figure 7. The original tongue model was based on the anatomy of a single subject using CT data and developed in the ANSYS environment representing the tongue with hexahedral finite elements and hyperelastic properties. Thanks to a collaboration with Buchaillard and colleagues, we were able to obtain data for the 3D tongue mesh and description of the lingual muscular fibers. The mesh and muscle geometry were imported into the ArtiSynth environment, using a large-deformation FEM framework, hexahedral elements with a density of 1040 kg/m 3 , and a fifth-order incompressible Mooney-Rivlin material with c 10 = 1037, c 20 = 486, and c 01 = c 11 = c 02 = 0 Pa. The deviatoric potential energy ˆ of this material is hence

ˆ = c 10 ( ĪC -3)+c 20 ( ĪC -3) 2 ,
where ĪC is the first invariant of the deviatoric component on the left Cauchy-Green tensor [START_REF] Bonet | Nonlinear Continuum Mechanics for Finite Element Analysis[END_REF]. Incompressibility was implemented using a constraint-based mixed u-P formulation (as discussed in Section 3.2). Rayleigh damping coefficients of = 40 s -1 and = 0.03 s were used to provide critical damping within the range of modal frequency from 3 to 10 Hz, consistent with the reference tongue model [START_REF] Buchaillard | A biomechanical model of cardinal vowel production: muscle activations and the impact of gravity on tongue positioning[END_REF].

Muscles are represented by sets of elements and implemented with node-to-node fiber forces distributed throughout the muscle elements along the principal direction of action. While the material is isotropic, the muscle fiber components provided an approximation of the transversely anisotropic nature of skeletal muscle. We chose to use a straightforward model for muscle activation with fiber forces directly scaled by input activation, as opposed to the -model used in the original tongue model. Our aim was to quantify the coupling between tongue and jaw and not to work on the -model motor control assumptions provided by the equilibrium point hypothesis [START_REF] Feldman | Once more on the equilibrium-point hypothesis (lambda model) for motor control[END_REF]. Tissue stiffening due to muscle activation was also modeled in the same way as Buchaillard and colleagues, i.e. a linear increase of c 10 and c 20 values ranging between (1037 Pa, 486 Pa) at no activation and (10370 Pa, 4860 Pa) at full muscle activation. As in the original model, each muscle's force capacity was a function of its cross-sectional area (see Table I in [START_REF] Buchaillard | A biomechanical model of cardinal vowel production: muscle activations and the impact of gravity on tongue positioning[END_REF]), with force capacity distributed across fibers weighed by the volume of their surrounding elements. For the jaw-hyoid model, we started from a previously published model developed in ArtiSynth that has been used to simulate free jaw movements [START_REF] Stavness | An Integrated, Dynamic Jaw and Laryngeal Model Constructed From CT Data[END_REF] and chewing [START_REF] Hannam | A dynamic model of jaw and hyoid biomechanics during chewing[END_REF]. The model included rigid-bodies for the skull, jaw, and hyoid bone, point-to-point Hill-type actuators for the jaw muscles, constraint surfaces for the temporomandibular joints, and planar unilateral constraints for teeth contact. To conform the disparate morphology of the two models, we adapted skeletal and muscle geometry of jaw-hyoid model to fit CT data (shown in Figure 1(a)) for the subject on whom the Buchaillard tongue model was based. The 3D jaw, skull, and hyoid surface meshes were morphed with a non-elastic mesh-based registration algorithm [START_REF] Bucki | Finite element speaker-specific face model generation for the study of speech production[END_REF] to a 3D skull surface segmented from CT data. Symmetry was attained by mirroring the left-side of the registered meshes. The inertia of the jaw and hyoid were computed from new mesh shapes, assuming uniform density of 3600 and 2000 kg/m 3 for the jaw and hyoid, respectively. Jaw muscle origin and insertions points were adapted with the same non-elastic transformation as was applied to the surface meshes and were manually verified to correspond to plausible anatomical landmarks. We removed the point-topoint geniohyoid and mylohyoid muscles from the jaw-hyoid model as these were included in the tongue model. The anterior and posterior digastric muscles were connected to the hyoid bone with the digastric sling modeled as a pulley. The same Hill-type muscle dynamics were used from the original jaw-hyoid model with force capacity proportional to maximum cross-sectional area (see Table I in [START_REF] Stavness | Predicting muscle patterns for hemimandibulectomy models[END_REF]). The hyoid bone was attached to a fixed larynx with a linear translational/rotational spring representing the hyothyroid membrane and ligament. The spring stiffness was set to be We couple the dynamics of the jaw, tongue, and hyoid models by defining attachment constraints between the FEM nodes of the tongue and the jaw and hyoid rigid-bodies, as described in Section 3.5. Point to rigid-body attachments can be made at arbitrary locations and is not required to be coincident with the rigid-body surface mesh. The attachment points in the model are shown in Figure 7. Tongue-jaw attachments include the insertion of genioglossus and geniohyoid onto mandibular geniotubercle and the insertion of mylohyoid along mandibular mylohyoid ridge. Tongue-hyoid attachments include the entire region around the anterior-superior surface of the hyoid bone, including insertions of geniohyoid, mylohyoid, and hyoglossus muscles. The posterior medial surface of tongue is not attached allowing the base of tongue to move relative to the hyoid bone. The soft palate and palatoglossus muscle are not included in the current model. The resulting combined model is pictured in Figure 9.

Simulation descriptions

We chose to simulate a set of tasks similar to those reported for the jaw and tongue models in isolation, including free jaw movements [START_REF] Stavness | An Integrated, Dynamic Jaw and Laryngeal Model Constructed From CT Data[END_REF], unilateral chewing [START_REF] Hannam | A dynamic model of jaw and hyoid biomechanics during chewing[END_REF], and tongue movements in speech [START_REF] Buchaillard | A biomechanical model of cardinal vowel production: muscle activations and the impact of gravity on tongue positioning[END_REF]. All of the tasks, with the exception of unilateral chewing, involve simple input muscle activations so that the passive dynamic coupling effects can be better isolated. Our objective was to analyze the effect of dynamic coupling using muscle activation to drive the coupled jaw-tonguehyoid model and observe differences in the movement of the 'active' body as well as movement induced on the 'passive' body.

Jaw movement tasks.

Jaw movement tasks used jaw muscle activation as input and were performed both with the jaw-hyoid model alone and with the jaw-tongue-hyoid model. All jaw tasks involved a simple pattern of muscle input (rest, ramp-up, hold, ramp-down, rest) as illustrated I.

in Figure 10. The durations of these standardized jaw movements (600 ms) are consistent with an average chewing cycle. The muscle sets and activation amplitudes for the jaw tasks are summarized in Table I.

A nominal jaw movement task is rest posture: the equilibrium position of the model with no muscle activation and under downward gravity. In relaxed humans, the jaw typically rests with a 4-6 mm incisal separation. We expect slightly wider gape in jaw-tongue-hyoid model than in the jaw-hyoid model alone, although a majority of the tongue body rests on hyoid, and therefore should result in minimal jaw lowering. Static jaw clenching was simulated with bilateral activation of jaw closing muscles. With no tongue muscle activation, we expect the tongue to remain stationary within the mouth during tooth clenching.

Opening is simulated by bilateral activation of lateral pterygoids along with the anterior belly of digastric. We chose to simulate a moderate opening gape, with 15% activation in each muscle. Maximum jaw opening in humans is 50 mm on average although both backward head rotation and hyoid positioning become important at wide gape [START_REF] Muto | Positional change of the hyoid bone at maximal mouth opening[END_REF] and would complicate the task. We also simulated hinge-like jaw opening with activation of digastric alone. We expect reduced opening with jaw-tongue-hyoid model due to passive compression at floor-of-mouth. We also expect the tongue to passively protrude from the mouth during jaw opening [START_REF] Lowe | Tongue and jaw muscle activity in response to mandibular rotations in a sample of normal and anterior open-bite subjects[END_REF], and require contact handling (see Section 3.6) between the tongue tip and the lower teeth.

Protrusion is simulated by bilateral activation of lateral pterygoid muscles. The effect of jaw protrusion on the tongue has important implications for OSA as a common therapeutic device is a dental appliance used to advance the jaw and tongue in order to open airway [START_REF] Schmidt-Nowara | Oral appliances for the treatment of snoring and obstructive sleep apnea: a review[END_REF]. We expect reduced jaw protrusion in jaw-tongue-hyoid case as the lingual elastic connection between the jaw and hyoid should provide some resistance to jaw movement. Also we expect forward translation of the base of tongue. Right laterotrusion is simulated by activation of the left-side lateral pterygoid muscles. As for protrusion, we expect reduced lateral deviation in the laterotrusion task with jaw-tongue-hyoid model as compared with the jaw-hyoid model alone.

Unilateral chewing involves a complex pattern of jaw muscle activity. We simulated right-sided chewing movement using the same muscle activation patterns and food bolus that were reported for our original jaw model [START_REF] Hannam | A dynamic model of jaw and hyoid biomechanics during chewing[END_REF]. An elastic, spherical food bolus (10 mm in diameter) was positioned between the right first molars, which provided resistance during the closing phase of the chewing stroke and collapsed when the applied force exceeded 35 N. Since our adapted jaw-hyoid model has a different bone and muscle geometry, we expect that its chewing movement will be altered, but still plausible, as compared with the original jaw-hyoid model. We also expect that the chewing movement for the jaw-tongue-hyoid model will be significantly altered as the muscle patterns were previously tuned by Hannam et al. to a model without a tongue.

Tongue movement tasks.

Tongue movement tasks used tongue muscle forces to move and deform the tongue within the mouth. Tongue tasks involved a ramp-up, hold, and ramp-down pattern of muscle input similar to the jaw tasks, but with faster transitions (50 ms, see Figure 10) for consistency with the speed of speech movements. Tongue retraction was simulated by activating styloglossus using the activation trajectory shown in Figure 10. We expect that, with the jaw at rest, a retracted tongue posture should induce backward movement of jaw. Tongue-palate contact is an important movement for speech. Tongue tip contact with the anterior hard palate was simulated by activation of superior longitudinal, posterior genioglossus, and transverse muscles (see Table I). We stabilized the jaw with low-level (0.5%) activation of jaw closing muscles to maintain a nearly closed jaw posture. We expect that tongue-palate contact will induce a downward movement on the jaw, causing the jaw to open wider. We also performed a maximum tongue-palate pressure simulation by ramping the superior longitudinal and transverse muscles to maximum activation. The ability to generate tongue-to-palate pressure is an important component of healthy swallowing function and we expect that the model's maximum tongue-palate pressure will be comparable with recorded maximum tongue pressure measurements. The tongue-palate contact simulation required contact handling between the tongue tip and the hard palate surface mesh, as discussed in Section 3.6.

Simulation results

Videos of our simulation results and information for acquiring source code for the model are available on our web site: http://www.artisynth.org/orofacial/ijnmbe2010/.

Jaw-tongue-hyoid coupling.

We observed a number of interesting influences of dynamic coupling on the simulated jaw-tongue-hyoid movements. Jaw movements were altered by the presence of the tongue and tongue movements were observed to induce jaw movement. We found a resting jaw posture with 5.6 and 6.6 mm incisal gape for the jaw-hyoid and jaw-tongue-hyoid models, respectively. Static clenching with the jaw-tongue-hyoid model simulated correctly with the tongue remaining stationary in the mouth.

The results of the jaw opening and hinge-opening simulations are shown in Figure 11. In both cases, the amplitude of jaw opening is reduced in the jaw-tongue-hyoid model due to compression of the lower portion of the tongue between the jaw and hyoid. Figure 11 also illustrates the 3D model at the wide gape position showing that jaw opening does indeed cause passive forward protrusion of the tongue such that the tongue tip is resting on the lower teeth.

Protrusion and laterotrusion movements are shown in Figure 12. The amplitude of protrusion was reduced in the jaw-tongue-hyoid model, likely due to stretching of tongue tissue between the jaw and hyoid, but the amplitude of lateral deviation was comparable. Interestingly, the tongue also induced significant downward movement of the jaw during protrusion and laterotrusion. The tongue model is pulled forward during jaw protrusion transferring more of its weight from the hyoid bone to the jaw. Right-side chewing movement produced by applying muscle patterns tuned for a different jaw geometry provided a plausible tear-drop-shaped incisor movement in the current jaw-hyoid model as shown in Figure 13. The movement compares well with the one produced by the original jawhyoid model (see Figure 2 in [START_REF] Hannam | A dynamic model of jaw and hyoid biomechanics during chewing[END_REF]). The amplitude of the chewing envelope is reduced with the jaw-tongue-hyoid model, which is consistent with the jaw opening and laterotrusion simulations.

Simulation of tongue retraction in the mouth is pictured in Figure 14 (upper panels), along with a plot of incisor displacement (lower left panel). Styloglossus activation initially causes an upward and backward movement of the incisor as the condyles move up the articular slope, followed by a backward and downward displacement as the tongue retracts farther. The tongue retraction simulation also demonstrates the large range of tongue movement capable with the model and motivates the need for a large-deformation FEM approach.

Figure 15 shows the mid-sagittal position of jaw, tongue, and hyoid for the tongue-palate contact simulation. Tongue lifting and palate contact causes a downward jaw movement as expected. The right-most panel plots the incisor point displacement, which starts before tongue-palate contact (as denoted by the × on the plot). The initial downward jaw movement is caused by tongue muscle activation and it increases as force is applied between the tongue and the palate.

Comparison with published data.

The jaw-tongue-hyoid model has been assembled from previously reported reference jaw [START_REF] Hannam | A dynamic model of jaw and hyoid biomechanics during chewing[END_REF] and tongue models [START_REF] Buchaillard | A biomechanical model of cardinal vowel production: muscle activations and the impact of gravity on tongue positioning[END_REF]. Qualitative evaluation of the coupled model shows similar levels of force and range of movement as exhibited by each individual model; therefore, the new implementation of these models in the ArtiSynth framework compares well with the previously published versions. As a step toward the validation of the jaw-tongue-hyoid model, we have made preliminary comparisons to published data on tongue kinematics and forces. We used the tongue retraction simulation as a means to compare tongue velocity generated in the model with measured tongue movement during speech. Payan and Perrier [START_REF] Payan | Synthesis of VV sequences with a 2D biomechanical tongue model controlled by the Equilibrium Point Hypothesis[END_REF] recorded the movement of surface tongue points with an electromagnetic articulograph during tongue retraction in an [y-o] speech utterance for the speaker on whom the tongue model was based. Figure 16 shows the recorded velocity profile reproduced from [START_REF] Payan | Synthesis of VV sequences with a 2D biomechanical tongue model controlled by the Equilibrium Point Hypothesis[END_REF] as well as the simulated anterior-posterior velocity of one node on the tongue's upper surface. The backward movement has a bell-shaped velocity profile with a 220 mm/s peak value and a 150 ms duration, which compares well with the tongue velocity during the recorded tongue retraction in speech.

We used the maximum tongue-palate pressure provided by the simulation as a means to evaluate whether or not the tongue muscle forces incorporated into the model are within a plausible range. Direct measurement of muscle force in vivo is not possible; therefore, we rely on external force measurements as an indirect means to evaluate the resultant force generation capability with the jaw-tongue model. In particular, Utanohara et. al [START_REF] Utanohara | Standard values of maximum tongue pressure taken using newly developed disposable tongue pressure measurement device[END_REF] used a balloon-type disposable oral probe to measure tongue pressure by having subjects compress it onto the palate with maximum voluntary effort. The authors recorded pressures for a large subject pool (850 subjects) and reported 40. the magnitude of the contact constraints by the area of the contact contours (see Section 3.6 for discussion of collision detection and handling in ArtiSynth). Maximum tongue-palate pressure simulated with the model was 38.2 kPa, which compares well with the values reported by Utanohara et al. [START_REF] Utanohara | Standard values of maximum tongue pressure taken using newly developed disposable tongue pressure measurement device[END_REF].

DISCUSSION

Simulation results

The simulations reported here are a proof of concept demonstrating the effectiveness of our hardsoft tissue simulation and motivating the need for including dynamic coupling in simulations of jaw-tongue-hyoid movements. The reported simulations demonstrate that a wide range of movement, large forces, and large tissue deformations are possible within the current simulation framework. We have provided preliminary qualitative comparisons of tongue velocity and pressure to illustrate that the model behaves within a plausible range of human movement and force production. Further validation is planned, as additional quantitative data were recorded from the subject on which the model is based [START_REF] Badin | Three dimensional modelling of speech organs: articulatory data and models[END_REF], although transducing upper airway movement and forces is a significant challenge. Biomechanical simulations of the type reported here are indeed challenging to create. However, the interactivity afforded by the fast simulation times in ArtiSynth allows for reasonably fast refinement and exploration of the model's capability. Complex movements require precise coordination among a large number of muscle input degrees-of-freedom making trial-and-error tuning of muscle inputs to generate simulations tedious and likely over-fitted to a particular model's geometry. We believe that optimization-based inverse dynamics approaches are a potential direction in this regard, with which we had some experience in ArtiSynth [START_REF] Stavness | Predicting muscle patterns for hemimandibulectomy models[END_REF].

Our preliminary results also point to a few promising directions. The results of the chewing simulation show that the muscle patterns of [START_REF] Hannam | A dynamic model of jaw and hyoid biomechanics during chewing[END_REF] are applicable to a different skull morphology, as they produced a very similar chewing pattern, suggesting that they are not overly sensitive to skeletal or muscle geometries. The results also show that the addition of passive tongue tissue has a significant effect on free jaw movements and the chewing movement, due to both the passive elastic connection between the jaw and hyoid made by the tongue, especially in compression, e.g. reducing jaw opening, as well as the additional mass of the tongue body, particularly in jaw protrusion.

We plan to investigate what changes in muscle patterns are required to improve the chewing stroke in the jaw-tongue-hyoid model, e.g. increasing jaw opener muscle activation. It is noteworthy that the original muscle patterns reported by Hannam et al. were reported as being low in amplitude; therefore, muscle activation amplitudes could be increased and remain plausible. We also plan to investigate activating the tongue muscles in concert with the jaw muscles to simulate tongue movements [START_REF] Hiiemae | Tongue movements in feeding and speech[END_REF] and palate contact pressure patterns [START_REF] Ono | Evaluation of tongue motor biomechanics during swallowing-from oral feeding models to quantitative sensing methods[END_REF] during the chewing cycle. Interestingly, this would add a dynamically changing inertia as opposed to the current passive tongue mass. We are also in the process of integrating our jaw-tongue-hyoid model with an FEM model of the face [START_REF] Nazari | Simulation of dynamic orofacial movements using a constitutive law varying with muscle activation[END_REF] to analyze the effect of facial passive soft-tissue forces on jaw movement.

The tongue tissue is currently modeled as isotropic with point-to-point actuators embedded within the material to model anisotropic muscle fiber forces. We are incorporating a transverseisotropic material [START_REF] Weiss | Finite element implementation of incompressible, transversely isotropic hyperelasticity[END_REF] into ArtiSynth, which will provide a more realistic representation of skeletal muscle mechanics; however, in the current model we increase the stiffness of elements associated with muscle activation as an approximation of skeletal muscle stiffening. The tongue muscles are particularly challenging to model because multiple muscle groups, with different principal fiber directions, converge and inter-digitate within the tongue body. For this reason we are investigating a formulation allowing for the superposition of multiple transverse-isotropic materials with different principal directions.

Integration error

The simulations reported in Section 4.3 were computed using a second-order Newmark integrator (see Section 3.2) with a 5 ms integration step size. In order to assess the numerical error, we compared with results computed with a 1 ms integration step. Differences between 5 and 1 ms integration steps were found to be small. Figure 17 plots lower incisor point and tongue tip displacements for the jaw opening and tongue-palate contact simulations computed with 5 and 1 ms step sizes. Error was computed as the difference between the displacement trajectories relative to the maximum displacement. The jaw opening simulation showed a very small difference between the 5 and 1 ms integration step conditions: the incisor displacement error had average and maximum values of 0.5 and 1.4% and the tongue tip displacement error had average and maximum values of 0.7 and 1.6%. Larger discrepancies were found in the tongue-palate contact task, which involves significant contact situations: the incisor displacement error had average and maximum values of 1.1 and 3.8% and the tongue tip displacement error had average and maximum values of 2.0 and 10.0%. The integration error in the contact simulation is likely due to the discontinuous nature of contact as well as the spatial discretization of the palate/tongue surface meshes on which collisions are detected and responses are generated. Even in the worst case, although, the integration error remains smaller than the range of uncertainty associated with the model's assumptions and approximations.

Computational performance

The jaw-tongue-hyoid model described above has two free rigid-bodies and 946 FEM nodes for a total of 2505 degrees-of-freedom. With respect to (6), the addition of point-based tongue attachments, incompressibility, and jaw joints results in an M that is 2505×2505 and a G that is typically 2505×740 (varying somewhat depending on the number of FEM contacts). In addition, a few unilateral constraints are used to implement bite contact. Solution times for (6) using the method described in Section 3.4 vary from around 130 to 200 ms (depending on whether unilateral constraints are in play) on a 2.6 GHz Core 2 Duo processor. Overall solution time (including collision detection and all the steps of Section 3.7) for a 600 ms jaw opening task with a time step size of 5 ms is around 40 s and a 400 ms tongue retraction task with a time step of 10 ms is around 20 s. Much of this involves Java code that could be significantly optimized. This improves on the computation time reported in [START_REF] Buchaillard | A biomechanical model of cardinal vowel production: muscle activations and the impact of gravity on tongue positioning[END_REF], where a 100 ms task for the same FEM tongue model (with jaw/palate contact) required 40 min of computing time using ANSYS. There are a number of ways we believe these solution times can be further improved. Some preliminary experiments with iterative solvers (including ones built into Pardiso) suggest possible speed improvements of 5 to 10 times, although how best to employ iterative techniques with unilateral constraints is not completely clear. Iterative techniques are also amenable to implementation on GPUs. Alternatively, Pardiso itself can be run on multicore technology.

Stability and contact handling

Most of the simulation tasks described above could be computed stably using a time step between 5 and 10 ms. It should be possible to increase this step size using either an adaptive or fully implicit integration scheme, although the latter will increase the computational cost per step.

When stability problems occur, they are often related to incompressibility, which is usually enforced using bilateral constraints to restrict the divergence of the FEM elements. Incompressibility will occasionally also cause a small amount of hourglassing [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF]. Both problems can sometimes be avoided by instead using the B-bar method [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF], in which the amount of incompressibility is controlled through the material's bulk modulus parameter [START_REF] Bonet | Nonlinear Continuum Mechanics for Finite Element Analysis[END_REF]. This is supported in ArtiSynth.

There are also some areas in which the contact handling can be improved. First, as indicated in Section 3.6, it would be useful to provide edge-edge-type contacts to better handle situations with low-resolution elements. Second, contact stability can occasionally be a problem: correcting a contact's interpenetration can cause it to separate so that it is not detected as a contact on the subsequent time step, even if the relevant surfaces are being pressed together. The result can be a small amount of jitter as the surfaces continuously interpenetrate and separate. ArtiSynth controls this by limiting the position correction so as to maintain a small amount of interpenetration, but this approach is somewhat ad hoc.

The contact jitter problem is a part of the larger issue of contact-coherence, or minimizing the change in contact constraints between time steps. Improved contact-coherence will not only reduce jitter, but also the extent to which the constraint configuration in [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF] changes between steps, in turn enabling faster solutions.

CONCLUSION

We believe that the simulation framework provided by the ArtiSynth platform is a promising approach to furthering the state-of-the-art for dynamic biomechanical models with combined soft and rigid structures. Our system is intended to bridge the capabilities provided by traditional FEM-based techniques on the one hand and multibody techniques on the other. We have shown that the dynamic simulation of such systems is complicated by the need to combine the velocity solve, which includes a non-diagonal M, with both bilateral and unilateral constraints, resulting in a large sparse mixed LCP.

The size of this LCP can be reduced by removing the degrees-of-freedom associated with interbody attachments. Collisions involving deformable bodies can be handled by using the intersection contour between their surface meshes to isolate interpenetrating FEM nodes, which are in turn used to generate unilateral contact constraints. The LCP itself is solved by (1) treating some contact constraints as bilateral for the duration of each step, hence reducing the number of unilateral constraints, (2) solving the linear part of the LCP using a direct solver (Pardiso), and (3) using a Schur complement to reduce the remaining unilateral constraints into a (small) dense LCP. The resulting solver performs considerably faster than many existing commercial FEM packages. Contributing reasons for this may include (a) our use of semi-implicit time solves and (b) handling contact constraints in a step-wise, impulse-based fashion. With further improvements, possibly based on iterative and/or multicore techniques, we are optimistic about achieving completely interactive simulation times for models of the sort we have described.

Our model of coupled jaw-tongue-hyoid dynamics is a significant step toward a complete model of orofacial biomechanics. We simulated a number of tasks, including unilateral chewing and tongue motions used in speech production, to test the model and evaluate the importance of jaw-tongue coupling. The results for tongue velocity and maximum tongue-palate pressure were consistent with the literature data and showed that coupling is indeed important; in particular, the presence of passive tongue tissue between the jaw and hyoid appears to significantly resist jaw movement.

Modeling dynamic interactions of a jaw-tongue-hyoid model required all the features of our system. For example, the soft tissue deformation of the tongue coupled to movements of the jaw and the impact on the hyoid required the tight coupling of FEM and rigid body simulation. The contact of the tongue tip to the palate used the contact detection mechanisms. The joint constraints depended upon the bilateral constraint implementation. As the jaw-tongue-hyoid model dynamics incur high velocities and large forces, maintaining a stable solve was critical. Further, the high dimensionality of muscle control in the model requires fast simulation times to enable effective exploration and experimentation. The ArtiSynth system embeds all these necessary functions in an interactive simulation environment. Thus, ArtiSynth is an effective simulation environment for problems that require coupled soft-tissue/hard-tissue biomechanics that include bilateral and unilateral constraints.

Figure 1 .

 1 Figure 1. (a) The upper airway anatomy shown in a mid-sagittal CT image and (b) screen shot of ArtiSynth showing the jaw-tongue-hyoid model along with a timeline for controlling input/output data and a control panel for adjusting model properties.

Figure 3 .

 3 Figure3. Log-log plot showing factor times for A as a function of the number of nodes (which is proportional to the size of A) for a series of 3D FEM problems with a uniformly increasing node density.The slope of the line indicates a complexity of O(n1.7 ).

Figure 4 .

 4 Figure 4. Contact handling between two deformable models (with the topmost rendered as a wireframe), showing the intersection contour (blue) and the contact normals (black lines) of interpenetrating vertices from the upper mesh.

Figure 5 .Figure 6 .

 56 Figure 5. Examples used for validation, shown in their final positions with a rainbow plot of the resulting Von Mises stresses and the locations of their respective reference nodes. (a) Beam example, ArtiSynth; (b) cube example, ArtiSynth (maximum stress 2787 Pa); and (c) cube example, ANSYS (maximum stress 2661 Pa).

Figure 7 .

 7 Figure 7. Front, back, and sagittal cutaway views of tongue model. Attachment nodes are also shown for the jaw (front view, red spheres) and the hyoid bone (back view, blue spheres). Muscle groups include Genioglossus (GGA, blue; GGM, green; GGP, red), Styloglossus (STY, cyan), Geniohyoid (GH, magenta), Mylohyoid (MH, orange), Hyoglossus, (HG, red), Vertical (VERT, green), Transverse (TRANS, blue), Inferior Longitudinal (IL, cyan), and Superior Longitudinal (SL, magenta).

Figure 8 .

 8 Figure 8. Front, oblique, and sagittal cut-away views of the jaw-hyoid model in ArtiSynth. Jaw-hyoid muscles modeled include the Anterior, Middle, Posterior Temporalis (AT, MT, PT), Deep and Superficial Masseter (DM, SM), Medial Pterygoid (MP), Superior and Inferior Heads of the Lateral Pterygoid (SLP, ILP), and Posterior and Anterior Bellies of Digastric (PD, AD).

Figure 9 .

 9 Figure 9. Front, oblique, and sagittal cutaway views of the coupled jaw-tongue-hyoid model.

Figure 10 .

 10 Figure 10. Input muscle activation pattern for jaw tasks (left) and tongue tasks (right). The activation amplitude for each task is given in TableI.

Figure 11 .

 11 Figure 11. The jaw-tongue-hyoid model pictured during rest posture (REST) and at peak jaw opening (OPEN) with grid spacing of 10 mm. With the tongue muscles at rest, jaw opening induced tongue protrusion such that the tip rested on the lower teeth. Right-most plot shows the mid-sagittal incisor point movement for opening and hinge-opening with the jaw-hyoid model alone (JH, blue) and the jaw-tongue-hyoid model (JTH, red) with point spacing of 10 ms.

Figure 12 .Figure 13 .

 1213 Figure 12. Incisor point movement in the sagittal and frontal planes during jaw-muscle activated protrusion and right laterotrusion for jaw-hyoid model alone (JH, blue) and the jaw-tongue-hyoid model (JTH, red) with point spacing of 10 ms.

Figure 14 .

 14 Figure 14. The jaw-tongue-hyoid model pictured during rest posture (REST) and at tongue retraction with styloglossus activation (RETRACTED) with grid spacing of 10 mm. The right-most panel plots the incisor displacement in the mid-sagittal plane with point spacing of 10 ms.

Figure 15 .

 15 Figure 15. The jaw-tongue-hyoid model pictured during rest posture (REST) and with tongue lifted into contact with palate (PALATE-CONTACT) with grid spacing of 10 mm. The right-most panel plots the incisor point displacement, which starts before tongue-palate contact (× denotes the beginning of tongue-palate contact) with point spacing of 10 ms.

Figure 16 .

 16 Figure 16. Backward velocity of a point on the upper tongue surface. (DATA) Backward velocity profile for a point on the upper back surface of the tongue recorded with an electromagnetic tracking system an during [y-o] speech utterance (reproduced from [12] page 15, Figure 9(a)). (MODEL) Backward velocity profile generated by the model during tongue retraction simulation.

Figure 17 .

 17 Figure 17. Comparison of output displacements computed with 1 ms (solid lines) and 5 ms (dotted lines) integration steps for the tongue tip (top) and lower incisor (bottom) during jaw opening (left) and tongue-palate contact (right) simulations.

27:367-390 DOI

		10 4
	solve time (ms)	10 2
		10 2 10 0	10 3	10 4
			number of nodes
				: 10.1002/cnm

2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Biomed. Engng. 2011;

Table I .

 I Percentage muscle activation used in jaw-tongue-hyoid tasks.

				Jaw tasks
	Rest		CLR * -	AD -	ILP -	SLP -
	Clench		10	-	-	-
	pen		-	15	15	15
	Hinge-open		-	15	-	-
	Protrude Right-lateral		--	--	15 15 †	15 15 †
				Tongue Tasks
	Retract	CLR * -	SL -	TRANS -	GGP -	GGM -	STY 25
	Palate	0.5	30	30	60	30	-
	Max-palate	1	100	100	80	30	10

* Jaw closing muscles (AT, MT, PT, MP, SM, and DM). † Only left-sided muscles are activated.

Copyright ᭧ 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Biomed. Engng. 2011; 27:367-390 DOI: 10.1002/cnm

ACKNOWLEDGEMENTS

This work was funded by the Natural Sciences and Engineering Research Council (NSERC) of Canada. The authors thank Alan Hannam for providing the jaw model and his expertise, Gipsa-Lab, Grenoble for providing CT image data (P. Badin) and the tongue model geometry and constitutive law (S. Buchaillard, J.M. Gérard, P. Perrier, Y. Payan), and the ArtiSynth team.