Vol. 35 No. 3, July 1998

Pages 327-334## Mathematical modeling of normal pharyngeal bolus transport: A preliminary study

## Michael W. Chang, MD, PhD; Brigette Rosendall, MS; Bruce A. Finlayson, PhD

Departments of Rehabilitation Medicine and of Chemical Engineering, University of Washington, Seattle, WA 98195

This material is based upon work supported, in part, by a biomedical engineering research grant from the Whitaker Foundation.

Address all correspondence and requests for reprints to: Michael W. Chang, MD, PhD, University of Washington, Department of Rehabilitation Medicine, Box 356490, Seattle, WA 98195; email: mwc@u.washington.edu.

Abstract—Dysphagia (difficulty in swallowing) is a common clinical symptom associated with many diseases, such as stroke, multiple sclerosis, neuromuscular diseases, and cancer. Its complications include choking, aspiration, malnutrition, cachexia, and dehydration. The goal in dysphagia management is to provide adequate nutrition and hydration while minimizing the risk of choking and aspiration. It is important to advance the client toward oral feeding in a timely manner to enhance the recovery of swallowing function and preserve the quality of life. Current clinical assessments of dysphagia are limited in providing adequate guidelines for oral feeding.

Mathematical modeling of the fluid dynamics of pharyngeal bolus transport provides a unique opportunity for studying the physiology and pathophysiology of swallowing. Finite element analysis (FEA) is a special case of computational fluid dynamics (CFD). In CFD, the flow of a fluid in a space is modeled by covering the space with a grid and predicting how the fluid moves from grid point to grid point. FEA is capable of solving problems with complex geometries and free surfaces.

A preliminary pharyngeal model has been constructed using FEA. This model incorporates literature-reported, normal, anatomical data with time-dependent pharyngeal/upper esophageal sphincter (UES) wall motion obtained from videofluorography (VFG). This time-dependent wall motion can be implemented as a moving boundary condition in the model. Clinical kinematic data can be digitized from VFG studies to construct and test the mathematical model. The preliminary model demonstrates the feasibility of modeling pharyngeal bolus transport, which, to our knowledge, has not been attempted before. This model also addresses the need and the potential for CFD in understanding the physiology and pathophysiology of the pharyngeal phase of swallowing. Improvements of the model are underway. Combining the model with individualized clinical data should potentially improve the management of dysphagia.

Key words:mathematical modeling, pharyngeal bolus transport, swallowing.

## INTRODUCTION

Effective clinical management of dysphagia requires insight about the transport properties of the bolus, forces applied to the bolus, and bolus movement. Common clinical strategies involved in dysphagia management include dietary modification, such as selective bolus volume and consistency, and/or modification of forces applied to the bolus, such as changing posture or utilizing certain maneuvers. By changing bolus consistency, volume, and/or altering forces applied to the bolus, bolus movement can be modified to aid in the development of a safe swallowing strategy.

Current methods for the clinical evaluation of dysphagia are still quite limited. Videofluorography (VFG), widely used in clinical dysphagia studies, is limited because it reveals only the movement of the bolus. VFG describes little about the forces applied to the bolus that propel it through the pharynx into the esophagus. Clinical management of dysphagia requires insight into these forces and how they can be modified to achieve safe swallowing. Manometry permits pressure (force/area) measurement for pharyngeal bolus transport with sufficient spatial and temporal resolution. However, these pressure measurements can be misleading. A distinction must be made between the hydrodynamic pressure, measured when the bolus is passing the transducer, and the contact pressure, measured when the pharyngeal wall is contracting upon the transducer. Other shortcomings also exist for manometry in the pharynx. The presence of the catheter causes discomfort and inherently disrupts the flow of the bolus (1). Also, laryngeal elevation and shortening of the pharyngeal wall result in movement of the catheter and uncertainty as to where the pressures are actually being measured. Some of these shortcomings can be overcome by performing concurrent pressure measurement and VFG, known as manofluorography (MFG). MFG can potentially provide a good clinical evaluation of dysphagia, because it provides insight about bolus movement under external forces, but it is still not well suited for a thorough evaluation of the various clinical options available in dysphagia management. Also, MFG is not readily available in most major medical centers for the evaluation of dysphagia.

