Inflammatory bowel disease (IBD) is a chronic inflammatory disease that primarily affects the gastrointestinal tract. IBD includes two main classifications: ulcerative colitis and Crohn's disease. Ulcerative colitis mainly affects the colon and rectum, leading to mucosal inflammation and ulcerative lesions in these areas. In contrast, Crohn's disease has the ability to damage any part of the gastrointestinal tract, extending from the oral cavity to the anal area, often causing transmural inflammation of the intestinal wall.

The manifestations of IBD generally encompass abdominal discomfort, diarrhea (occasionally accompanied by hematochezia), weight reduction, fatigue, and nutritional deficiencies. These clinical manifestations may exhibit variability, sometimes correlating with periods of disease exacerbation (flare-up phase), and at other times during remission, yet both phases can precipitate substantial distress for patients. Regrettably, the precise etiology of IBD remains inadequately elucidated and is posited to be intricately associated with genetic predispositions, aberrations in immune responses, and environmental influences. Presently, a definitive cure is elusive; however, symptom management and enhancement of patients' quality of life can be attained through pharmacological interventions, modifications in lifestyle, and surgical procedures.

Consequently, we have formulated a therapeutic strategy aimed at treating diseases via the secretion of pharmaceuticals through engineered probiotics, and have established mathematical foundations rooted in fluid dynamics, biostatistics, and pharmacokinetics, thereby providing comprehensive elucidations and feasibility analyses pertaining to the entire continuum of drug dissemination post-administration, the proliferation dynamics of engineered microorganisms within the host, the therapeutic mechanisms for disease management, and ultimately, the modulation or activation of the "self-destruction" pathway of the engineered probiotics. The comprehensive model is segmented into several distinct components as delineated below:

1. Model 1 ：To closely replicate the actual medication intake scenario of patients, it is assumed that the pills are taken with water on an empty stomach. The digestion process of the drug after entering the gastrointestinal tract is then simulated using fluid dynamics. This includes modeling the motion of the pill in the stomach, its transit through the intestines, the release of E. coli, and the flow and retention of E. coli in the intestines. These analyses provide a deeper understanding of the pill’s behavior within the digestive system and its influence on E. coli release and function.

Figure 1. Flow chart

2. Model 2：The proliferation of engineered bacteria in the gastrointestinal tract: To simulate the survival and reproduction of Escherichia coli in the intestines of patients with inflammatory bowel disease (IBD), we focused on a 5 cm segment of the transverse colon and abstracted it as a cylindrical model. The model not only takes into account the flow of the non-Newtonian fluid but also factors in how E. coli concentration changes over time due to the fluid’s shear forces, the specific pH levels typical in the intestines of IBD patients, temperature (set at 37.5°C), anaerobic conditions, and adhesion mechanisms. To simplify both understanding and computation, we applied a grid-based discretization, allowing for random diffusion of E. coli across different grid cells. The growth and death of E. coli are described by corresponding rate equations.

3. Model 3：Drug Release and Action:To simulate the process of AvCystatin acting on the human body after its production, we first established four phases—absorption, distribution, metabolism, and excretion—to describe the drug's release and absorption. The AvCystatin secreted by engineered bacteria is released in the intestines, with the release rate controlled by the concentration of thiosulfate and absorbed into the body following a first-order kinetic mechanism. Once inside the body, AvCystatin is distributed according to a two-compartment model, transferring between the blood (central compartment) and tissues (peripheral compartment). Metabolism mainly occurs in the liver, following Michaelis-Menten kinetics, regulated by the maximum metabolic rate and the enzyme saturation constant. When the drug concentration is high, the metabolic rate tends to become saturated. The excretion process occurs through the kidneys and other pathways, following first-order kinetics. From the literature, it is known that AvCystatin primarily acts on macrophages after entering the body, inducing the secretion of IL-10 and IL-12/23p40. Using IL-10 as an example, literature indicates that its production is regulated through a negative feedback mechanism in the ERK signaling pathway, which in turn regulates IL-10 concentration.

Oral administration is one of the safest and most cost-effective drug delivery methods, and it is also the most widely used route in humans. However, due to its inherent complexity, it is considered one of the most challenging drug delivery methods in terms of drug behavior within the body. This complexity primarily arises from the fact that drug efficacy depends not only on the drug itself and the patient’s absorption capacity but is also influenced by the contents of the gastrointestinal tract and its motility patterns. Additionally, the delivery method we designed is not a traditional pharmaceutical formulation but rather involves engineered bacteria. Therefore, understanding the movement and retention of these engineered bacteria in the gastrointestinal tract is crucial, as these factors directly affect the bacteria's ability to secrete therapeutic agents effectively, which is essential for disease treatment. Furthermore, due to the specific characteristics of inflammatory bowel disease (IBD), which can impair gastrointestinal motility, drug absorption is indirectly impacted. Thus, there is a strong justification for modeling this process.

On the one hand, it can improve the awareness of care workers about inflammatory bowel disease, so that they can better take care of elderly patients in their daily care work. On the other hand, promoting inflammatory bowel disease knowledge to the elderly can help improve their understanding and prevention awareness of the disease. Through these activities, we can raise social awareness of IBD, help patients achieve chronic disease self-management, and improve their quality of life.

