Abstract
Balloon dilation catheters are often used to quantify the physiological state of peristaltic activity in tubular organs and comment on their ability to propel fluid which is important for healthy human function. To fully understand this system's behavior, we analyzed the effect of a solitary peristaltic wave on a fluid-filled elastic tube with closed ends. A reduced order model that predicts the resulting tube wall deformations, flow velocities, and pressure variations is presented. This simplified model is compared with detailed fluid–structure three-dimensional (3D) immersed boundary (IB) simulations of peristaltic pumping in tube walls made of hyperelastic material. The major dynamics observed in the 3D simulations were also displayed by our one-dimensional (1D) model under laminar flow conditions. Using the 1D model, several pumping regimes were investigated and presented in the form of a regime map that summarizes the system's response for a range of physiological conditions. Finally, the amount of work done during a peristaltic event in this configuration was defined and quantified. The variation of elastic energy and work done during pumping was found to have a unique signature for each regime. An extension of the 1D model is applied to enhance patient data collected by the device and find the work done for a typical esophageal peristaltic wave. This detailed characterization of the system's behavior aids in better interpreting the clinical data obtained from dilation catheters. Additionally, the pumping capacity of the esophagus can be quantified for comparative studies between disease groups.
1 Introduction
Peristaltic flow through cylinders and other geometries has been studied extensively since the mid 1960s [1–4]. The investigations have ranged from analyses of specified wall motion on Newtonian [1,5], non-Newtonian [6,7], and particulate fluids [8–10] to the effects of a prescribed forcing on the coupled fluid–structure system [11–13]. Both infinite and finite geometries [13,14] have been considered, and the quantitative effects of the channel geometry (cylindrical versus rectangular channels) on pumping characteristics have been established [3,15]. However, one problem that has received little attention is the effect of peristalsis on fluid-filled elastic tubes that are closed at both ends. In such a setting, fluid inside the tube can neither enter nor leave the system for the entire duration of peristalsis. This configuration is commonly found in long, slender balloon catheters used to evaluate mechanical properties and the contractile response of blood vessels (similar to catheters used during angioplasty). It is also found in functional lumen imaging probe (FLIP) catheters used to characterize the contractility of anatomical sphincters and esophageal motility.
A major assumption utilized to facilitate the mathematical analysis of peristaltic flow through a tube is that the problem is identical in both the fixed (lab) frame of reference and in the reference frame attached to the peristaltic wave. In the latter, the wave is stationary, and the shape of the tube walls does not change with time. When the tube length is finite, this assumption is not valid [14]. Thus, in order to build a simplified model for our problem, we turn to the approach taken to analyze flow driven by valveless pumping, which occurs in a finite, periodic domain, and modify the forcing method and boundary conditions to reflect the operating conditions for our problem. Interestingly, the configuration we aim to study first appears in Ref. [1] but is subsequently ignored. To the best of our knowledge, this configuration does not appear again in any other studies involving flow in deformable tubes.
The device we intend to focus our analysis on is the aforementioned FLIP [16]. It consists of a long flexible, hollow catheter, at the end of which is mounted a polyurethane bag [17]. The section of the catheter enclosed by the bag incorporates paired impedance planimetry sensors that can measure the regional cross-sectional area (CSA) along the bag's length as seen in Fig. 1. Various versions of the device exist where the bag length is either 8 or 16 cm; our focus is on the latter version. Also within the bag at the distal end, there is a pressure sensor that captures the fluid pressure at any given time. When deployed in a patient, the device is positioned such that most of the bag rests within the esophageal lumen. The bag is then incrementally filled with saline, and the esophageal contractile response is evaluated by monitoring the internal fluid pressure and CSAs along the length of the bag. When fully distended with saline, the bag diameter is 22 mm. The external end of the catheter is connected to a computer that stores the data collected by these sensors. The bag can be filled or drained by pumping fluid through the catheter, and fluid enters or leaves the bag through a small hole on the catheter. Figure 1 shows a simplified representation of the bag-catheter assembly and shows the locations of the various sensors. For additional details on the device's construction, data resolution, accuracy and frequency of data collection, see Refs. [16] and [18].