A mathematical model of the pharyngeal phase of swallowing can be used as an adjunct tool for clinical evaluation of the effects of changing bolus volume, bolus properties, and/or forces without exposing the individual to uncomfortable procedures and excessive radiation. Mathematical modeling can also make up for the limitations of current clinical tools for the evaluation of dysphagia. Hence, mathematical modeling of the fluid dynamics of pharyngeal bolus transport can provide another method for the evaluation of dysphagia.

Concepts of mathematical modeling have been widely used in science and engineering; however, the application of these concepts to clinical medicine is currently quite limited. Li and Brasseur have used mathematical modeling to successfully describe esophageal bolus transport under various constraints (2). To the best of our knowledge, a mathematical model has never been applied to the pharyngeal phase of swallowing. The inherent difficulty of the pharyngeal bolus transport is its complex and moving geometry.

In this work, we present a preliminary model for a normal bolus movement in the pharynx. This model uses pharyngeal and UES wall movements obtained from VFG (3,4). Transport properties for the standard e·z·hd (250 percent w/v) barium sulfate mixture are provided in the literature (5). This model demonstrates the feasibility of simulating pharyngeal bolus transport and potential clinical applications in dysphagia management. Improvements of the current model are underway to incorporate more clinical parameters into consideration (6).

Background

Mathematical Modeling

Mathematical modeling is a technique widely used in science and engineering; however, direct application in clinical medicine for client care is still quite limited. Steps involved in formulating a mathematical model applied to clinical medicine are summarized inFigure 1.

Figure 1.Three steps involved in mathematical modeling: the formulation of mathematical equations from clinical problems; the solution of mathematical equations; and the interpretation of clinical answers from mathematical solutions.When problems are encountered in clinical practice, it may be possible to pose them with a set of equations (Step 1). Once that is accomplished, the clinical problem, in essence, is in the mathematical world, and there are many mathematical techniques available to solve these equations (Step 2). With developments in numerical analysis and rapid advancement in computer technology, the task of solving equations has been greatly simplified. Once solutions are obtained, they must be interpreted in the context of the real world, namely clinical problems, to apply the solutions to client-care issues (Step 3).

Physics of Fluid Flow

Equations that describe the physics of fluid movement include the continuity equation and the equation of motion (7). The continuity equation, Equation 1, is a mathematical statement of the conservation of mass.

ais the divergence operator, r is the density of the fluid,vis the velocity vector, and t is time. The equation of motion, Equation 2, is a mathematical statement of the conservation of linear momentum, or Newton's Law applied to fluids.

p is the pressure,

gis vector representing the body force per unit mass (most commonly gravity), t , and is the stress tensor. Stress is a force per unit area.For a Newtonian fluid, the relevant transport properties of the bolus are the density (r ) and the viscosity (µ). The viscosity appears in the equation of motion through a constitutive relation which relates the stress response (t ) to the rate of deformation of the fluid (

d). For a Newtonian fluid, the stress is linearly proportional to the rate of deformation. The constitutive relation for a Newtonian fluid is given in Equation 3.

A non-Newtonian fluid is simply a fluid that does not follow Newtonian behavior: the stress is not linearly proportional to the rate of deformation. A bolus to be swallowed is often classified as a viscoelastic fluid. A viscoelastic fluid is a non-Newtonian fluid that exhibits behavior typical of both a fluid (viscous) and a solid (elastic).

Due to the transient nature of pharyngeal bolus transport, these equations must be solved in their time-dependent form. The first term on the right hand side of the Equation 2 represents inertial forces. It is shown below that the inertial force cannot be ignored in the pharyngeal phase of bolus transport.

Equations 1 and 2 describe the physics of the bolus movement under external forces. The constitutive equation relates bolus consistency to stress. Writing down these equations accomplishes Step 1 of the mathematical modeling (see

Figure 1). Solving these equations helps to reveal the relationship between bolus transport properties and the bolus movement under external forces. Computational fluid dynamics (CFD) is a method for the numerical solution of these equations.

Finite Element Analysis

