A human stomach model is generated by using the “Virtual Population” model (Gosselin et al., 2014) – which has detailed, high-resolution anatomical full body models derived from magnetic resonance image data. The 3D model of the stomach is created by segmenting the stomach lumen from the Virtual Population 1.3 “Duke” model of this dataset (Figure 1A). The segmentation is done by using the software, Materialize MIMICS, and the smoothing and surface mesh generation are done by using the Materialize 3-matic.

The stomach wall motility, especially the antral contraction wave (ACW) is modeled following the previous study (Ferrua et al., 2011). The ACW wall kinematics are modeled as a pulse-wave propagating down toward the pylorus with its axis along the centerline of the antrum (see Figure 1B). The traveling cosine wave of the contraction strain, λ_a is modeled by

where s is the distance along the antrum centerline, n_p is the pulse count, V_p is the pulse propagation speed, T_p is the pulse interval, W_p is the pulse width, and h(s) is a spatial amplitude modulation function. In the present study, the parameter values are set to V_p = 2.3 mm/sec, T_p = 20 sec, W_p = 20 mm, and λ_a,max = 0.7 based on the published data (Ferrua et al., 2011). The centerline is divided into 4 segments as shown in Figure 1B; s₀–s₁: no contraction, s₁–s₂: the ACW amplitude grows, s₂–s₃: a constant amplitude ACW, and s₃–s₄: the ACW diminishes, and this is controlled by the amplitude modulation function, h(s) given by

The lumen wall contraction is then prescribed by

where x→w is the position vector of the lumen wall, x→w,0 is the lumen wall position at the initial, undeformed state, and r→w is the vector from the wall to the antrum centerline. Eqs. (1)–(3) generate peristaltic, traveling antral contraction wave, and the wall motility can be controlled by adjusting the model parameters. Note that there could be the substantial variability in the stomach size and shape as well as wall motility, which would affect the gastric flow strength and pattern.

The stomach model described in Section “Human stomach model” is immersed into a Cartesian volume grid to perform the fluid-structure interaction simulation using an in-house immersed boundary method (Mittal et al., 2008). The Cartesian domain (see Figure 2B) size used are 18 cm in x-, 8 cm in y-, and 12 cm in z-direction. This volume is discretized into a total of 360 × 160 × 240 (about 14 million) uniform Cartesian cells with a grid spacing of 0.5 mm. This resolution is based on the grid refinement study for the present model (see Appendix A).

In the current study, we focus on modeling drug dissolution accompanied by an aqueous solution whose properties are those of water. The equations governing the flow of the fluid are the Navier-Stokes equations:

where u→ is the flow velocity, p is the pressure, ρ is the density of the fluid, μ is the viscosity, and g→=-G⁢z^ is the gravitational acceleration vector and G = 9.8 m/s². The direction of gravity is determined from the orientation of the imaging data assuming a standing position. The incompressible Navier-Stokes equations (Eq. 4) are advanced in time using a fractional-step based algorithm (Chorin, 1969). A hybrid second-order upwind and central differencing is used for the convection term, and a second-order central-scheme is used for the diffusion terms. Time-integration is performed implicitly with the second-order Crank-Nicolson method and the non-linear convection term is updated iteratively. The pressure Poisson equation solver employs an efficient, immersed boundary based stabilized Bi-conjugate gradient (BiCGSTAB) method (Zhu et al., 2017).

The boundary conditions on the immersed boundary is applied using the sharp-interface immersed boundary method (Mittal et al., 2008). In this method, the surface of the immersed body is represented by an unstructured mesh, which consists of triangular elements. Using the surface mesh, the Cartesian cells in the fluid and solid regions are identified, and the flow equations are solved only in the fluid region. In the present immersed boundary method, no additional forcing term is used to apply the boundary conditions. The boundary conditions on the fluid-solid interface are imposed by prescribing the variable values on the nearest Cartesian cell in the solid region. The details about the present flow solver can be found in the multiple references (Mittal et al., 2008; Seo and Mittal, 2011; Zhu et al., 2017). This flow solver has been extensively validated for a variety of laminar/turbulent flows (Mittal et al., 2008), biological (Zheng et al., 2013a,b) and cardiovascular flows (Vedula et al., 2014), and fluid-structure interaction problems (Bhardwaj and Mittal, 2012; Shoele and Mittal, 2014), and employed for a wide range of studies of biomedical fluid flows (Seo and Mittal, 2013; Seo et al., 2014, 2016; Vedula et al., 2015, 2016; Zhu et al., 2018; Bailoor et al., 2021).