Details of the FLIP bag and catheter assembly (positioned within the esophagus). The proximal end of the bag points toward the subject's throat and the distal end, toward the stomach. The length of the bag is 16 cm, and diameter at full distension is 22 mm.
1.1 Operating Details and Simplifying Assumptions.
Representative FLIP data obtained from a subject can be visualized as shown in Fig. 2. During examination, the subject is in a supine position. The left panel shows the pressure and volume inside the bag as a function of time. It should be noted that the volume change is controlled by the physician conducting the procedure, and that the pressure change is a consequence of esophageal contractility in response to distension. The right panel of the figure shows the axisymmetric profile of the tube at the instant in time indicated by the black vertical line on the left panel. The dotted lines superimposed on the tube profile show the maximum and minimum diameters that can be accurately measured by the planimetry sensors. The esophageal contraction is highlighted with the three black bands where the CSA is the least. On the pressure curve, we see four peaks, indicating that four peristaltic waves have passed over the bag for the duration of this plot. The image of the bag profile is oriented such that peristalsis begins at the top and travels downward. Measurement of area is not continuous along the bag's length. Rather, area is measured at 16 locations at 1 cm intervals, and intermediate values are interpolated to construct the lumen profile [19].
![Visualization of clinical data collected from the FLIP. Graph on the left shows the bag volume and pressure variation over 40 s of time (10 Hz recording frequency). Figure on the right shows the profile of the esophageal lumen at the instant marked by the vertical black bar on the left graph (at recording number 200). Dotted lines show the range of measurable cross-sectional area (21–380 mm2, 5.2–22 mm diameter) [16].](https://asmedc.silverchair-cdn.com/asmedc/content_public/journal/biomechanical/143/7/10.1115_1.4050284/1/m_bio_143_07_071001_f002.png?Expires=1747750001&Signature=jB6SqlesWE71LoT1aKK5tLqu5hNz~KzvBmWXlOtQmvlO5H20FqpKouSSdu0Wi48mZrL6MVBVYXf-X26war~0tFoxMBfd~ZKy3DZAho5zF45KKJls-hLCY4F5izKAyYmoOS~ZKYjs1nBOZhc0Un3ixAWfpOZkMdzZtFB4f2m-VeaV4AhwnGewpfToTvpqYw6NkgaVxGEZ1ylhBlZ4-3DirqN40WoewEl2unAugORtvo8fXonUQgLbfaG3LovLAGETOHsT1PpVW63LVUqK4dsC5-xJrBqbKu1xQlxhcjXvLHzt05gVO5FJ4Ckp8vH3KgOpZEDotiP7qdsx5ZnUYAWYfQ__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Visualization of clinical data collected from the FLIP. Graph on the left shows the bag volume and pressure variation over 40 s of time (10 Hz recording frequency). Figure on the right shows the profile of the esophageal lumen at the instant marked by the vertical black bar on the left graph (at recording number 200). Dotted lines show the range of measurable cross-sectional area (21–380 , 5.2–22 mm diameter) [16].
![Visualization of clinical data collected from the FLIP. Graph on the left shows the bag volume and pressure variation over 40 s of time (10 Hz recording frequency). Figure on the right shows the profile of the esophageal lumen at the instant marked by the vertical black bar on the left graph (at recording number 200). Dotted lines show the range of measurable cross-sectional area (21–380 mm2, 5.2–22 mm diameter) [16].](https://asmedc.silverchair-cdn.com/asmedc/content_public/journal/biomechanical/143/7/10.1115_1.4050284/1/m_bio_143_07_071001_f002.png?Expires=1747750001&Signature=jB6SqlesWE71LoT1aKK5tLqu5hNz~KzvBmWXlOtQmvlO5H20FqpKouSSdu0Wi48mZrL6MVBVYXf-X26war~0tFoxMBfd~ZKy3DZAho5zF45KKJls-hLCY4F5izKAyYmoOS~ZKYjs1nBOZhc0Un3ixAWfpOZkMdzZtFB4f2m-VeaV4AhwnGewpfToTvpqYw6NkgaVxGEZ1ylhBlZ4-3DirqN40WoewEl2unAugORtvo8fXonUQgLbfaG3LovLAGETOHsT1PpVW63LVUqK4dsC5-xJrBqbKu1xQlxhcjXvLHzt05gVO5FJ4Ckp8vH3KgOpZEDotiP7qdsx5ZnUYAWYfQ__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Visualization of clinical data collected from the FLIP. Graph on the left shows the bag volume and pressure variation over 40 s of time (10 Hz recording frequency). Figure on the right shows the profile of the esophageal lumen at the instant marked by the vertical black bar on the left graph (at recording number 200). Dotted lines show the range of measurable cross-sectional area (21–380 , 5.2–22 mm diameter) [16].
With the device located in the esophagus, the passage of a peristaltic wave causes the (axisymmetric) profile of the bag and fluid pressure within it to change. During peristalsis, the esophageal wall may or may not approximate (i.e., come in contact with) the catheter. It is difficult to predict the relationship between flow rate and pressure drop across the contraction zone. The Reynolds number (Re) of the system is estimated to be 660 when using the following values: density , viscosity , tube diameter , and peristaltic wave speed [20]. If the flow inside is assumed to be similar to pipe flow at every location, then the flow is laminar. But the Reynolds number in the contraction zone is difficult to estimate, so the nature of flow at the neck is unknown. As such, we will analyze two flow types, (corresponding to low and high Re): (1) parabolic flow everywhere and (2) assuming that that a simple friction factor can be used to relate flow rate and pressure drop in the entire domain. As the subject is in a supine position, we assumed that the effect of gravity on the device is negligible. Depending on the positioning of the esophagus in the subject, the catheter can become slightly curved during measurement. The impedance planimetry sensors measure only local lumen area. Thus, curvature information is not available during clinical examination or subsequent analysis. However, fluoroscopic imaging [21] shows that the curvature is not significant; thus, we assume a straight geometry to model this problem. During the procedure, bag volume is changed from 0 to 70 mL in 10 mL steps. The analysis presented in this work considers peristaltic activity occurring at 30 or 40 mL of distension where the esophagus is gently stretched, and circumferential wall strains are not high. Due to the inflated nature of the bag, the lumen profile tends to be largely circular.
Our goal was to build a simple model to analyze this system and understand the relationships between the tube profile, internal bag pressure, esophageal wall stiffness, and the intensity of the peristaltic contraction. With sound mathematical foundations, the work done by the esophageal walls to pump fluid within the bag can be defined. This fast, simple model will eventually be deployed at point-of-care to estimate pressure distributions and work done in real-time. Armed with this knowledge, we can quantify peristaltic work of the esophagus, which may prove valuable in the evaluation of dysphagia.
2 Mathematical Details of the One-Dimensional Model
In the context of peristalsis-driven flow in the esophagus, fluid velocities are much smaller than the velocity at which a disturbance propagates in the esophageal wall. Thus, we expect no shocks or discontinuities in the problem, which allows us to use the nonconservative forms of the continuity and momentum equations. For smooth solutions, both conservative and nonconservative versions of the governing equations are equivalent and will lead to the same numerical solution.
Thus, the addition of this term helps stabilize the numerical solution. A similar manipulation was also used in Ref. [22] to eliminate p. In our analysis, the damping is kept small so as to ensure that the linear part of the tube law is the major contributor to fluid pressure for all regimes of operation. Note that A(x, t) and are not equal. The former is an unknown whose value must be found as part of the solution. Before proceeding, it is important to make note of the viscous shear stress term in the momentum equation. Several approximations can be made to write τR as a function of tube area and local fluid velocity. For our analysis, we choose to study the system when (1) the flow is parabolic everywhere with viscosity μ and, (2) a viscous term that uses a constant friction factor f is used to compute flow stresses. In the former case, the viscous resistance term is equal to and is for the latter.
2.1 Nonimensional Version of the Governing Equations.
where the following nondimensional numbers , , and emerge. Here, η is a nondimensional tube damping parameter, and β is a flow resistance coefficient. Eventually, the nondimensional pressure P in the momentum Eq. (9) will be replaced with the right-hand side of Eq. (10). Thus, the total number of coupled equations to be solved is two with α and U as the unknowns. When working with the friction factor version of the momentum equation, the nondimensional number for flow resistance β becomes and ψ remains unchanged, and the viscous resistance term in the momentum equation takes the form .
2.2 Boundary Conditions and Peristaltic Wave Input.
as Neumann-type boundary conditions for α at the tube ends. For nonzero damping, this BC is equivalent to changing the viscous resistance by an amount at the ends. This additional resistance at the boundary will have negligible effect for small η. With the additional regularizing terms and applied BCs, the solution differs from the hyperbolic version at a thin region near the boundary. In Sec. 2.3, using the Method of Manufactured Solutions, we show that the effect of damping on the overall solution is negligible. It should be noted that more complicated forms of the tube law can be used which account for longitudinal curvature, bending, and tension [29,30] in the tube wall. Such a tube law can include a double derivative, Axx or a fourth order derivative term Axxxx of the tube area. Under this setting, additional boundary conditions for area and pressure can be applied in a straightforward manner without leading to inconsistencies in the problem definition. The inclusion of these terms will lead to higher order spatial derivatives into the governing equations. However, in our current analysis, we choose to work with the simplest version of the tube law as the esophagus' material properties and behavior for these deformation modes are unknown.
The nondimensional width of the peristaltic wave is denoted by with W denoting the dimensional contraction width. As mentioned earlier, the CSA at the tube ends in the device tapers off to a point on the catheter, so a small correction is made to the activation term to account for the reduced CSA at the ends. The final form of the activation wave is plotted in Fig. 3. The value of θ0 specified represents the reference area of the tube at the strongest part of the contraction. Physiologically speaking, θ0 is a parameter that determines the strength of the contraction. It represents the ratio by which the zero-stress luminal area is reduced (i.e., the reference luminal area A0) due to muscular contraction. The smaller the value of θ0, the smaller is the activated reference area and the greater is the contraction intensity. When , there is no contraction and the system is fully at rest. Since the velocity scale was chosen to be the same as the speed of the peristaltic wave, the nondimensional wave speed is 1.0.
2.3 Numerical Solution of the System of Equations.
A snapshot of the manufactured tube shape and fluid velocity are shown in Fig. 4. The left segment of the tube has reached the maximum area of at . At this instant, the right segment of the tube is beginning to expand ,and the left segment is contracting, causing fluid to flow into the expanding section of the tube. This corresponds to a positive fluid velocity as shown by the dotted curve. Eventually, the right segment of the tube reaches its maximum area and begins to contract and the cycle continues with fluid traveling from one segment of the tube to the other in a periodic manner.
Using the nondimensional form of the momentum Eq. (9), the source term corresponding to the chosen solutions was found with values of β and ψ set to 103 and ω = 10. Using the regularized version of the governing equations given by Eqs. (19) and (20), along with and , the system was then solved to obtain α and U. Relative tolerance for the solver was set to . A comparison of the manufactured solution to the computed solution is given in Figs. 5 and 6. The agreement between the chosen and numerically computed solutions was found to be satisfactory, indicating that the equations were being solved correctly, and that the effect of the regularizing terms on the overall solution was negligible.

Comparison of fluid velocity obtained from the numerical solution of the regularized equations to the manufactured solution given by Eq. (23) at

Comparison of fluid velocity obtained from the numerical solution of the regularized equations to the manufactured solution given by Eq. (23) at
3 Dynamics Displayed by the One-Dimensional Model
3.1 Comparison With Three-Dimensional Immersed Boundary Simulations.
Before proceeding with a discussion on the regimes and computing peristaltic work, we set out to further validate our one-dimensional (1D) model with detailed fluid–structure interaction simulations using the immersed boundary (IB) method. It has already been stated that the nature of flow within the device is not definitively known to be either laminar or fully developed turbulent flow. As explained in Sec. 1.1, the Reynolds number for flow in the tube is estimated to be around 660, indicating that the flow is likely to be laminar. For healthy individuals, the value of ψ is of the order 100 and β is of the order 1000 based on known values of peristaltic wave speed [31], tube stiffness coefficient K [26], and material properties of water. The material parameters for the three-dimensional (3D) IB simulation were chosen such that they were equivalent to ψ = 100, β = 621, and , and flow was observed to have a parabolic profile everywhere within the tube. Following peristaltic activation, the resulting tube shapes and fluid velocity variations obtained from both models were then compared.
The geometry of the structure used in the IB simulations consists of a long cylindrical tube with open ends. The two ends are closed using “caps” that are simply short flat cylinders of specified thickness. The entire structure is then meshed using hexahedral elements to take advantage of the Adaptive Anisotropic Quadrature that was developed by Kou et al., in Ref. [32] and implemented in the open-source immersed boundary framework IBAMR [33]. The structure is represented using finite elements, and the combined fluid–structure equations are solved on a Cartesian grid using the Finite Volume Method. Mathematical details on the interaction equations, discretization of the Lagrangian and Eulerian domains, and temporal schemes for solving the coupled system can be found in Ref. [34]. Computations were performed on the Northwestern University “Quest” cluster and on SDSC Comet. The mechanical properties of the tube consist of two materials. Just like the esophagus, the structure in our IB simulations consists of fibers embedded in an isotropic matrix material. The matrix component is represented by a simple Neo-Hookean strain energy function, and the fiber's effect is implemented using a bilinear strain energy function that computes the stress based on strain in the circumferential direction. The combined effect of these layers leads to an approximately linear relationship between pressure and area, corresponding to the tube law used in the 1D model. Additional details on the material properties and strain energy functions corresponding to each layer can be found in Refs. [20,32,35]. Similar to the activation method used in Ref. [32], the peristaltic wave is applied as a controlled reduction in the fiber rest length along the tube. This reduction leads to a greater generation of circumferential stresses in the tube at the activation location, the consequence of which is a reduction of the tube area resembling peristaltic contraction.
The results of our simulations are summarized in Figs. 7 and 8. The figures show the resultant tube shapes generated by the peristaltic contraction wave as it travels from left to right along the tube length for both models. As the contraction moves forward, it displaces fluid, causing an increase in tube area to accommodate the excess fluid. When the pressure in the segment ahead of the contraction becomes high enough, the contraction opens up, allowing fluid accumulated in this segment to flow in the opposite direction and refill the segment of the tube behind the contraction. Figures 7(e) and 7(f) show the tube profile and fluid velocity after the wave has passed. At each instant, we see satisfactory agreement between the velocity variations observed in the 3D immersed boundary simulations and the nondimensional velocity profiles predicted by the 1D model. The tube shapes are also observed to match well with the shape of the structure computed by the 3D simulations during and after the activation wave has traveled over the tube's length.

Comparing tube shapes and velocity fields computed from the 1D model with results from equivalent 3D immersed boundary simulations. At each instant, the top figure shows the tube mesh along with the internal velocity fields along an axial slice obtained from the IB simulations and the bottom figure shows the corresponding variations predicted by the 1D model. In the latter, we plot the top and bottom profile of the tube using the nondimensional radius , and the variation of velocity U in the tube at each instant as a function of χ. The dotted line represents the tube's axis. Values corresponding to are on the left axis and values for U are on the right axis. Negative values of U indicate fluid is traveling in the left direction. Contraction wave began at τ = 0 and leaves the domain at τ = 1. Instants 5 and 6 correspond to and , respectively. Figure 8 compares areas and the area-averaged axial velocity with U from the 1D model for instant 2. (a) Instant 1, (b) instant 2, (c) instant 3, (d) instant 4, (e) instant 5, and (f) instant 6.

Comparing tube shapes and velocity fields computed from the 1D model with results from equivalent 3D immersed boundary simulations. At each instant, the top figure shows the tube mesh along with the internal velocity fields along an axial slice obtained from the IB simulations and the bottom figure shows the corresponding variations predicted by the 1D model. In the latter, we plot the top and bottom profile of the tube using the nondimensional radius , and the variation of velocity U in the tube at each instant as a function of χ. The dotted line represents the tube's axis. Values corresponding to are on the left axis and values for U are on the right axis. Negative values of U indicate fluid is traveling in the left direction. Contraction wave began at τ = 0 and leaves the domain at τ = 1. Instants 5 and 6 correspond to and , respectively. Figure 8 compares areas and the area-averaged axial velocity with U from the 1D model for instant 2. (a) Instant 1, (b) instant 2, (c) instant 3, (d) instant 4, (e) instant 5, and (f) instant 6.

Quantitative comparison of nondimensional tube area and fluid velocity obtained from 3D IB simulations with the numerical solution of the 1D system of equations. This snapshot corresponds to Instant 2, shown in Fig. 7 where the contraction is approximately at .

Quantitative comparison of nondimensional tube area and fluid velocity obtained from 3D IB simulations with the numerical solution of the 1D system of equations. This snapshot corresponds to Instant 2, shown in Fig. 7 where the contraction is approximately at .
Encouraged by this agreement between the 1D model and the 3D immersed boundary simulations, we continued the analysis of the system by probing the 1D model for a wide range of operating conditions. Although the agreement was found only for parabolic flow, we also investigated the system's response when the viscous term is friction-factor based. In Sec. 4, we present the results for both of these scenarios in a “unified” regime map. One of the primary goals for this comparison with the 3D IB simulations was to check if the introduction of the tube damping and smoothing terms significantly affected the response of the system. Figures 7 and 8 show that that the introduction of the regularizing terms not only leaves the physics of the system unharmed but also aids in stabilizing the numerical solution which directly leads to more rapid solution times.
It should be noted that the 3D immersed boundary simulations did not show any modes of asymmetric collapse or solid–solid contact between the tube walls for the operating parameters that were tested. The cross section of the tube remained circular across its entire length for the entire duration of the run. The speed at which a disturbance propagates in the system is given by [29] and is approximately 1 m/s for the esophagus. The peristaltic wave speed and internal fluid velocities are of the order of 1 cm/s indicating that the flow speeds are always well below the critical value needed for collapse or choking. At the location of the contraction, the tube is not passive and the there is no collapse of the structure. The reduction in area is solely due to active contraction brought upon by a controlled reduction of the rest length of the fibers in the tube's material model. It should also be emphasized that the coupled fluid–structure simulations we have developed are perfectly capable of capturing any event of collapse if there was a possibility for such an event to occur. However, for the range of parameters that are relevant to the operating conditions of this device, collapse is improbable due to the tube's geometry and the significant amount of fluid in the bag (30 mL or greater). Another thing to note is that as the FLIP only measures area, comparison of simulation results with realistic three-dimensional structures is not possible at the moment. In the future, we aim to obtain detailed geometric information about the esophageal structure from fluoroscopy and 4D-MRI imaging [36] to compare numerical results with experimentally measured values.
3.2 Types of Tube Deformations Observed During Peristalsis.
With all the parts of the 1D model in place, we can now investigate the system's response to the applied activation as a function of the operating parameters, θ0, ψ, and β. Changing these values is equivalent to investigating the effect of tube wall stiffness, wave speed, contraction strength, fluid density and flow resistance on the tube wall deformation and internal flow patterns during peristalsis. After studying the system's behavior for a broad range of operating values, we found four distinct, physiologically relevant patterns of peristaltic pumping based on the way the tube walls and the fluid inside respond to the applied activation. These regimes are shown in Fig. 9. The progression of these regimes from numbers 1 to 4 can be interpreted simply as the response of the system as fluid viscosity is continually increasing. It should be noted that the fourth regime displays very little deformation of tube area. The deformation dynamics for this regime are unremarkable, and the tube shape remains quite similar to the shape of the initial condition for the entire duration.

Tube deformations corresponding to regimes occurring for physiologically relevant parameter ranges. Definitions used to classify a solution as belonging to a particular regime are given in Table 1. (a) Regime 1, (b) regime 2, (c) regime 3, and (d) regime 4.

Tube deformations corresponding to regimes occurring for physiologically relevant parameter ranges. Definitions used to classify a solution as belonging to a particular regime are given in Table 1. (a) Regime 1, (b) regime 2, (c) regime 3, and (d) regime 4.
The occurrence of these regimes is dictated by a competition between elastic forces generated due to deformation of the tube wall and resistance to flow through the narrowest part of the contraction. A relatively stiff tube wall resists deformation and forces the fluid that was displaced due to the advancing wave to flow back through the contraction. In this scenario, we see the walls deform as shown in regime 1. An equivalent way of looking at this is to assume a low resistance to flow through the contraction. The energy required to expand the tube walls is significantly greater than the energy required to overcome viscous resistance across the contraction, and thus, the pumping is almost quasi-steady with the tube walls having the same diameter on either side of the contraction. On the other hand, if the resistance to flow is high, it is favorable for the system to expand the tube walls to accommodate the fluid that is being displaced by the advancing peristaltic wave. At moderate flow resistances, this exact situation occurs when the stiffness of the tube walls is low. The tube deformation pattern in this scenario is denoted as regime 2.
Regime 3 has been investigated in great detail by Takagi and Balmforth [11]. In their work, the lubrication approximation is used in an infinitely long domain with open ends. In our configuration, we see the formation of a “blister” as predicted by their model when the resistance to flow is high, and the tube walls are fairly compliant. The formation of a blister is due to the fact that the amount of time it takes for the fluid to redistribute is comparable to the time it takes for the wave to travel along the tube length. Regime 4 occurs when the resistance to flow is much higher than that in regime 3. In this case, the amount of time it takes for the fluid to “move” and respond to the peristaltic contraction is much higher than the amount of time it takes for the wave to travel over a section of the tube. The fluid velocity in this case is extremely low leading to a small value in as well. In this region of operation, the fluid essentially behaves as a semisolid entity, and the contraction causes only a small local deformation as it travels forward.
When the value of ψ (which is the inverse of the Mach number squared) is small, the wave-like nature of the system begins to emerge. In this region, the speed of the peristaltic wave is comparable to the speed at which a disturbance in the system travels at. For the esophageal wall, this disturbance speed is of the order of 1 m/s, but the peristaltic wave speed rarely exceeds 10 cm/s. As such, we will not investigate the dynamics of the system at these unphysiologic wave speeds.
4 System Behavior Represented as a Regime Map
In Sec. 3, four patterns of observed tube wall deformation depending on the system's operating conditions were presented (Fig. 9). Our goal was to combine the various operating parameters to summarize the response of the system in the form of a regime map. We wished to find the least number of parameters that could be used to describe the response of the system to the peristaltic contraction wave. The Buckingham-Pi approach leads to the several nondimensional numbers already shown above. Plotting the system's response for each of these quantities is cumbersome. To that effect, we utilized the wave-frame approach that has often been employed by other researchers to analyze peristaltic flow to find the most convenient combination of operating parameters to visualize the system's response. Detailed derivation of these parameters is provided in the Supplemental Materials on the ASME Digital Collection.
The system's response was investigated for a wide range of physiologically relevant parameters. Values of ψ were selected from the following set: {0.1, 0.25, 0.6, 1, 2.4, 5, 10, 24, 50, 100} and β had the following values: {0.1, 0.25, 1, 5, 10, 50, 100, 250}. Two values of the initial condition αIC = {1.25, 2.0}, corresponding to the volume of fluid in the tube, were chosen. The following values for contraction intensity θ0 were selected: {0.05, 0.1, 0.15, 0.2}. The nondimensional width w was set to 0.25 for all runs. All possible combinations of these parameters lead to 640 cases each for the parabolic and friction factor-version of the momentum equation. The damping parameters and solver tolerance were set to the same values that were used during verification using the Method of Manufactured Solutions. Volume of fluid in the tube was monitored as a measure of solution accuracy. Percent change in tube volume across all runs was of the order of . The tube shape computed from each of these runs was analyzed when the contraction was at the halfway point, i.e., at . At this instant, the ratio was computed along with the maximum tube area achieved () and area of the contraction (Ac). Based on visual inspection of their values for several cases, the tube shape was then classified as regime 1, 2, or 3/4 using the definitions provided in Table 1. Complete parametric information and the associated regime for each case plotted in the regime map are provided in a data sheet available in the Supplemental Materials on the ASME Digital Collection.
Definitions used to classify each solution as belonging to a particular regime
Regime 1 | |
Regime 2 | |
Regime 3 | |
Regime 4 |
Regime 1 | |
Regime 2 | |
Regime 3 | |
Regime 4 |
Note: areas A1 and A2 are the distal and proximal areas (Fig. 10), and is the area of the bulge (observed in regime 3 only) at its widest point. For regimes 1 and 2, . The area at the contraction is denoted by Ac.
Armed with this simplification, we present the occurrence of each regime as a set of points on a versus plot. Figure 11 shows the occurrence of the various regimes at different values of the two s. We combine all the regime points from the two different flow types and present them in a single plot. We see that the regime data points from both flow types fall within the general vicinity of each other. The set of points belonging to each regime is represented by a specific marker.

Regime map showing system's response for a variety of input conditions for both parabolic and friction factor-based flow resistance terms. The parameters and are computed using Eqs. (24)–(26) for each data point. The x-axis represents a cumulative stiffness parameter, and the y-axis represents a cumulative flow resistance parameter. Legend: regime 1 (), regime 2 (
), regimes 3, 4 (
).

Regime map showing system's response for a variety of input conditions for both parabolic and friction factor-based flow resistance terms. The parameters and are computed using Eqs. (24)–(26) for each data point. The x-axis represents a cumulative stiffness parameter, and the y-axis represents a cumulative flow resistance parameter. Legend: regime 1 (), regime 2 (
), regimes 3, 4 (
).
The wave frame approach that was taken to find the combination of parameters does not account for the formation of a bulge or a blister. As both regimes 3 and 4 occur in the region of high fluid viscosity, they are combined into a single entity for plotting purposes. From this plot, we clearly observe the presence of a general boundary between each set of regime points. The boundary represents the change of the system from one solution to another. This analysis can be used as a starting point to quantitatively derive the location of the boundaries between the regimes shown in Fig. 11. In addition to the parameters mentioned above, it has been observed that there is active relaxation ahead of, and stiffening behind, the contraction wave that can change the effective material properties of the tube to improve bolus transport through the esophagus [37,38]. In a future work, we will be incorporating these nonuniform changes in effective material properties into our analysis and study the effect on the observed regimes and quantitatively predict the regime boundaries as a function of all physiologically relevant parameters.
With the help of the regime map, one can predict the shape assumed by the tube during peristalsis by simply computing the relevant and from the given operating conditions and finding the region which these points belong to. The regime map also helps in understanding the change in tube shape that would occur as one of the operating conditions is changed. For instance, as the tube stiffness increases, the system will deform in such a way that the area ratio tends to 1.0. Increasing the fill volume or decreasing the peristaltic wave velocity leads to a similar outcome. When it comes to increasing the fluid's viscosity, we see a transition from regime 1 to regime 2 as predicted by the regime map when the value of is increased. When the area contraction factor θ is reduced, i.e., the amount of wave “squeezing” is increased, we again observe the system transition from regime 1 to regime 2.
Each of these regimes corresponds to a specific deformation pattern observed using the FLIP device in different subjects and/or diseases. For instance, assume that the device is calibrated, and the operating conditions are benchmarked in such a way that it shows pumping patterns corresponding to regime 2 in a healthy individual. Under the same operating conditions, a patient with peristaltic dysfunction will display pumping patterns corresponding to either regime 1 or regime 3/4. Depending on the specific behavior displayed, the regime map helps us identify the cause of the abnormality. For instance, a patient displaying regime 1 might have a stiffer esophagus due to fibrosis, and a patient with regime 4 has dysphagia due to ineffective peristalsis. Thus, the quantification of the device's behavior in the form of these regime maps directly assists in better interpreting the various shapes seen in patients during the FLIP procedure.
5 Defining and Quantifying Pumping Effort
Using this breakdown of fluid pressure, we can compute the values of each of the terms in Eq. (33) and visualize the variation of work done or energy dissipated over a single peristaltic event for each regime. The above analysis can be repeated for the momentum equation with the friction factor stress term, but the resulting work variations show no qualitative differences between the two scenarios.
5.1 Work Curves From the Reduced Order Model.
which shows that for large values of β and ψ (as is the case for flows relevant to esophageal peristalsis), the term corresponding to rate of change of fluid kinetic energy is very small compared to the other terms in the energy balance equation.

