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.
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.
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.
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.
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.
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:
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]:
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:
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:
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:
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:
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:
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.
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).
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.
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.
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:
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:
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.