The time step size used in the simulation is Δt = 0.002 s, which resolves one ACW period (T_p) with 10,000 time points. Assuming a fully closed pylorus, no-slip and no-penetration boundary condition are applied on all stomach wall except the top part of the fundus. An open boundary condition with zero-gradient pressure and velocity, which allows small amount of inflow and outflow, is applied on the top part of the fundus to satisfy the mass conservation, and to model the accommodation of the fundus.

The shape of the pill is modeled as a cylinder with spherical cap ends (see Figure 2A). The length and diameter are set to 1 cm and 0.5 cm, respectively. For the composition, the pill is modeled as a non-disintegrating pill made of salicylic acid following the previous study (Bai and Armenante, 2009). In the present study, two different pill densities (ρ_p) are considered for the specific gravities, SG = ρ_p/ρ = 1.0 and 1.2. The pill motion in the stomach interacting with the fluid flow is simulated by solving a six-degrees of freedom (6DOF) model of Newton’s Second Law for its linear and angular motions. The governing equations are given by

where m is the mass and I is the moment-of-inertia of the pill, v and ω are translational and angular velocities, and F_f and M_f are the force and moment induced by the shear and pressure of the surrounding fluid, and F_c and M_c are ones by contact with the stomach wall. The fluid force and moment are computed by

where S_p is the surface of the pill, τ is the viscous shear stress, and r is the vector from the center-of-mass of the pill to the point on the surface.

For modeling contact between the pill and the stomach lumen, we employ a non-linear spring-based model. A schematic is shown in Figure 3. First, at each point on the pill surface, the contact point on the stomach wall is determined by finding the nearest point on the wall. Then the distance vector d→ is obtained, and the non-linear spring and damping force are calculated by using this distance vector:

where v_b is the pill surface velocity, v_w is the wall velocity, f_b is a non-linear function, δ_min is the minimum distance parameter, k_max is the spring constant, k_d is the damping constant, and f_c is the contact stress. δ_min controls the distance range of the contact force activation. Since f_b decays very rapidly for δ > δ_min, the contact force is only active very near the wall, i.e., δ ≤ δ_min. In this study, δ_min is set to 1.5 mm. Finally, the contact force and moment are obtained by

The parameters for the present contact model are k_max, k_d and δ_min. k_max represents the maximum contact stress, and to satisfy the no-penetration condition, a scaling analysis can show

where Δρ is the density difference between the body and the surrounding fluid, L_p is the length scale of the body which can be obtained by dividing the volume by the frontal area of the body, and C_k is an O(1) constant. In the present study, C_k is set to 2.5. The damping constant needs to be set to suppress unphysical rebound after the contact, and this can be represented by using a relaxation time scale, t_R,

In this study, we set the relaxation time scale relative to the simulation time step size using the constant, C_d, which is set to 5 in the current simulations.

The release of API from the pill and the transport of the API by the flow are simulated by solving the convection and diffusion equation for the mass transport:

where C_API is the concentration of the API and D is the mass diffusivity. For the salicylic acid pill, the Schmidt number for the mass diffusion, ν/D is about 400 (Bai and Armenante, 2009), where ν is the kinematic viscosity of the water. The dissolution of the drug is modeled by a surface erosion, and the mass transfer rate is given by