Representative work curves for regimes that show nontrivial tube deformation. Each regime shows a unique signature for variation of work during peristalsis. At each instant, the sum of passive and viscous work is equal to the active work done. (a) Regime 1, (b) regime 2, and (c) regime 3.
The work curves for each of these regimes show unique identifying features that confirm what is visually observed through the tube deformation patterns. For regime 1, once the wave has created a zone of reduced area, the active work goes into overcoming the viscous resistance, and the potential energy stored in the walls remains unchanged. Upon approaching the right boundary, the contraction leaves the domain causing the tube to relax, and the stored elastic energy is recovered. For regime 2, however, we observe a gradual rise in stored elastic energy with time indicating that the activation wave is continuously doing work on the tube wall. Unlike regime 1, when the wave approaches the end and allows the tube walls to relax, the stored energy is lost via fluid dissipation. This is observed by following the viscous work curve which sharply rises around the mark to meet the active work curve. It is also interesting to note that in spite of large tube wall deformations associated with this regime, the majority of active work done is still lost via viscous dissipation. Work curves for regime 3 reflect the low wall deformations observed from the tube shape plots, which only show a bulge ahead of the contraction. The passive work done is small compared to the active work done, which is almost entirely lost via viscous dissipation. When the contraction approaches the end, it acts on the bulge in the tube ahead of it as it cannot go any further. Contraction of this bulge then reduces its area and causes the accumulated fluid to flow in the opposite direction. This leads to the sharp rise in viscous work at τ = 1 as seen in Fig. 12(c). For regime 4, there is a simple, linear increase in viscous dissipation with time. Passive work done is zero for all time due to negligible wall deformation.
This detailed look into the variation of the active, passive, and viscous work done during a peristaltic event gives us the tools to identify what a healthy pumping wave looks like and the conditions under which one might observe reasonable tube wall deformations. A significant change in the tube area for this configuration due to peristalsis (as seen with the healthy peristaltic wave in Fig. 2) indicates that the wave has some ability to move fluid forward in a normal setting, which involves bolus transport following normal swallowing.
5.2 Work Curves From Functional Lumen Imaging Probe Patient Data.
In Sec. 5.1, we have defined and provided a strong mathematical foundation for computing work done during peristalsis in the current configuration and the various sinks that consume the energy generated by muscular activation. The work calculations were performed in the context of the reduced order model presented in Sec. 2. In this section, we wish to extend the utility of the model to compute work curves using data obtained from the FLIP device when used in human subjects. The goal was to get a sense of the magnitude of peristaltic work done in normal subjects and to see what characteristics are observed in work curves obtained from patient data.
It is clear that the 1D model in Sec. 2 takes as input the activation wave and predicts the resultant tube wall deformation and the associated fluid velocity and pressure fields. The FLIP data on the other hand contains the variation of tube wall deformation as a function of time and the value of pressure at the (approximate) location . Simply speaking, and p(L, t) are already known. The values of fluid velocity and pressure at all the other locations along the tube are unknown. The device measures the lumen profile and stores diameter readings from each of the 16 sensors (as seen in Fig. 1). The recorded diameters are used to compute lumen area after accounting for the catheter's presence and are corrected to ensure volume conservation in the tube. Equation (1) is then used to calculate U(x, t) from planimetry data A(x, t) collected by the device. With the available values of A(x, t) and U(x, t), the momentum Eq. (2) is used to compute the pressure gradient. Using p(L, t), the pressure is then found as a function of x and t for the entire duration of peristalsis. Implementation details of this step using Eqs. (8) and (9) are provided in Ref. [39].
With these steps in place, we isolate a contraction wave to analyze and generate work curves for. An example of this can be seen in Fig. 2. A window of readings is chosen that includes a single pressure peak. At the start of computation, the wave begins to travel over the tube length and results in a pressure rise. The contraction travels along the esophagus as time goes on and by the last reading the wave has finished traveling over the entire tube length and the computation ends. In Fig. 13, we show the work curves and the pressures at (referred to as the proximal pressure) and χ = 1 (referred to as the distal pressure) as a function of time for a typical peristaltic contraction. Again, we emphasize that the former is predicted by our model, and the latter is measured by the device and is applied to the model to fix the pressure levels inside the tube.