Finite element analysis (FEA) is a special case of CFD, in which the flow of a fluid in a space is modeled by covering the space with a grid and then predicting how the fluid moves from grid point to grid point. CFD is used extensively to model the flow of air past aircraft and automobiles. The results are used in conjunction with wind tunnel experiments to deduce the effect of geometric changes on the wind resistance. Description of the FEA technique is beyond the scope of this article, but reviews of the method and its application to fluid dynamics can be found in some standard textbooks (8-12).To simulate bolus movement, it is necessary to model the flow of viscoelastic fluids, which is more complicated than modeling the flow of air. A viscoelastic fluid exhibits elasticity, like a rubber band, but flows, like water. When there is a free surface (a gas-liquid interface like that around a bolus), the elastic effects are particularly important. FEA is a computer-aided mathematical technique for obtaining approximate numerical solutions to the abstract equations of calculus that predict the response of bolus movement subjected to external influences. FEA is also capable of solving problems with complex geometry such as occurs when modeling the human oral-pharyngeal anatomy, and of modeling problems with a moving boundary such as peristalsis in the pharynx. Not only viscoelastic properties, but also surface properties of the food bolus can be modeled in FEA. No attempt has been encountered in the literature to model all these phenomena in pharyngeal bolus transport. Such problems with moving boundaries in intricate geometries are complex, and the FEA that solves them is necessarily complicated.

## METHOD

We present here a preliminary model for pharyngeal bolus transport using a relatively simple geometry. The major assumptions of this model are:

- Axisymmetric geometry. Cross-sectional areas are taken from literature. Radius is calculated assuming cylindrical geometry. Pharyngeal and UES diameter data obtained from Kahrilas et al. and Cook et al., respectively (3,4). Length data are estimated from Kahrilas et al. (13).
- Time zero is taken as the opening of the glossopalatal junction (GPJ). GPJ closure at 0.54 s.
- UES opening at 0.34 seconds, closure at 1.04 s.
- Laryngeal elevation and pharyngeal shortening are complete at the start of pharyngeal bolus transport.
- Minimum radius corresponding to pharyngeal closure is equal to 0.1 cm. Residual volume at final time equals to 0.16 mL.
- Bolus size is 20 mL. Pharynx is initially filled with 2.7 mL of fluid.
- Single-phase, incompressible, Newtonian fluid.
- Fluid properties for standard e·z·hd mix (250 percent w/v) from Li et al. (5).

r = 1.8 g/cm ^{3}

µ = 1.5 g/cm·sA good model for the transport of a liquid bolus from the mouth to the esophagus needs to incorporate the changing geometry of the bolus. A model of the esophageal phase of swallowing accomplished this by assuming a certain shape of the bolus that remains constant throughout the esophagus. This assumption is acceptable in the esophagus since its geometry is simple: at rest, the diameter is the same at every length, and while distended, the radius is constant over the majority of the wavelength. This is not true in the pharynx. The best description of the shape of any cross section of the pharynx is ellipsoidal. While the bolus is passing, at any length, the magnitude of either semi-axis changes rapidly with time. Kahrilas et al. and Cook et al. have digitized data from biplanar fluoroscopy studies to determine the shape of the pharynx at various locations as a function of time (3,4). Their interpolated data are shown in

Figure 2.The shapes of these curves vary significantly between lengths of the pharynx; thus, an algebraic expression for a simple wave form cannot be accurately derived. This is where using FEA is advantageous. The motion of the nodes on the boundary between the bolus and the wall of the pharynx can be prescribed as a function of time. This is what has been done in the preliminary model.

Figure 2.Preliminary mathematical model in axisymmetric geometry, with the centerline of symmetry at the right border. Initial and boundary conditions are prescribed at the glossopalatal junction (GPJ), pharyngeal wall, and the upper esophageal sphincter (UES).