where m_d is the dissolved mass and n is the surface normal direction. The boundary condition for Eq. (11) on the pill surface is C_API = C_s, where C_s is the surface concentration or the saturated concentration, given by the solubility of the pill. For the salicylic acid, the solubility constant, ρ_p/C_s is about 400 (Bai and Armenante, 2009). The estimated surface erosion velocity due to mass transfer is about 2.5 × 10^–6 cm/s, and thus this effect is neglected. A zero mass-flux condition is applied on the stomach wall. Equation (11) is discretized by using a second-order finite difference scheme in space and advanced in time by a four-stage Runge-Kutta method. The same sharp interface immersed boundary method with the flow solver is also used for the mass transport equation.

Since the density ratio of the solid (pill) to the fluid medium (water) is close to 1, the coupled fluid and solid equations (Eqs. 4 and 5) are solved with an implicit scheme. The implicit coupling is done iteratively as follows:

where the function q represents the Navier-Stokes equations (Eq. 4), the function f is for the calculation of the forces (Eqs. 6–8), k is the iteration index, α is the relaxation parameter, and x→k and v→k are the position and velocity of the solid, respectively. In the iteration process, the flow velocity and pressure (u→k+1,pk+1) are updated first by solving the Navier-Stokes equations. The forces on the solid body (F→k+1) is then calculated using the updated flow variables. Finally, the velocity and position of the solid body (v→k+1,x→k+1) are updated using this force. For the stability and better convergence, we employed under-relaxation during the solid velocity update. The iteration continues till |x→k+1-x→k|<ε and then the solutions at k+1 are taken for the next physical time step n+1. In this study, we used ε = 0.0001 mm and α = 0.5.

The simulation is first performed without the pill to investigate the flow patterns in the stomach. The gastric flow patterns are Eulerian time averaged over one antral contraction cycle (T_p = 20 s) and shown in Figure 4. Figure 4A shows the streamlines of the time-averaged gastric flows. One can see the energetic flows consist of jet and recirculation in the antrum region, while the flows in the body and fundus are very weak. This is because the gastric flow is primarily driven by the antral contraction wave. The flow velocity contours on the cross-section in the antrum region are plotted in Figure 4C. It shows a strong jet at the center of the antrum, which is surrounded by a flow toward the pylorus. The jet at the center of the antrum is called the “retropulsive jet,” and has been reported in many previous studies (Pal et al., 2004; Ferrua and Singh, 2010; Ferrua et al., 2011; Trusov et al., 2016) as one of the characteristics of gastric flow during digestion. The strength of the retropulsive jet depends on the ACW strength and the gastric fluid properties (Pal et al., 2004; Alokaily et al., 2019). In the present study, the peak jet velocity is about 3 cm/s. The velocity contours on the cross section at the body region are also shown in Figure 4B the retropulsive flow toward the fundus is observed here. However, the overall velocity magnitude is O(10^–3) cm/s, which is very weak compared to the flow in the antrum region.

To investigate the time-dependent flow pattern in the antrum region and the development of the retropulsive jet, the snapshots of the flow velocity and relative pressure on the antrum cross-sectional plane (see Figure 2B) are plotted in Figure 5. At 0/4 T_p the ACW is propagating toward the pylorus and this traveling contraction with the closed pylorus increases the pressure downstream. The pressure drop across the ACW develops the retropulsive jet. As the ACW is propagating closer to the pylorus the downstream pressure is further increased and the retropulsive jet becomes stronger (1/4 T_p). The ACW is then diminishing as it reaches the pylorus, and this releases the downstream pressure and also weakens the jet (2/4 T_p). When the ACW reaches the pylorus the retropulsive jet disappears and no pressure drop across the antrum is observed (3/4 T_p). The process is repeated as the next ACW is coming and this generates the pulses of the jet (i.e., retropulsive jet). Since the retropulsive jet makes the strong and energetic flows in the antrum region, it should have a significant impact on the pill motion and drug dissolution.