Estimated proximal pressure variation and work computations from FLIP data corresponding to a single peristaltic contraction occurring at 30 mL in a healthy control. The esophagus does not return to its exact starting shape at the conclusion of peristalsis leading to a small nonzero value of passive work at the end. (a) Contraction 1—pressure variations and (b) contraction 1—work versus time curves.

Estimated proximal pressure variation and work computations from FLIP data corresponding to a single peristaltic contraction occurring at 30 mL in a healthy control. The esophagus does not return to its exact starting shape at the conclusion of peristalsis leading to a small nonzero value of passive work at the end. (a) Contraction 1—pressure variations and (b) contraction 1—work versus time curves.
The first thing to note from the analysis of the patient's secondary peristaltic contractions is the pressure predicted by the model at the proximal (left) end of the tube at . When the wave is incoming, the pressure in the entire system rises, but once it passes over the left end, there is a sharp drop in pressure. This is due to the temporary displacement of fluid due to the peristaltic pumping action of the wave leading to the tube having reduced area at the left end. The displaced fluid then stretches the walls of the distal end of the tube at χ = 1 and this is marked by the continuous rise in pressure at that location. Once the wave has passed, the fluid accumulated at the end flows back and the tube attains a uniform shape and pressures at both locations equalize at the end of wave travel. This development of a proximal pressure trough has also been observed in a FLIP prototype where there was an additional pressure sensor near the proximal end. Due to inaccuracy in the sensor, the readings have not been plotted. Qualitatively speaking, however, the readings showed a drop in pressure which is supported by our calculations.
The active work is then estimated by subtracting the passive work from the total fluid pressure work. These estimated values of passive and active work are then plotted in the work curves shown in Fig. 13(b). The passive work is again seen to be smaller than the magnitudes of active and viscous work. One observation from these curves is that the total magnitude of work done for a single healthy secondary peristaltic contraction is of the order of 10 millijoules (mJ). Studies conducted using traction force sensors have been able to quantify the propulsive force generated by the esophagus during peristalsis [31]. On the lower end of the spectrum, the esophagus was able to generate 0.11 N of force. During FLIP examinations, around 12 cm of the bag lies within the segment of the esophagus that undergoes peristaltic contraction. Using these values of force and distance, a rough estimate of the work done was found to be 13.2 mJ, which is of the same order of magnitude as predicted by our model. In a separate clinical study, we compute the work done for several controls and patients belonging to four disease groups to understand differences in peristaltic work done [40].
6 Model Limitations
In this section, we discuss a few limitations of the 1D model developed to study the FLIP device and its response to peristaltic contraction.
During particularly strong peristaltic activity, circular muscle contraction can obliterate the esophageal lumen and completely cut off the proximal and distal segments of the tube from each other. In this scenario, both velocity U and area α in the contraction region are zero, and the model becomes singular. However, in a majority of peristaltic contractions studied with the FLIP, flow is observed to occur between chambers indicating nonzero α and U in this region.
Another avenue of exploration is the effect of nonlinear tube laws on the quantity of work done during peristalsis. The complex fiber architecture of the esophageal wall is one reason for this nonlinear behavior. The other reason is the polyurethane bag used in the device. For bag volumes beyond 40 mL, there is greater stretching of the bag walls. When fully stretched, they sharply resist further dilation. This sudden increase in stiffness might have a significant effect on the estimate of work done. The esophagus' material properties vary along its length as well. The exact variation of these properties is unknown, so as a first step, we assumed a simple linear tube law to estimate work. Additionally, tube laws that incorporate viscoelasticity will also significantly increase the amount of work done as the energy lost within the walls due to hysteresis will be considered.
Finally, the peristaltic activation wave does not have a constant velocity. The wave can speed up or even completely stop depending on the amount of obstruction it senses as a form of feedback regulation. In our configuration, this can become quite important as the ends are closed, and continuous pumping will lead to a larger obstruction which can alter the speed and intensity of contraction in subjects, directly affecting the work done by the wave. Another detail that has not been considered is the precise positioning of the bag within the esophagus. When deployed by the physician, part of the bag straddles the Esophagogastric junction, and the tip lies within the stomach. The Esophagogastric junction is a physiologically complex region that behaves quite differently compared to the rest of the esophagus. In our work, we have assumed that the entire bag rests within the esophagus.
7 Concluding Remarks
In this work, we have analyzed a hitherto unstudied configuration of fluid flow in an elastic tube. A simple reduced-order model is presented that is able to predict fluid flow and the resultant tube wall deformation when a peristaltic wave passes over a closed cylindrical tube. The model has been compared to detailed three-dimensional immersed boundary simulations and is shown to have satisfactory agreement with them. The system's response under this set of operating conditions has been thoroughly quantified and visualized as a regime map.
Finally, the system's utility is demonstrated by applying it to enhance the data collected with balloon dilation catheters. Paired with this tool, the device can now be thought of as having multiple pressure and velocity sensors housed on the catheter. With the help of these quantities, an appropriate mathematical foundation was laid to quantify the peristaltic work done. This step has led us to estimate the work done by the esophagus during a healthy contraction under these operating conditions. The device can now be used to find work done in other peristaltic waves and comment on their capability to propel fluid.
The analysis has a few limitations due to the inherent complexity of the combined organ-device system. However, the model presented in this work provides a foundation on which future studies can be based and used to further our understanding of esophageal pumping processes.
Acknowledgment
This research was supported in part through the computational resources and staff contributions provided for the Quest high performance computing facility at Northwestern University which is jointly supported by the Office of the Provost, the Office for Research, and Northwestern University Information Technology.
This work also used the Extreme Science and Engineering Discovery Environment (XSEDE) cluster Comet, at the San Diego Supercomputer Center (SDSC) through allocation TG-ASC170023, which is supported by National Science Foundation Grant No. ACI-1548562 [41].
Funding Data
National Institutes of Health (NIDDK Grant Nos. R01 DK079902 and P01 DK117824; Funder ID: 10.13039/100000062).
National Science Foundation (OAC Grant Nos. 1450374 and 1931372; Funder ID: 10.13039/100000105).
Disclosures
Dustin A. Carlson, Peter J. Kahrilas, and John E. Pandolfino hold shared intellectual property rights and ownership surrounding FLIP panometry systems, methods, and apparatus with Medtronic, Inc. No conflict of interest exists for Shashank Acharya, Sourav Halder, Wenjun Kou and Neelesh A. Patankar.