Figure 2also shows a diagram of the finite element (FE) mesh set up for the preliminary model of the pharynx. The geometry is assumed axisymmetric with the cross-sectional area of the cylinder at any length equal to the cross-sectional area of the ellipse from the literature values. The area included in the model extends from the GPJ to the UES, with the distance between these two locations taken as 5 cm. This length assumes laryngeal elevation and pharyngeal shortening have been completed before the initiation of the pharyngeal bolus transport. The position of the nodes on the wall boundary are at the time when the bolus is held in the mouth prior to the initiation of the oral phase. At this time, the GPJ and the UES are closed. Complete closure of any portion of the pharynx cannot be modeled: the mesh cannot collapse onto itself. A minimum radius, corresponding to complete closure, is chosen as 1 mm. This value can be varied to obtain desired pressure signatures. A similar procedure was used by Li et al. in their model of the esophagus (2).The position of the nodes on the wall boundary are taken from literature sources for a 20 mL bolus at 4 locations shown in

Figure 2.The geometry at the entrance to the UES was taken as the same as the exit, but the wall movement is started at an earlier time. This shift in time was determined using a velocity found by integrating the cross-sectional area at the exit over time and dividing this result into the volume of the bolus. This gives a velocity at 18.4 cm/s. Since the length of the UES is 1 cm, this gives a shift in time of 0.0542 s. Data for the entire lumen was found using a b-spline interpolation through these 5 data points.The boundary conditions are also outlined in

Figure 2.A symmetry condition is applied at the centerline. The radial velocity of the nodes on the wall are prescribed as discussed above. The axial velocity of the nodes on the wall are 0 (no slip).At the GPJ, the tongue is forcing the fluid from the mouth into the pharynx (t=0). This driving force is modeled by specifying a force/area (normal stress) at the GPJ during the time it is open. The maximum opening is reached at 0.32 s. The GPJ reaches the minimum diameter again at 0.54 s and remains closed during the rest of the simulation. The stress normal (pressure) to the inlet is initially set at 22.5 mmHg and is reduced linearly to zero at time 0.32 s. This is the oral phase of the swallow. The bolus is being transported from the mouth to the pharynx. If it is assumed that there is a constant force pushing the bolus through the GPJ, the stress should decrease over this time, since the cross-sectional area is increasing. From time 0.32 to 0.54 s, the normal stress remains 0 (free). From time 0.54 to 1.08 s, the axial velocity is set to 0 (no flow).

From time 0 to 0.33 s, the UES is closed. It begins to open at 0.33 s and reaches the maximum opening around 0.6 s. The UES is closed again when the bolus has passed into the esophagus at time 1.04 s. The axial velocity at the outlet is 0 from time 0 to 0.33 s (no flow). The UES opens at time 0.33 s. From time 0.33 s to the end of the solution, the normal stress is 0 at the outlet (free).

The mathematical model is solved using FEA, with geometry information and boundary conditions provided to the computer program FIDAP&trade (v. 7.06, Fluid Dynamics International, Inc., Evanston, IL), running on a Silicon Graphics Indigo Workstation XZ 4000. Eighty quadrilateral elements were used in this model, with 369 nodes and 896 degrees of freedom.

## RESULTS AND DISCUSSION

Figure 3shows FE meshes representing pharyngeal chambers at various time steps from 0 to 1.04 s. All pharyngeal chambers are represented by axisymmetric geometries, with centerlines of symmetry at the right border. Clearly the axisymmetrical geometry represents an oversimplification of the pharyngeal wall motion. Movements of the anterior pharyngeal wall are currently "lumped" together with the posterior pharyngeal wall motion to represent a peristaltic pharyngeal wall action. Future models, incorporating our own VFG data, will include separate consideration for the anterior and posterior pharyngeal wall movements. The preliminary model incorporates tongue actions by specifying stresses at the GPJ (seeFigure 2). Future models will directly include tongue-based movements as boundary conditions.Figure 3demonstrates the ability of FEA to represent the changing geometry of the pharynx while the bolus is passing. Whereas this assumed axisymmetric geometry in the preliminary model is simple, the use of FEA will allow much more complex geometries, including the three-dimensional.

Figure 3.Finite element meshes representing the geometry of the pharyngeal chamber at various time steps. All pharyngeal chambers are represented by axisymmetric geometries, with centerlines of symmetry at right borders. Important time events are marked as: t=0: onset of GPJ opening; t=0.37 s: onset of UES opening; t=0.54 s: GPJ closure; t=1.04 s: UES closure.Including inertial effects in the solution means that nonlinear term in Equation 2 cannot be set to zero. These nonlinear terms result in a mathematical problem that cannot be solved analytically (on paper). The importance of inertial effects can be estimated by calculating the Reynolds number. The Reynolds number for peristaltic transport is defined below (14).