In order to investigate the pill interacting with the gastric flow and its dissolution in the stomach, fluid-structure interaction simulations are performed employing the pill model described above. The simulation for one ACW cycle (20 s) took about 138,240 CPU h (240 h wall-clock time with 576 CPU cores) on the TACC Stampede 2 cluster.¹ In this study, we considered two different pill densities to study the effects of gravity as well. It is also important to note that the pill density can be modulated in the oral drug design to control the residence time of the pill in the GI system (Gupta et al., 2009). The first case, the specific gravity (SG) of the pill is set to 1 so that the pill is neutrally buoyant and the gravity does not play a role on the pill motion, and, for the second case, SG is 1.2, which is the actual specific gravity of salicylic acid pill (Bai and Armenante, 2009). The pill is initially placed at the entrance to the antrum region (see Figure 2A) and the simulations are performed for 60 s (3 ACW cycles).

The trajectories of the pill for 60 s obtained from the fluid-structure interaction simulations are plotted in Figure 6 for SG = 1.0 and 1.2. For the neutrally buoyant case (SG = 1.0), the pill floats around the entrance region of the antrum due to its interaction with the gastric flow. It is observed that the pill is always located upstream of the ACW because it continues to be pushed out by the retropulsive jet. The heavier pill (SG = 1.2), on the other hand, quickly settles down by the action of gravity (Figure 6B, 0–1), and displaces out of the core of the retropulsive jet. The pill itself is therefore not pushed out by the retropulsive jet but is pushed toward the pylorus directly by the antral wall contraction (Figure 6B, 1–3). Once the wall contraction is released, it is again settled down by the gravity (Figure 6B, 3–4). In this way, the pill is always located downstream of the ACW and near the pylorus.

The dissolution of the pill (surface erosion) is simulated by directly solving the convection and diffusion equation for the API concentration (Eq. 11) for the duration of 60 s, and the volumetric distributions of dissolved API concentration are shown in Figure 7 for SG = 1.0 and 1.2. Further details of the interactions between the gastric flow and the transport of API concentration are shown in Figure 8 for SG = 1.0 and 1.2. These figures show the velocity vectors and API concentration contours on the antrum cross-sectional plane.

For SG = 1.0, the pill is located superior to the ACW and dissolved there. The dissolved API is, however, entrained by the recirculating flow around the retropulsive jet and transported into the antrum. At the same time, the API is also transported toward the body of the stomach by the retropulsive jet. As a result, the API is distributed over the entire antrum except very near the pylorus. For SG = 1.2, the pill is located very near to the pylorus. The dissolved API concentration is then transported by the retropulsive jet to a region superior to the ACW and toward the body. The API is then mixed in the antrum by the retropulsive jet and the re-circulating flow. It is interesting to note that although the locations of the pill are very different for SG = 1.0 and 1.2, after a while (t = 50 s), the overall distributions of the API are very similar for both cases except very near the pylorus. This is because of high quality mixing driven by the ACW and the retropulsive jet.

The dissolution rate or mass transfer coefficient is one important metric to quantify drug dissolution (Bai and Armenante, 2009). For a non-disintegrating drug, the dissolution occurs primarily by the surface erosion. A well-known model of the dissolution is the Nernst-Brunner equation (Nernst, 1904):

where m_d is the dissolved mass, A_p is the surface area of the drug, D is the diffusion coefficient, δ_d is the apparent diffusion boundary layer thickness, C_s is the solubility, C_∞ is the bulk concentration in the medium, and k_m = D/δ_d is the mass transfer coefficient. The bulk concentration can be estimated by C_∞ = m_d/V, where _Vis total volume of the medium. At the early stage or for the case where the volume of the medium is much larger compared to the pill, C_∞ is much smaller than C_s and thus negligible, and the solution of (Eq. 14) yields the linear profile of the dissolved mass (Healy et al., 2002):

At the later stage or for the case where the volume of the medium is small, C_∞ = m_d/V has to be included and the solution exhibits an exponential profile.