First, a 3D model of the human stomach needs to be constructed, with the cardia and fundus serving as the model's entry point, and the pylorus along with a small segment of the post-pyloric duodenum serving as the exit (Figure 2). In this model, the stomach is primarily divided into four regions: the Fundus, Corpus, Antrum, and Pylorus (Figure 3).

The motion of the gastric wall is often defined as an antral contraction wave (ACW)[1],[2]. The starting point of ACW is considered to be a point \( s_{1} \) located at the proximal end of the antrum, near the corpus. This point is regarded as the pacemaker of ACW. These peristaltic waves propagate along the gastric wall toward the pylorus and terminate at \( s_{4} \), with the line connecting \( s_{1} \) and \( s_{4} \) defined as the ACW axis, also known as the central line of the antrum. The traveling cosine wave model of the gastric wall contraction strain \( \lambda_s \) is represented by the following equation[2]:

Figure 2. The 3D model of stomach

Figure 3. The flat model of stomach

In the equation above, \( s \) represents the distance along the central line, and \( \delta h(s) \) is the amplitude modulation function, used to vary the waveform intensity along the central line. \( s_{0} \) denotes the position within the half-width range of the wave. \( V_p \) is the pulse propagation speed, and \( T_p \) is the pulse interval. \( F(s, s_{0}) \) is a high-order filtering function that restricts the deformation within the half-width around the wave center, and \( W_p \) represents the width of the waveform.

In the above equation, the waveform of the antral contraction wave (ACW) is primarily controlled by \( F(s, s_{0}) \) and \( W_p \). However, considering that tracking wave propagation in physical space is critical in practical simulations and computational complexity should be minimized, the complexity of the \( \lambda_s \) formula can be significantly reduced by adjusting \( s_{0} \), while maintaining sufficient accuracy to describe the behavior of ACW. Thus, we have the simplified form:

Where \( \vec{s} \) is the direction vector of ACW along the central line after reaching a certain point, and \( \vec{s'} \) is the direction vector to another point within the half-width range of the wave center after the ACW has reached that point, as shown in Figure 4.

Figure 4. The flat model of stomach(Improved)

This simplification transforms a complex wave equation into a more physically meaningful vector representation, eliminating several physical quantities such as \( V_p \), \( T_p \), and \( W_p \). Instead, it restricts the range of ACWs using only vector forms. Additionally, since the original model is based on a cosine wave function, the normal vector of each point within the wave width (vibration direction) must be known, as shown in Figure 5. However, calculating these normal vectors for every point involves a substantial computational cost.

In the improved model, the vibration direction of each point within the wave width can be simplified by pointing towards the wave peak corresponding to the midpoint of the wave, assuming a uniform vibration amplitude (Figure 6).

As illustrated in the figure, when the time step is sufficiently small and the wave width \( W_p \) is narrow enough, the magnitude and direction of vibration between the two models show no significant difference, while the computational load is greatly reduced. This enables the model to achieve sufficient accuracy with less computational cost.

Figure 5. The flat model of stomach(Analysis)

Figure 6. The flat model of stomach(Improved and Analysis)

The amplitude modulation function and filtering function are given by the following equations:

The above equations indicate that when the ACW wave propagates, the gastric wall outside the half-width of the wave is not affected by it.

The amplitude modulation function \( \delta h(s) \) is defined as:

The points $s_1 - s_4$ denote different divisions along the central line (Figure 3). Specifically, \( s_1 \) is the pacemaker point of the ACW in the corpus, \( s_2 \) is the junction between the corpus and the antrum, \( s_3 \) is the boundary between the antrum and the pyloric canal, and \( s_4 \) is located at the pylorus.

To fully characterize and quantify the post-medication gastric motility patterns for refining the model, it is also necessary to determine the propagation speed \( V_{x} \) and period \( T_{p} \) of the ACW. The stomach model in Figure 1 is immersed in a Cartesian volume of 25 \times 9 \times 16 \, $cm^{3}$ along the x-, y-, and z-directions, as shown in Figure 5. The stomach model extends 15 cm in the x-direction from \( s_1 \) to \( s_4 \). The ACW initiates at \( s_1 \) every 20 seconds, and the propagation speed of the ACW along the x-direction is constant at \( V_{x} = 2.5 \, mm/s \). Thus, the lifetime of an ACW is approximately one minute. Given that the ACW period is \( T_{p} = 20 \, s \), this implies that up to three antral contraction waves may coexist on the stomach model at the same time, as depicted in Figure 6.

Figure 7. The diagram of ACWs

However, due to the curvature of the stomach, it is evident that the curvature of the upper and lower sides of the 2D stomach model is not uniform, and the curvature of the wave segments divided by different boundary points also varies. This variation in curvature results in non-uniform and non-constant propagation speeds of ACW across different segments and sides. Consequently, waves emitted at different times may overlap within the same segment.(Figure 7) Therefore, it is necessary to modify the traveling cosine wave equation into a summation form. Let the wave count be denoted as \( n_p \), then the equation becomes:

Where `\( \vec{s_n} \)` is the direction vector of the \( n \)-th ACW along the central line after reaching a certain point, and `\( \vec{s_n^{'}} \)` is the direction vector to another point within the half-width range of the wave center for the \( n \)-th ACW after reaching that point.

Due to the secretion of gastric acid, the pH value in the stomach ranges from approximately 1.5 to 3.5. Combined with the intense motility of the gastric wall, which leads to higher fluid pressure within the stomach, this environment is extremely unfavorable for the survival and proliferation of the engineered bacterial carrier (EcN 1917) used in this model. If the engineered bacteria are formulated as conventional capsules or disintegrating tablets that dissolve and release the bacteria in the stomach, they would struggle to survive under acidic and high-pressure conditions (as shown in Figure 8), significantly diminishing the therapeutic efficacy. Therefore, using non-disintegrating tablets or enteric-coated capsules as carriers is a more suitable option.

Figure 8. The reproductive curve of EcN in the stomach

Additionally, the Avcystatin protein produced by the engineered bacteria is primarily intended for anti-inflammatory effects in the intestines. However, this protein is highly sensitive to moisture, light, and oxygen. To prevent these factors from reducing drug efficacy and considering that the optimal single-dose amount of the engineered bacteria cannot yet be determined, the model employs No. 1, No. 2, No. 3, and No. 4 capsules as carriers for the engineered bacteria. All four capsule types are modeled as cylindrical shapes with spherical caps at both ends (as shown in Figure 9). This design effectively protects the engineered bacteria and the drug, ensuring their release in a suitable environment to achieve optimal therapeutic effects.

Figure 9. Capsule Size

In practice, patients usually take capsules with water, and this study assumes that the capsules are taken on an empty stomach. Therefore, the gastric contents at this time consist only of a mixture of fluids and gastric acid (assumed to be entirely in the liquid phase), whose properties are similar to those of water. The governing equation for the flow of an incompressible Newtonian fluid is the Navier-Stokes equation[3]:

The first equation indicates that the fluid is incompressible, meaning its density remains constant. The second equation is the Navier-Stokes equation, where \( \rho \) represents the fluid density, \( \vec{u} \) is the fluid velocity vector, \( \frac{\partial \vec{u}}{\partial t} \) is the time derivative of velocity, \( \vec{u} \cdot \nabla \vec{u} \) is the convective term, \( \nabla p \) is the pressure gradient force, \( \mu \) is the dynamic viscosity, \( \nabla^2 \vec{u} \) is the Laplacian of velocity, representing viscous diffusion, and \( \vec{g} \) is the gravitational acceleration. Considering that the patient is assumed to be standing while taking the capsule, the gravitational acceleration is oriented vertically downward.

The initial boundary condition of the gastric wall is set as a no-slip physical boundary. That is, if the gastric wall has not started moving, the velocity of the gastric fluid is zero. Therefore, a zero-gradient boundary condition for pressure is applied on the gastric wall, i.e., \( \vec{\nabla} p \cdot \hat{n} = 0 \). The gastric wall is also impermeable to the fluid. As for the boundary conditions at the fundus (inlet) and the pylorus (outlet), which are the two openings of the gastric container, it is assumed that the pylorus remains open, and the velocity across the pyloric cross-section is uniform. Thus, a zero-gradient pressure boundary condition is set at the pylorus. However, since the fundus and pylorus alternate between open and closed states in reality, keeping the pylorus open continuously would lead to a rapid reduction of gastric contents. Therefore, periodic boundary conditions are applied at both the fundus and pylorus, ensuring that the boundary conditions at the fundus match those at the pylorus. This approach ensures that the volume of gastric contents remains constant over a short period.

Since there are two-phase flows (capsules and fluid) within the stomach, a fluid-structure interaction (FSI) model is considered for solution. To describe capsule motion, a six-degree-of-freedom (6-DOF) model is selected[2]. This model characterizes the capsule's motion using six different fundamental modes of movement, as shown in Figure 10:

Figure 10. 6DOF

The equations governing the interaction between the capsule and the fluid are as follows:

Where \( m \) is the mass of the capsule, \( I \) is the moment of inertia of the capsule, \( \vec{v} \) and \( \vec{\omega} \) are the translational and angular velocities of the capsule, respectively. \( \vec{F}_f \) and \( \vec{M}_f \) represent the forces and moments induced by shear forces and pressure from the surrounding fluid, while \( \vec{F}_c \) and \( \vec{M}_c \) are the forces and moments generated by the contact between the capsule and the gastric wall. \( \vec{g} \) denotes gravitational acceleration.

The formulas for calculating the fluid-induced forces and moments are:

where \( S_p \) is the surface of the capsule, \( p \) is the fluid pressure, \( \hat{n} \) is the surface normal vector, \( \vec{\tau} \) is the viscous shear stress of the fluid, and \( \vec{r} \) is the position vector from the capsule's center of mass to a certain point on its surface.

Considering that the engineered bacteria are not released in the stomach, the capsule does not need to remain in the stomach for an extended period. Previous studieshave shown that the higher the specific gravity of the capsule (\( SG = \frac{\rho_{capsule}}{\rho_{water}} \))[2], the less it is affected by the retrograde jet flow in the stomach. This means that the capsule tends to stabilize near the pylorus, reducing its residence time in the stomach. Therefore, in this study, the capsule's specific gravity is set to SG = 1.3.

Since the engineered bacteria are released in the intestines rather than in the stomach, it is necessary to construct an intestinal model to investigate the subsequent movement of the capsule in the intestines and the diffusion behavior of the bacteria upon release.

Considering that the total volume and scale of the intestines are significantly larger than that of the stomach, it is impractical to simulate the entire intestinal tract for the patient. Therefore, we first segment the intestines for analysis. Once the capsule enters the small intestine, where the pH value shifts to neutral or even alkaline, the enteric capsule begins to dissolve and release the engineered bacteria. Since the dissolution time of the capsule is not prolonged, this study selects a segment of the transverse small intestine as the subject of analysis and assumes that the release of the engineered bacteria occurs within this segment. After the engineered bacteria are fully released, they proceed through the remaining part of the small intestine.

In this model, we abstract the small intestine as a transverse cylinder, with a length of approximately 15 cm, an inner diameter of 2.5 cm (i.e., a diameter of approximately 2.5 cm), and a wall thickness of 3 mm. Additionally, certain uneven surfaces are added to the inner wall to simulate folds and villi, as illustrated in the following schematic diagram:

The movement of the intestinal wall is similar to that of the gastric wall, propagating in the form of waves. However, since the small intestine is abstracted as a segment of a cylinder, with the centerline represented as a straight line, the wave form is relatively simpler. In this study, we adopt a sinusoidal wave pattern , which can be expressed as follows[4]:

Figure 11. The model of intestines

Among them, \( w \) represents the boundary coordinates of the wall, the superscripts \( u \)and \( d \) denote the upper and lower walls, respectively, while the subscripts \( u \) and\( l \) indicate the proximal and distal regions, respectively. Additionally, \( A \) represents the amplitude, and \( t \) denotes time. The period of peristaltic motion, \( T \), is set to 7.0 seconds, and the wavelength is set to 60 mm. The schematic of the peristalsis is illustrated as follows:

Figure 12. waves of intestines

For the contents of the small intestine, we treat them uniformly as a Newtonian fluid, so the intestinal contents also satisfy the Navier-Stokes equation:

The proximal opening is defined as the inlet, where the inflow velocity is set to 1 cm per minute, and the pressure boundary condition is set as a zero-gradient boundary condition. The distal outlet is defined as the outlet, similar to the gastric model. Both the inlet and outlet are configured as periodic boundaries to maintain a constant content volume within the model. The intestinal wall is set with a no-slip boundary condition.

Under these conditions, for fluid-structure interaction simulation, assuming that the capsule begins to dissolve and completely dissolves in this segment of the small intestine, according to the Nernst-Brunner equation, we have:

In this context,\( m_d \) is the mass of the dissolved substance, \( A_p \)is the surface area of the drug,\( D \) is the diffusion coefficient,\( \delta_d \) is the thickness of the apparent diffusion boundary layer, \( C_s \) is the solubility, \( C_\infty \) is the volume concentration in the medium. The mass transfer coefficient is given by\( k_m = D/\delta_d \).The volume concentration can be estimated by the equation\( C_\infty = m_d / V \) ,where \( V \) is the total volume of the medium.

Since EcN consists of cells, its diffusion coefficient in intestinal fluid is very small, typically around $10^{-12}$. Furthermore, engineered bacteria do not have solubility, so we set the colony-forming units (CFU) solubility \( CFU_s \) in the small intestine model (with a length of 15 cm, an inner diameter of 2.5 cm, and a volume of 73.59 mL) to replace the solubility in the equation when all engineered bacteria are fully dissolved. Correspondingly, the volumetric concentration of the medium \( C_\infty \)and the left-hand side term \( m_d \)in the equation are replaced by \( CFU_\infty \) in the medium and dissolved colony-forming units \( CFU_d \)respectively. Since the capsule type taken varies,\( CFU_s \)will also differ accordingly. Assuming the bacterial powder density is $10^{7}$ CFU/mg, with a bulk density of 0.6–0.8 g/mL, a plot of the mass dissolution of different capsule types over time can be generated as follows:

Figure 13. CFU Dissolution of Different Capsules Over Time

From the figure, it can be observed that the larger the surface area and volume of the capsule, the faster the dissolution rate. After the capsule dissolves, the engineered bacteria will flow along with the intestinal fluid, and when they reach the diseased segment of the intestine, they will begin to secrete the drug. Additionally, due to the folds and villi of the small intestine and large intestine, a considerable portion of the EcN bacteria will not flow forward with the intestinal fluid until being excreted by the patient, but will remain in the intestines to continue secreting the drug. Therefore, at this point, it is necessary to establish a reproduction model for the engineered bacteria in order to quantitatively analyze the survival conditions of the bacteria and the efficacy of the drug.

Due to the fact that the colon accounts for a large proportion of the large intestine and that engineered bacteria primarily proliferate in the colon, our focus is mainly on the colon section. For the colon, the modeled cylindrical length is approximately 5 cm, with an inner diameter of 6 cm (i.e., the diameter is approximately 6 cm), and a wall thickness ranging from 2 to 4 mm. Regarding the contents of the colon, we only consider their effects in terms of erosion and nutrient supply (i.e., fluid flows from one end of the cylinder to the other, eroding and carrying away part of the engineered bacteria while replenishing nutrients). In the intestines of most patients with inflammatory bowel disease (IBD), the pH value of the gut remains largely unchanged [10]. We set the pH value of the intestines to 7, the temperature to 37.5°C, and assume an anaerobic environment. Nutrient levels decrease over time and are replenished when the contents of the small intestine pass through. Based on the above information, we proceed with constructing a model of the colon engineered bacteria in vivo. We simulate the scenario where the intestinal contents are present and flow slowly through the large intestine.

First, we consider the fluid dynamics equations governing the flow of chyme. Based on the work of Ibitoye S. E. et al. [8], we obtain the basic physical properties of intestinal chyme. We assume that the viscosity of chyme within the human body is similar to the experimental data measured in their study, and that the flow is steady, without the influence of gravity. Due to the increased intestinal peristalsis in IBD patients, the absorption of nutrients is weakened. Therefore, we consider that the water content of the chyme is relatively higher, allowing us to simplify it as a Newtonian fluid. The physical properties of the chyme flowing into various parts of the intestine remain essentially unchanged, and the inflow velocity is set at 1 cm per minute. In a three-dimensional cylindrical coordinate system, the flow characteristics of a Newtonian fluid are described by the Navier-Stokes equations. The Navier-Stokes equations are expressed as:

Figure 14. Pressure field

Figure 15. Velocity field

Where: $\mathbf{U}=(U_r,U_\theta,U_z)$ represents the velocity field components (the velocity components along the $r$, $\theta$, and $z$ axes, respectively); $\rho$ is the fluid density; $\mu$ is the viscosity.

Where, $\rho$ is the fluid density, $v$ is the flow velocity, $d$ is the pipe diameter, and $\mu$ is the fluid viscosity.

After performing the calculations, it satisfies the condition for laminar flow, i.e., the Reynolds number $Re$ is less than the critical value (2300), ensuring that the flow is steady and laminar rather than turbulent. We apply a no-slip boundary condition at the fluid boundary: $\mathbf{U} | _{\text{boundary}} = 0$, with a constant average inflow velocity, and zero pressure at the outflow boundary. Additionally, the fluid is incompressible, homogeneous, and has constant density and viscosity. Therefore, we can regard the flow as Poiseuille flow. After simulation, the resulting pressure field and velocity field are as follows:

For the proliferation of engineered bacteria, we refer to the experimental conclusions of Lucija K. et al. [9], assuming that their proliferation in LB medium is similar to that in the intestinal environment. The concentration distribution of the engineered bacteria is influenced by the flow of the non-Newtonian fluid, diffusion, external inflow, and the concentration of nutrients. By incorporating the convection-diffusion equation, we obtain the evolution equation as:

Figure 16. Engineered bacteria random movementn

Where: $C(r,\theta,z,t)$ is the concentration of the engineered bacteria at time $t$ and position $(r,\theta,z)$; bf{U}=(U_r,U_\theta,U_z)$ is the velocity field obtained from the solution of the Navier-Stokes equations; $D$ is the diffusion coefficient of the engineered bacteria; $\alpha$ is the interaction parameter between the engineered bacteria and nutrients; $\beta$ is the interaction parameter between the engineered bacteria (such as competition effects).

According to the work of Wang Z. et al. [7], for the diffusion coefficient, we have:

Where: $V$ is the velocity, with $E[V] = 24.14\mu m/s$, and assuming its standard deviation is small, it can be approximated as a constant. $\tau$ is the run duration, which follows an exponential distribution with a mean value of $E[\tau] = 0.84s$. $\theta$ and $\phi$ are the polar angle and azimuthal angle, both following uniform distributions:

Finally, we obtain the diffusion coefficient as $6.3430 \times 10^{-8}\text{m/s}$.

For the nutrients in chyme, we only consider substances that promote the growth of Escherichia coli, and based on the work of Lucija K. et al. [9], we assume that the nutrient concentration is consistent with the LB medium containing $N_M = 2.1 \times 10^9\text{CFU/ml}$. We assume that the engineered bacteria exhibit similar proliferation behavior in both liquid and non-Newtonian fluids. For the nutrient concentration, we assume it is proportional to the concentration of the components peptone and sodium chloride (NaCl), and that the engineered bacteria consume the nutrients in proportion to the initial composition of the nutrients. For the two components, we use the Wilke-Chang equation to calculate the diffusion coefficient for NaCl. For peptone, based on available data, the diffusion coefficient in water is $10 \times 10^{-6} \text{cm}^2/\text{s}$. We assume that this diffusion coefficient is inversely proportional to the viscosity of the solvent and directly proportional to temperature, thus estimating its diffusion rate in chyme. By applying a weighted average based on the initial nutrient composition for these two components, we obtain the average diffusion coefficient. Since the nutrient concentration is proportional to the concentration of its constituent substances, the average diffusion coefficient can be approximately considered as the diffusion coefficient of the nutrients. The distribution of the nutrient concentration is affected by fluid flow, diffusion, external inflow, and interactions with the engineered bacteria, and will be consumed over time. Its evolution equation in cylindrical coordinates is given by:

Where: $N(r,\theta,z,t)$ is the concentration of nutrients at time $t$ and position $(r,\theta,z)$; $\mathbf{U}=(U_r,U_\theta,U_z)$ is the velocity field obtained from the solution of the Navier-Stokes equations; $D_N$ is the diffusion coefficient of the nutrients.

For the concentration of engineered bacteria, we set the initial condition as $C(r,\theta,z,t)$. At time $t=0$, the engineered bacteria are only present in a certain region of the fluid, and no additional bacteria are introduced as the fluid continues to flow:

Assuming that a portion of the engineered bacteria adheres to the wall surface , we employ the adhesion model:

Where $a$ is the adhesion coefficient. For the nutrient concentration $N(r,\theta,z,t)$, it is assumed to be uniformly distributed in the inflowing fluid at the initial time, and since chyme passes through the colon over a period of 24 to 36 hours, we continuously replenish the nutrients during subsequent simulations. For the adhesion coefficient, let $C_1$ represent the concentration at $r = R$. Using the previously defined adhesion rate formula, we obtain:

Figure 17. The curve of total bacterial concentration variation

For a single colonic epithelial cell, we assume its surface area to be $9.42 \times 10^{-10}m^2$. Assuming that the attached bacteria are distributed within a thin film layer on the surface of the epithelial cells, with a thickness of 1 micron, the surface volume concentration can be approximated as $C_1 \approx 1.70 \times 10^{16}\text{CFU/m}^2$. Assuming that the bacterial concentration variation primarily occurs in a small region near the wall, with the region thickness $\Delta n$ set to 0.01 cm (1 millimeter), for the concentration gradient, we have:

Where $C_\text{2}$ is the concentration of the suspension far from the epithelial cells. Based on the experimental data from the paper, we take $C_1=5\times 10^{15}\text{CFU/cm}^3$. Substituting into the formula, we obtain the adhesion coefficient as $a \approx 4.48 \times 10^{-2}\mathrm{mm}^{-1}$.

We simulated four different drug administration methods as described above, using the dissolution rate of the capsule at 60 minutes as the initial concentration of engineered bacteria. The behavior of the engineered bacteria in the gut was simulated over a five-hour period. We spatially discretized the model and set the time step to 0.01 minutes. For the engineered bacteria released from the four different types of capsules, we first simulated the total amount of engineered bacteria in the colon and the consumption of nutrients, as shown in the figure below:

From the figure, we can observe that the quantity of each type of engineered bacteria reaches a peak almost simultaneously during the first hour (the first 6000 time steps). Subsequently, due to nutrient consumption, the flushing effect of chyme, and competition among the engineered bacteria, their numbers rapidly decline and level off after about one and a half hours. Furthermore, despite the differences in initial quantities, the decline trends are nearly identical for all four cases. In other words, regardless of the initial quantity, the population always stabilizes after a certain period of time. The nutrient depletion trend also follows a nearly identical pattern. Therefore, we will focus on the bacterial growth during the first hour in the following analysis.

Figure 18. Curve of Changes in the Total Concentration of Nutrients

Figure 19. Curve of Changes in the Total Concentration of Bacteria (Half an Hour)

From the figure, we can observe that for Capsule Size01, it reaches the maximum quantity the fastest, and the maximum quantity it achieves is also the highest among the four capsule types. There is a positive correlation between the initial quantity of engineered bacteria released by the four capsules and the maximum quantity produced through reproduction in the gut. This information suggests that Capsule Size01 should generally be preferred. Although the decline in engineered bacteria in the gut slows after two hours, due to the depletion of nutrients, the flushing effect of intestinal contents, and the diffusion of bacteria, a replenishment is required after several hours, which aligns with real-world conditions. Next, we simulate the reproduction of engineered bacteria in the gut. In the following four figures, due to the cylindrical symmetry, we selected a cross-section at a specific angle for the first and third figures. In the second and fourth figures, we visualized the entire cylindrical region, using concentration to represent the color and transparency of each point (the lower the concentration, the more transparent the color).

Figure 20. Reproduction Diagram of Engineered Bacteria in the Gut (Angle 180)

Figure 21. 3D Reproduction Diagram of Engineered Bacteria in the Gut

From the figure, we can observe that at the beginning, the engineered bacteria are rapidly transported forward due to convection and diffusion, and accumulate on the intestinal wall because of the adhesion effect. During the first one to two hours, the flow velocity at the edges is relatively slow, and the accumulation of substances results in a higher concentration in these areas. In contrast, the concentration of engineered bacteria in the central region is lower due to the faster flow velocity and the influence of surrounding substances. After three hours, the concentration at the edges remains comparatively higher but gradually decreases over time, while an anomalous concentration appears in the center. We can better analyze this phenomenon through the following figure.

In this figure, observing the second image, we can see that the engineered bacteria, influenced by the flow of chyme, move rapidly at the beginning, causing the region with the highest concentration of bacteria on the intestinal wall to gradually shift forward. At the middle of the entrance, after two hours, there remains a certain amount of engineered bacteria due to the continuous supply of nutrients at the entrance and the weaker adhesion effect of the intestinal wall. After the initial process, the bacteria at this location reach a stable state, although their quantity is relatively low after the changes in the first two hours. At this point, we hypothesize that they have entered a "pseudo-lag phase." In the final hour, there is a slight increase in bacterial concentration, suggesting that the engineered bacteria at the middle of the entrance may enter a "pseudo-logarithmic phase" of growth after a certain period of reproduction.

Figure 22. Nutrient Variation Diagram (Angle 180)

From the figure, we can observe a similar trend with nutrients. Due to bacterial consumption and the effects of convection, nutrients also accumulate on the intestinal wall and gradually decrease as they are consumed over time. Since nutrients are continuously replenished at the entrance, the decrease in nutrient concentration at the entrance is less significant. However, because a portion of the nutrients has already been consumed by the engineered bacteria, the nutrients in the middle section along the z-axis receive less replenishment, making them more affected by convection and diffusion. The reduction in nutrient concentration on the intestinal wall is relatively small, but since engineered bacteria accumulate on the intestinal wall, the consumption of nutrients is higher in that region. This indicates that nutrients also tend to accumulate on the intestinal wall.

Based on the above images, we can observe that after two hours, the nutrient concentration in most areas has significantly decreased, and by the third hour, it is almost entirely depleted. In summary, we conclude that it is preferable to select engineered bacteria with better adhesion to the intestinal wall to ensure that they can accumulate in large quantities on the colonic wall, thereby producing more AvCystatin. However, how do these substances interact with the human body to alleviate symptoms? This requires us to model their absorption and mechanism of action.

Figure 23. 3D Nutrient Variation Diagram

Based on the literature review and experimental results, we know that the secretion rate of AvCystatin by the engineered bacteria can be expressed by the following equation:

Where: $C_{\mathrm{AvCyst}}$ represents the concentration of AvCystatin, $S_{\mathrm{thio}}$ represents the concentration of thiosulfate, $r_{\mathrm{max}}$ is the maximum production rate of AvCystatin, $K_S$ is the half-saturation constant for thiosulfate, representing the concentration at which thiosulfate affects AvCystatin production.

For AvCystatin secreted into the intestine, we consider its absorption, distribution, metabolism, and excretion (ADME) processes. The absorption phase represents the process by which the protein enters systemic circulation from the intestine. We assume the absorption follows first-order kinetics. Let $A(t)$ represents the amount of AvCystatin in the intestine at time $t$ (in mg), $k_a$ is the absorption rate constant (in h$^{-1}$), $k_{deg}$ is the degradation rate constant in the intestine (in h$^{-1}$), $F$ is the bioavailability (dimensionless). The equation is:

For a first-order linear differential equation, the general solution can be expressed as:

Where $A(0)$ is the initial amount of AvCystatin in the intestine at time $t=0$. Substituting into the formula for the drug concentration after absorption, we have:

Substituting$A(t^\prime)=A(0)\cdot e^{-(k_a+k_{deg})t^{\prime}}$into the equation, we get:

The solution is:

The distribution of absorbed proteins in the body can be described using a two-compartment model. We define the central compartment as the blood (plasma) and the peripheral compartment as the tissues. Let $V_c$ represent the volume of the central compartment (units: L), $V_p$ represent the volume of the peripheral compartment (units: L), $C_c(t)$ represent the concentration in the central compartment (units: mg/L), $C_p(t)$ represent the concentration in the peripheral compartment (units: mg/L), $k_{cp}$ represent the transfer rate constant from the central compartment to the peripheral compartment (units: h$^{-1}$), and $k_{pc}$ represent the transfer rate constant from the peripheral compartment to the central compartment (units: h$^{-1}$). Thus, the equations are:

Define the state vector $\mathbf{C}(t)=\begin{pmatrix}C_c(t)\\C_p(t)\end{pmatrix}$.Thus, the system of equations can be expressed as:

Where:

For solving the homogeneous system of equations:

We solve for the eigenvalues and eigenvectors of matrix $\mathbf{A}$

The characteristic equation is：

The eigenvalues are:

We obtain:

The concentration in the peripheral compartment is:

For metabolism, the metabolic process primarily occurs in the liver, and we assume it follows first-order kinetics. Considering the enzyme saturation effect, Michaelis-Menten kinetics is applied. Thus, we define: $k_{met}$ is the metabolic rate constant (unit: h$^{-1}$), $V_{max}$ represents the maximum metabolic rate (unit: mg/h), $K_m$ is the Michaelis constant (unit: mg/L). The equation is:

The excretion of proteins involves trace amounts being excreted through pathways such as the kidneys, which we assume follows first-order kinetics. Let: $\bullet$ $k_{ex}$ be the excretion rate constant (unit: h$^{-1}$). The excretion kinetics equation is:

We obtain:

The specific process is outlined as follows:

Figure 24. Pharmacokinetic Flowchart of AvCystatin

By reviewing the works of Figueiredo A. S. et al.[5] and Christian K. et al.[6], we learn that AvCystatin, after being absorbed into the human body, typically acts on macrophages, leading to the production of IL-10 and IL-12/23p40. Taking IL-10 as an example, we assume that AvCystatin binds to receptors on macrophages without being consumed, and according to the literature, IL-10 activates phosphatases that negatively regulate ERK signaling via feedback. Based on this, the following hypotheses are formulated: We assume that IL-10 regulates ERK dephosphorylation in the ERK signaling pathway through the activation of phosphatases. Specifically, extracellular IL-10 (denoted as IL-10e) binds to IL-10 receptors on the surface of macrophages, thereby activating a certain phosphatase. This phosphatase promotes ERK dephosphorylation and inhibits its activity. This implies that ERK phosphorylation is initially triggered by AvCystatin , but as IL-10 is secreted and accumulates, ERK is rapidly dephosphorylated, forming a negative feedback loop. We hypothesize that the activation patterns of the ERK and p38 signaling pathways are distinct. ERK activation is transient, as it is regulated by the negative feedback of IL-10, whereas p38 activation is sustained and is not directly influenced by IL-10 feedback. Furthermore, we hypothesize that IL-10 is induced not only by external stimuli but also regulated by an autocrine mechanism that further influences its own production. In macrophages, IL-10 activates a feedback mechanism through its receptor, affecting ERK activity and consequently regulating IL-10 expression levels. Based on the above content, and choosing Model 2 from Figueiredo A. S. et al.'s paper[5], we can derive the following system of differential equations:

The relationships and effects between the substances are described as follows:

Figure 25. AvCystatin signaling pathway diagram

Figure 26. Time-concentration curve

The first figure illustrates the interactions between multiple signaling pathways and substances within the model. In the upper left section, ERK and its phosphorylated state (erkp) form a feedback loop, maintaining ERK activation through a closed circuit. In the upper right section, p38 undergoes feedback regulation through its phosphorylated form (p38p), sustaining its activity. The formation of IL-10 involves multiple steps, transitioning from its mRNA form (IL10m) to intracellular protein (IL10i), and finally to extracellular IL-10 (IL10e). These steps are influenced by different feedback pathways that also affect external variables (e.g., extVar). Ap is a regulatory factor affecting many reaction steps, and through interactions with other substances, it further modulates the system. The species X (X0, X1, X2) influence each other through a series of feedback mechanisms, amplifying or attenuating the signals. The second figure shows the time-course dynamics of several variables in the model. The most notable change occurs in IL10i, which peaks around time point 3 and then gradually declines. This reflects the rapid production of IL10i followed by negative feedback regulation, as shown in the first figure. X0 and A have relatively high and stable levels (approximately 12 units), indicating that they remain in an active state throughout the secretion process. The fluctuations in erk and IL10e are minor, and the decrease in erkp corresponds closely with the increase in erk. IL10m experiences a slight rise followed by a decline. Both p38 and p38p exhibit some early variation but remain relatively stable thereafter, consistent with the self-sustaining behavior described in the first figure. X1 and X2 show almost no change, indicating they are not dominant substances in this system. The increase in Ap mirrors the decrease in A, with both showing minimal variation. Together, the two figures illustrate the process by which AvCystatin promotes the production of IL-10 in macrophages. IL-10, as an anti-inflammatory cytokine, reduces intestinal inflammation and alleviates symptoms in IBD patients. Additionally, IL-10 promotes the production of NOS, leading to a reduction in NO levels, which in turn induces apoptosis in engineered bacteria[11]. This mechanism effectively controls the proliferation of engineered bacteria in the human intestine, ensuring the safety and effectiveness of the treatment. The model provides a comprehensive understanding of the interactions between AvCystatin and the human body, offering insights into the therapeutic effects of engineered bacteria on IBD.

[1] Kuhar S, Lee J H, Seo J H, et al. Effect of stomach motility on food hydrolysis and gastric emptying:
Insight from computational models[J]. Physics of fluids, 2022, 34(11).

[2] Seo J H, Mittal R. Computational modeling of drug dissolution in the human stomach[J].
Frontiers in Physiology, 2022, 12: 755997.

[3] Palmada N, Hosseini S, Avci R, et al. A systematic review of computational fluid dynamics models
in the stomach and small intestine[J]. Applied Sciences, 2023, 13(10): 6092.

[4] Yasumasa Ito, Koichiro Asato, Inhyeok Cho, Yasuhiko Sakai, Koji Iwano, Takahisa Tainaka, Hiroo Uchida,
Intestinal flow after anastomotic operations in neonates,Computers in Biology and Medicine,Volume 118,2020,103471,ISSN 0010-4825,

[5] Figueiredo A. S., Höfer T., Klotz C., et al. Modelling and simulating interleukin-10 production and regulation
by macrophages after stimulation with an immunomodulator of parasitic nematodes[J]. FEBS Journal, 2009, 276(13): 3454-3469.

[6] Christian K ,Thomas Z ,Sofia A F , et al. A helminth immunomodulator exploits host signaling events to regulate
cytokine production in macrophages.[J].PLoS pathogens,2011,7(1).

[7] Wang Z., Kim M., Rosen G. Validating models of bacterial chemotaxis by simulating the random motility coefficient[C].
Proceedings of the 2008 8th IEEE International Conference on BioInformatics and BioEngineering (BIBE), Athens, Greece, 2008: 1-4.

[8] Ibitoye S. E., Adegun I. K., Omoniyi P. O., Ogedengbe T. S., Alabi O. O. Numerical investigation of thermo-physical properties
of non-newtonian fluid in a modelled intestine[J]. Journal of Bioresources and Bioproducts, 2020, 5(3): 211-221.

[9] Lucija K ,Matilda Š ,Ana M , et al. A simple interaction-based E. coli growth model.[J].
Physical biology,2019,16(6):066005.

[10] Yang Ying, Chen Longdian, Chen Zairong, et al. Study on the Therapeutic Effects of Nitric Oxide Synthase Inhibitors
on Experimental Colitis in Rats [J]. Gastroenterology, 2007, (01): 27-30.

[11] Chen Yuzhou, Wang Dongkai. Physiological Characteristics of the Colon and Colon-Targeted Drug Delivery in Inflammatory Bowel Disease
[J]. Chinese Pharmaceutical Journal, 2006, (02): 86-90.