r and µ are bolus density and viscosity, respectively; v is the velocity; h is the average radius of the bolus, and a is the wavenumber (a=h /l ). The wavenumber (a ) is the ratio of the average radius of the bolus to the wavelength of peristalsis (length of the bolus). If the Reynolds number for peristalsis is less than unity, the inertial terms in the equation of motion can be neglected (14). Since the wave form for peristalsis in the pharynx is variable, it is difficult to calculate a Reynolds number defined by Equation 4. Time averaging the cross-sectional areas at various locations in the pharynx results in average radii of about 2.0 cm at the tongue base and the level of the valleculae. For a 20 mL bolus and assuming cylindrical geometry, this gives a peristaltic wavelength of approximately 6.4 cm. The velocity of lumenal closure in the pharynx measured from fluoroscopy studies is approximately 15 cm/s (14,15). The density and viscosity for e·z·hd barium sulfate suspension is 1.8 g/cm

^{3}and 1.45 g/cm;cds, respectively (5). Using these parameters and Equation 4, the Reynolds number is approximately 12; thus, inertial effects may be important.

Figure 4is a time history plot of the pressure at the UES for simulations with and without the inertial terms included in the equations of motion. Notice the pressure calculated with inertial effects can be more than twice the pressure calculated without them. This demonstrates that pharyngeal bolus transport requires a much more complicated mathematical model (our approach) than esophageal transport. Note that these pressure tracings should not look like those measured with manometry. As stated previously, manometry is a clinically important but technically difficult measurement. This difficulty arises from the presence of a relatively large magnitude of contact pressure obtained from direct pharyngeal wall contraction on the manometer. The contact pressure has little implication on bolus transport. It is the hydrodynamic, or intrabolus pressure that drives the bolus transport. Without concurrent VFG monitoring, no distinction of these two pressures can be possibly made. Other problems associated with the intrabolus pressure measurements are relative motion of the pharynx versus the transducer during the act of swallowing, and the pressence of the catheter within the pharynx may interefere with bolus transport. In our model, only the hydrodynamic pressures can be reported because we are modeling the fluid, not the pharyngeal walls. The pressures shown only represent the pressures that would be measured when the bolus is passing the transducer.

Figure 4.Comparison of pressure histories at the UES for models with and without inertial effects.To estimate the accuracy of the preliminary model, the cumulative volume through the UES was calculated.

Figure 5describes the volumetric flow rate of the bolus and cumulative bolus volume through the UES. Recall that data used to generate the geometry for that figure was for bolus volumes of 20 mL. The total volume passing through the UES from time 0.34 to 1.04 s is 18 mL. This is relatively good agreement, considering only four points on the mesh were obtained from VFG measurements and considering that flow rates were not specified in this model. Also note that the maximum volumetric flow rate corresponds to maximal UES opening at 0.6 s.

Figure 5.Bolus transport through the UES.

Figure 6shows the hydrodynamic pressure histories at three locations: the GPJ, 2 cm distal to GPJ, and the UES. Again note that these pressure tracings should not look like those measured with manometry. In using clinical manometry, both the hydrodynamic pressure and the contact pressure are measured. Without concurrent VFG monitoring, no distinction of these two pressures can be made. In our model, only the hydrodynamic pressures can be reported because we are modeling the fluid, not the pharyngeal walls. Certain similarities can be noted, however. The initial positive and large pressure peak is due to the tongue forcing the bolus into the pharynx. The large negative pressure peak at UES occurs when the UES opens. As the cross-sectional area expands, the fluid pressures are reduced. Conversely, when the cross-sectional area contracts, the pressure increases.

Figure 6.Hydrodynamic pressure history at the GPJ, 2 cm distal to GPJ (i.e., near laryngeal introitus), and the distal UES.Transphincter pressure drop is a technically difficult but clinically important measurement. Prior to UES opening, there is a drop of up to 17 mmHg; the pressure drop then decreases to slightly above zero after the opening of UES. The closest comparable data encountered in the literature at maximum UES opening (from proximal UES to esophagus at ~0.6 s) was 2.5 ± 2.8 mmHg, which is in reasonably good agreement with the result from the model (4).