In the present study, the dissolved mass in the stomach is monitored in time, and the profiles of the normalized dissolved mass are plotted in Figure 9 for both SG = 1.0 and 1.2. For the duration of initial 60 s, the time profiles of the dissolved mass follow the linear function (Eq. 15) as expected. The heavier pill (SG = 1.2) follows the linear trend better, while the neutrally buoyant pill (SG = 1.0) shows slight oscillatory deviations, of which the period is close to the ACW cycle (∼ 20 s). The mass transfer coefficient, k_m, can be evaluated from the slope of the best-fit line. The evaluated mass transfer coefficient is about 0.0011 cm/s for SG = 1.0 and 0.0012 cm/s for SG = 1.2. These values are apparently smaller than the ones measured by Bai and Armenante (2009) in the experiments with the USP-II device (k_m = 0.005–0.01 cm/s) (Bai and Armenante, 2009).

The mass transfer from the surface of the pill is governed by mass diffusion, and the mass transfer rate is directly proportional to the wall normal gradient of the concentration (Eq. 12). The fluid flow over the pill surface can increase the concentration gradient and thus, the mass transfer rate. The strength of the fluid flow over the surface is related to the wall shear stress or wall shear rate, and this suggests that the concentration gradient may be correlated to the wall shear rate. We can rationalize this by considering the diffusive mass flux to be analogous to the diffusive momentum flux (shear stress) from the classical Reynolds analogy in heat transfer (Geankoplis, 2003). In order to understand the effects of gastric flow on the mass transfer rate further, we investigated the correlation between the wall shear rate and the concentration gradient on the surface of the pill for both SG = 1.0 and 1.2.

The wall shear rate magnitude, WSR, and the normalized API concentration gradient, CG, are calculated on the pill surface by

where n denotes the direction normal to the pill surface. Note that the wall shear stress magnitude is given by μ × (WSR) and the local mass flux is given by DC_s.(CG). These quantities are averaged over time in a Lagrangian manner for the duration of the pill motion shown in Figure 6. The spatial distributions of the time averaged WSR and CG are plotted in Figure 10 for SG = 1.0 and 1.2. The average WSR for SG = 1.0 is higher on one side than the other. The overall average value (time and space) is 2.19 s^–1. For SG = 1.2, the overall average value of the WSR is 7.65 s^–1, which is noticeably higher than SG = 1.0. This indicates that the relative flow strength over the pill is stronger for SG = 1.2 than SG = 1.0. It is because the neutrally buoyant pill (SG = 1.0) is convected by the flow, and thus the relative velocity between the pill surface and the flow is small. The average WSR and CG for SG = 1.2 show the similar pattern: higher values at the sides and lower in the middle. Overall average value of CG for SG = 1.2 is 19.2 cm^–1 which is slightly higher than one for SG = 1.0, 19.0 cm^–1, and this explains the slightly higher mass transfer coefficient for SG = 1.2.

To check the local correlation between the WSR and CG on the pill surface, the spatial correlation coefficient between two quantities is computed by

where bar denotes spatial average, ϕ₁ = WSR, and ϕ₂ = CG. For SG = 1.0, the COR is -0.01, which means that two quantities are not correlated. For SG = 1.2, however, the WSR and CG are correlated well with COR = 0.7. The results show that the WSR (or wall shear stress) and the concentration gradient (or mass flux) are locally correlated only when the relative flow strength over the pill is substantially larger.

The API is eventually meant to be absorbed in the intestine and the flux of the API through the pylorus into the duodenum determine the rate of availability of the API in the intestines. Although the emptying through the pylorus is not included in the current mode, the emptying rate can be estimated by the API concentration or dissolved mass in the pyloric region, since the emptying rate is proportional to the concentration at the pylorus.

The dissolved mass of the pill in the pyloric region is obtained by integrating the API concentration over the volume within 2 cm from the pylorus (see Figure 11A), and the time profiles of the dissolved mass in the pyloric region are plotted in Figure 11B for SG = 1.0 and 1.2. Overall, the amount of dissolved mass in the pyloric region is larger for the heavier pill (SG = 1.2), primarily because the pill stays near the pylorus. For the neutrally buoyant pill (SG = 1.0), the dissolved mass arrives at the pyloric region after about 30 s and the peak value is about 3 times smaller than the SG = 1.2 case. For both cases, the profiles fluctuate with the ACW period (∼20 s) because of the retropulsive jet. The result indicates that the heavier pill should have the faster drug emptying rate.