## CONCLUSIONS

The purpose of this current model is to demonstrate the ability to model swallowing in the pharynx using a boundary with a prescribed motion. CFD using FEA is needed for the solution of the pharyngeal model, due to the complex geometry and inertial effects. Based on the preliminary model, it has been argued that the effects of inertia cannot be ignored in the pharyngeal phase of swallowing. This model demonstrates the need and the potential for CFD in understanding the physiology and pathophysiology of the pharyngeal phase of swallowing. It can potentially circumvent the need for clinical manometric measurements. The current model is limited, because only four locations on the moving surface were taken from experimental data reported in the literature. Future models will incorporate more data by using digitized information from VFG directly, considering two-phase flow and non-Newtonian behavior of the bolus (6). More importantly, future models will be constructed in three dimensions and will need to be tested against combined kinematic and manometric measurements obtained through MFG. Whereas this preliminary model does show the transport of nearly 20 mL of fluid through the UES, the accuracy of future models should be better. Investigation of future models will include sensitivity analysis of the effects of changing the boundary conditions at the GPJ and UES. The model will be applied to analyze pathological swallowing. Long-term goals of this work include studying variables such as bolus properties, bolus volume, head position, and maneuvers, which all influence pharyngeal bolus transport in nonimpaired and dysphagic individuals.

## ACKNOWLEDGMENTS

We thank Drs. Miller, Yorkston, Rohrmann, and Pope for their valuable discussions and suggestions.

## REFERENCES

- Brasseur JG, Dodds WJ. Interpretation of intraluminal manometric measurements in terms of swallowing mechanics. Dysphagia 1991;6(2):100-19.
- Li M, Brasseur JG, Dodds WJ. Analyses of normal and abnormal esophageal transport using computer simulations. Am J Physiol 1994;266(4 Pt 1):G525-43.
- Kahrilas PJ, Lin S, Logemann JA, Ergun GA, Facchini F. Deglutitive tongue action: volume accommodation and bolus propulsion. Gastroenterology 1993;104(1):152-62.
- Cook IJ, Dodds WJ. Dantas RO, Massey B, Kern MK, Lang IM et al. Opening mechanisms of the human upper esophageal sphincter. Am J Physiol 1989;257(5 Pt 1):G748-59.
- Li M, Brasseur JG, Kern MK, Dodds WJ. Viscosity measurements of barium sulfate mixtures for use in motility studies of the pharynx and esophagus. Dysphagia 1992;7(1):17-30.
- Rosendall BM. Mathematical modeling of the pharyngeal phase of swallowing (dissertation). Seattle, WA: University of Washington; 1996.
- Bird RB, Stewart WE, Lightfoot EN. Transport phenomena. New York: John Wiley & Sons; 1960.
- Pironneau O. Finite element method for fluids. Chichester, NY: Wiley; 1989.
- Gunzberger MD. Finite element methods for viscous incompressible flows: a guide to theory, practice, and algorithms. Boston: Academic Press; 1989.
- Schwartz HR. Finite element methods. London: Academic Press; 1988.
- Burnett DS. Finite element analysis. Reading, MA: Addison-Wesley Publishing Co.; 1987.
- Girault V, Pierre-Arnaud R. Finite element methods for Navier-Stokes equations. New York: Springer-Verlag; 1986.
- Kahrilas PJ, Dodds WJ, Hogan WJ. Effect of peristaltic dysfunction on esophageal volume clearance. Gastroenterology 1988;94(1):73-80.
- Shapiro AH, Jaffrin MY, Weinberg SL. Peristaltic pumping with long wavelengths at low Reynolds number. J Fluid Mech 1969;37(4):799-825.
- Jaffrin MY. Inertia and streamline curvature: effects on peristaltic pumping. Int J Engng Sci 1973;11:681-699.
- Kahrilas PJ, Dodds WJ, Dent J, Logemann JA, Shaker R. Upper esophageal sphincter function during deglutition. Gastroenterology 1988;95(1):52-62.

Submitted for publication July 22, 1996. Accepted in revised form May 16, 1997.

Go to TOP

^{↑}