As the core site for digestion and absorption of nutrients in the human digestive system, the peristaltic and segmental motions of the small intestine play a crucial role in promoting the mixing of intestinal contents, accelerating substance transport, and improving absorption efficiency. Accurate simulation of the dynamic environment of the small intestine and the product transport process provides guidance for predicting product usage and action sites. This study develops a multi-physics coupled numerical model of the human small intestine. The specific objectives are:
(1) Construct a geometric model of the small intestine (including the duodenum, jejunum, and pylorus) based on physiological parameters to simulate the synergistic effect of peristaltic waves and segmental motions;
(2) Integrate fluid flow, dilute substance transport, and wall adsorption-desorption kinetics into the model to realize coupled simulation of multi-physical processes;
(3) Verify the reliability of the model through theoretical analysis and comparison with indirect experimental data, providing a basis for subsequent studies on bacterial colony adhesion and microplastic degradation.

(1) Axisymmetric Simplification: The digestive tract (stomach, pylorus, duodenum, jejunum) is constructed as a 2D axisymmetric model, ignoring circumferential (φ-direction) changes and only considering radial (r) and axial (z) motions and transport.

(2) Peristaltic Wave Simplification: The radial contraction of peristaltic waves in the jejunum is described using a Gaussian distribution function. The gastric peristaltic waves are simulated by superimposing a sine-squared function with a smooth step function. The opening and closing of the pylorus and changes in the duodenal radius are achieved through piecewise functions to ensure smooth spatial and temporal transitions.

(3) Segmentation Wave Assumption: The segmentation of substances in the jejunum and duodenum is realized through wall displacement in the form of a sine wave, and the segmentation region is limited to a specific axial range (constrained by a smooth function).

(1) Incompressible Laminar Flow: The fluid density (ρ = 1040 kg/m³) and dynamic viscosity (μ = 0.0014 Pa·s) are constant^[3]. The flow satisfies the Navier-Stokes equations, and turbulence (turbulence model is turned off) and gravity effects are ignored.

(2) Dilute Substance Assumption: The initial concentration of the product bacterial solution is low (initial concentration C₀ = 1 mol/m³), which does not affect the fluid viscosity and density. Substance transport only considers convection and molecular diffusion (turbulent diffusion is ignored). Since Colony-Forming Units (CFU) are not the dependent variable unit used in COMSOL, mol is used instead of CFU in this study. This is only a difference in measurement units, and there is no numerical difference between [mol/m³] and [CFU/m³] in numerical calculations.

(3) Wall Adhesion Kinetics: The adsorption-desorption process of the bacterial solution on the intestinal wall conforms to the Langmuir adsorption model, and the desorption rate is enhanced by the wall shear stress (regulated by the shear sensitivity coefficient α).

(1) Moving Mesh Compatibility:The computational mesh required for solving all control equations of the fluid flow model is created using the built-in mesh generation code in COMSOL Multiphysics® software. The mesh deforms synchronously with the wall peristalsis, and the "deformed domain" physics field is used to update the mesh displacement. The mesh smoothing type is Yeoh, and the boundary layer mesh is refined to capture the wall shear effect. There are 16,972 mesh vertices, including 22,279 triangular elements, 4,578 quadrilateral elements, and 26,857 domain elements.

(2) Time Step and Convergence:An adaptive time step scheme is adopted to automatically adjust the solver time step size to maintain an absolute tolerance of ≤ 0.001. The total duration of the transient calculation is 600 s, and the output results are stored with a time step of 0.1 s. The PARDISO linear solver is used, and Anderson acceleration is applied for nonlinear iterations to stabilize convergence.

(3) Physics Field Coupling:Four physics fields are set up: laminar flow, boundary ordinary differential and algebraic equations (describing the adsorption-desorption process), dilute substance transport, and moving mesh. All four physics fields are involved in coupling during numerical calculations.

This numerical model is constructed using COMSOL Multiphysics® 6.3 software, covering four core physics fields: moving mesh (simulating intestinal wall motion), laminar flow (simulating fluid flow), dilute substance transport (simulating substance diffusion and convection), and boundary ordinary differential equations (simulating wall adsorption-desorption kinetics). The model adopts 2D axisymmetric geometric simplification (ignoring circumferential changes) to balance computational efficiency and simulation accuracy. The key parameters are detailed in the following tables.

Parameters for Geometry, Fluid Properties, Laminar Flow, Moving Mesh, and Dilute Substance Transport

Name	Expression	Unit	Value	Description
L	1.225	[m]	1.225 m	Tube length (changed to 1.80 [m] for the ileum)
rho	1040	[kg/m^3]	1040 kg/m³	Fluid density
mu	0.0014	[Pa*s]	0.0014 Pa·s	Dynamic viscosity
amp	0.73	-	0.73	Relative contraction amplitude
tc	7	[s]	7 s	Period
ta	1.75	[s]	1.75 s	Contraction gating start time
tb	5.25	[s]	5.25 s	Relaxation gating start time
thalf	1.75	[s]	1.75 s	Gating smoothing time
z0	L/2	[m]	0.6125	Peristalsis center (0.90 [m] for the ileum)
up	7.50E-04	[m/s]	7.5E-4 m/s	Peristaltic propulsion velocity
sigma2	2.00E-04	[m^2]	2E-4 m²	Gaussian variance (= 2 cm²)
R0	0.00865	[m]	0.00865 m	Intestinal radius
zc	0.05*L	-	0.06125 m	Axial position z of the bolus center
rc	0	[m]	0 m	Radial position r of the bolus center (m); set to 0 when placed near the axis
az	0.02	[m]	0.02 m	Semi-axis of the bolus in the z-direction (m)
ar	0.8*R0	-	0.00692 m	Semi-axis of the bolus in the r-direction (m)
sm	0.5*min(az,ar)	-	0.00346 m	Smoothing half-width to avoid oscillations caused by hard edges
lambda	0.07	[m]	0.07 m	Segmentation wave wavelength
xigema	0.25*R0	-	0.0021625 m	Segmentation wave amplitude
f	10	[1/min]	0.16667 1/s	Segmentation frequency
a	8.5*lambda-d	-	0.59 m	Segmentation wave smoothing start point
b	9*lambda+d	-	0.635 m	Segmentation wave smoothing end point
d	0.005	[m]	0.005 m	Segmentation smoothing transition half-width
dPdz	-0.04	-	-0.04	Target axial pressure drop [Pa·s]
target	5.64E-08	[m^3]	5.6415E-8 m³	Target cumulative volume per second
L1	0.245	[m]	0.245 m	Duodenal length
L2	0.159	[m]	0.159 m	Stomach length
R1	0.0386	[m]	0.0386 m	Stomach fundus width
L1.5	0.015	[m]	0.015 m	Pyloric region length
A	0.9	[cm]	0.009 m	Gastric peristaltic wave amplitude
v	2.34	[mm/s]	0.00234 m/s	Gastric peristaltic wave propagation velocity
L4	2.04	[cm]	0.0204 m	Gastric peristaltic wave width
delta	0.1*L4	-	0.00204 m	Spatial smoothing parameter, L is half the period length
delta_t	0.5	-	0.5	Temporal smoothing parameter
tanth	L2/(R1-R0)	-	5.3088	Tangent value
sinth	abs(tanth)/sqrt(1+tanth^2)	-	0.98272	Sine value
R2	1.05	[cm]	0.0105 m	Duodenal radius
t0	0	-	0	Time when the duodenum starts to change
tau_on	2	[s]	2 s	Temporal smoothing activation duration
zw	0.05*L1	-	0.01225 m	Spatial transition length at both ends
L3	0.002	[m]	0.002 m	Thin domain length

Adhesion Parameters

Name	Expression	Unit	Value	Description
D	1.00E-13	[m^2/s]	1E-13 m²/s	Bulk diffusion coefficient
alpha	1.00E-07	[1/Pa]	1E-7 1/Pa	Shear sensitivity coefficient
kon	1.00E-05	[m/s]	1E-5 m/s	Adsorption coefficient
koff0	1.00E-04	[1/s]	1E-4 1/s	Basal desorption rate
Gmax	1.00E-05	[mol/m^2]	1E-5 mol/m²	Maximum adsorption capacity
C0	1	[mol/m^3]	1 mol/m³	Initial concentration in the bolus

Geometric Model

For numerical purposes, the geometry of the small intestine is assumed to be a perfect cylinder with flexible walls, and the entire geometry is assumed to be filled with liquid intestinal contents. The flow in the cylindrical tube is assumed to be axisymmetric, and the gastric antrum is assumed to have a funnel-shaped feature. All geometric parameters are based on real physiological data. The geometric model includes four parts: the gastric antrum (only including the starting point of gastric peristalsis to the pylorus), the pylorus, the duodenum, and the jejunum (the ileum model is constructed separately for comparison). The "Form Union" operation in COMSOL is used to construct a unified 2D axisymmetric computational domain^[1][3][5].

(1) Gastric Antrum: The diameter of the upper surface is 1.73 cm, the diameter of the lower surface is 7.71 cm, and the length is 15.9 cm. Its volume in 3D is consistent with the typical volume of the gastric antrum.

(2) Pyloric Region: The diameter is 1.73 cm, and the length is 1.5 cm.

(3) Duodenum: The diameter in the relaxed state is 2.1 cm, and the length is 24.5 cm.

(4) Jejunum: The diameter is 1.73 cm, and the length is 120 cm^[1][3][5].

Fig 1. 2D and 3D Axisymmetric Model Geometry (1 - Stomach, 2 - Pylorus, 3 - Duodenum, 4 - Jejunum)

Model Assumptions

The intestinal contents are assumed to be incompressible Newtonian fluids. Since the expected Reynolds number (Re) is ~3, the velocity field of the intestinal contents in the numerical models of the jejunum and ileum generated by peristaltic waves is assumed to be laminar. The intestinal contents are assumed to be Newtonian fluids with a dynamic viscosity of 0.0014 Pa·s and a density of 1040 kg/m³^[3]. Due to the complex spatial position of the small intestine in the human body, the effect of gravity is relatively small compared to the peristaltic driving force, so the effect of gravity is ignored. The flow is axisymmetric, and the circumferential (φ) velocity component and concentration gradient are zero.

Governing Equations of the Laminar Flow Model

The continuity equation and Navier-Stokes equations are used to solve the laminar fluid flow in the entire domain. The opening side of the gastric antrum and the end of the jejunum are set as open boundaries. A moving boundary condition is applied to all boundaries representing the digestive tract wall except the axis of symmetry boundary.

(1) Navier-Stokes Equations

The Navier-Stokes equations solve the gradients of velocity with respect to pressure, viscous forces, and gravity.In Equation，t is time, p is pressure (Pa), μ is dynamic viscosity (Pa·s), ∇u is the velocity gradient tensor, (∇u)t is the transpose of ∇u (s⁻¹), and g is the gravity vector (m/s²). In this study, the effect of gravity on the flow is ignored. The intestinal contents are assumed to be Newtonian fluids with a dynamic viscosity (μ) of 0.0014 Pa·s and a density (ρ) of 1040 kg/m³.

(2) Open Boundary Condition

The open edges at both ends of the axisymmetric geometry are defined as open boundaries to describe the boundaries in contact with a large amount of fluid. The given open boundary condition is applied to these two edges, indicating that the sum of pressure and viscous forces is zero. The fluid can flow freely in and out of the domain without external stress. In the equation, represents the boundary normal pointing outward of the domain.

(3) Wall Condition

This is the most commonly used wall condition in fluid mechanics. It assumes that at the contact surface between the fluid and the solid wall, the relative velocity of the fluid is zero. That is, the fluid molecules close to the wall will "stick" to the wall without relative sliding with the wall.

Moving Mesh

The dynamic mesh module is used to simulate the deformation of the digestive system wall caused by peristalsis and segmentation. The motion modes simulated in the model include: gastric peristalsis, intermittent opening and closing of the pylorus and gastric emptying, duodenal segmentation, jejunal peristalsis, and jejunal segmentation.

Introduction to Gastric Peristalsis

In general, the motion pattern of the gastric wall can be characterized by two types of muscle contractions, as described below. The first type of motor activity originates from the upper part of the stomach and propagates. It is characterized by slow and weak wave activity, which only produces slight indentations in the gastric wall. The second type of motor activity is characterized by a series of regular antral contraction waves (ACWs) that originate from the middle of the stomach and propagate circumferentially toward the pylorus. These contractions are born as shallow indentations but deepen as they move downward (nearly blocking the gastric lumen when they approach the pylorus). Their frequency is approximately 3 cycles per minute, and they are controlled by electrical stimulation generated by the gastric pacemaker, which is located in the middle body region of the greater curvature of the stomach.

In the current modeling stage, we have been able to identify most of the parameters based on literature data. The peristaltic pattern of the gastric wall is characterized as follows: ACWs are initiated every 20 seconds at a location 15.9 cm from the pylorus, with a lifespan of 61 seconds^[1].

Fig 3. Schematic Diagram of Gastric Peristaltic Waves

The relative occlusion degree is defined as RO=A/D, where D is the distance from the gastric wall to the central axis, and A is the amplitude of the gastric peristaltic wave. The variation of RO is specified as follows: from 0% to 40% linearly from 0 s to 19.5 s, remains constant at 40% from 19.5 s to 35.5 s, increases linearly again from 35.5 s to 60 s, reaches a maximum of 80% when advancing toward the pylorus, and finally decreases to 0 gradually from 60 s to 61 s.

Governing Equations of Gastric Peristalsis
(1) The generation of gastric peristaltic waves is controlled by setting the displacement of the gastric wall. The normal displacement of the gastric wall is:

Fig 4. The four figures are respectively the Aof - wave amplitude function, the Sa - occlusion degree function, the spatiotemporal diagram of gastric peristaltic waves, and the schematic diagram of 2D planar gastric peristaltic waves. It can be seen that the amplitude of the peristaltic wave reaches the maximum in the middle of the gastric antrum.

Introduction to Pyloric Mechanism
As the outlet of the stomach, under the set parameters, the motor activities of the gastric antrum, gastric outlet, and duodenum are synchronized. The pylorus acts as a "soft valve" that opens and closes periodically over time, acting only in a short axial range (pyloric segment). When the valve contracts, the effective radius of this pyloric segment decreases, almost blocking the passage from the stomach to the duodenum; when the valve relaxes, the radius returns to a larger value, and the gastric contents flow out in a short period of time. This cycle repeats, forming a typical "stepwise" gastric emptying rhythm. The pyloric cycle is divided into 4 phases:
(1) Contraction Rising Edge: Gradually shrinks from fully open to the target minimum radius.
(2) Contraction Holding Phase: Maintains a small radius for a long time (the valve is "almost closed"), with minimal emptying and strong inhibition of reflux.
(3) Relaxation Falling Edge: Quickly returns from a small radius to a larger radius (fast opening).
(4) Relaxation Holding Phase: Maintains an open state for a short time, allowing a wave of contents to pass through quickly.

Fig 6. The z-component velocity field of the model at t = 499 s. It can be seen that there is an obvious positive large-velocity flow at the pylorus. The results shown by the model are consistent with the classic function of physiological gastric emptying.

Governing Equations of the Pylorus We model the pylorus as a local radius modulation with a length of L1.5, expressed by a sine-squared function. Its instantaneous radius is simulated by defining the movement of the boundary mesh to realize the opening and closing of the pylorus:

The contraction ratio c(a) and its smooth step are:

The axial spatial weight is:

The initial state of the pylorus is relaxed, with a period of 20 s: it contracts gradually from 0 to 4 seconds, maintains contraction from 4 to 12 seconds, relaxes gradually from 12 to 18 seconds, and maintains relaxation from 18 to 20 seconds. Under the set parameters, the motor activities of the gastric antrum, gastric outlet, and duodenum are synchronized^[2]. The geometry of the computational domain changes periodically.

Jejunal Peristaltic Waves
Three different motion patterns are observed in the small intestine^[6] to promote food digestion and absorption. Peristalsis is responsible for the propulsion of chyme along the intestine; segmental contractions involve local contraction and relaxation without overall propulsion along the intestine, and the model studies their role in the mixing of chyme with surrounding fluids; pendular contractions involve the "back-and-forth" motion of the intestinal wall. The model mainly simulates the peristalsis and segmental contraction processes in the small intestine.
Governing Equations of Jejunal Peristaltic Waves
The peristaltic waves generated in the jejunum are simulated by defining the displacement of the intestinal wall. The jejunal peristaltic parameters are modeled as a Gaussian distribution function. The intestinal wall displacement and the step function used are defined as follows:

The relative amplitude amp of the peristaltic curve and the contraction period tc are determined based on literature data^[3].

Segmental Motion of the Jejunum and Duodenum Segmental motion refers to the quasi-steady, nearly local reciprocating contraction of the small intestinal wall on a short axial scale^[5]. Its main functions are to stir and mix the luminal contents and enhance mass transfer near the wall. Unlike peristaltic waves that mainly promote propulsion, segmental motion contributes little to the net propulsive flow but can significantly affect local shear, boundary layer thickness, and absorption/reaction environment. In the model, segmentation only acts in a preset short axial interval, and segmentation is regarded as a small radial displacement (the negative sign indicates contraction); its spatial shape is described by a sine function in the axial direction and a harmonic function in time; the intensity is clipped by a soft window to ensure smooth transitions at the ends and avoid numerical oscillations. In the model settings, segmental motion occurs simultaneously in the jejunum and duodenum, while peristaltic waves only occur in the jejunum.
(1) Governing Equations of Segmental Motion
Similar to peristaltic waves, the governing equation of segmental motion is the wall displacement defined on the wall. Segmental motion is regarded as a small radial displacement field oscillating sinusoidally with time in the interval^[5]:

Where Zs is the reference phase origin of the segmental wave in this segment, and W(z) is a spatial smoothing function (soft window). The soft window is defined by the difference between smooth steps at both ends:

The image characteristic of the segmentation function is a column of standing waves with phase changing uniformly with time.

Fig 9. Spatiotemporal diagram of duodenal segmental motion and schematic diagram of the change in the z-component velocity field caused by segmental motion in the 2D symmetric plane.

Model Introduction
The dilute substance transport model is the core module for quantifying the convection-diffusion behavior of the bacterial solution (using mol instead of CFU as the measurement unit, with consistent numerical calculation logic) in the intestine. Its core goal is to couple fluid flow (laminar flow) with the dynamic deformation of the intestine (moving mesh) to simulate the process of the bacterial solution from initial release, propulsion with intestinal contents, to radial diffusion. The model is based on the "dilute substance assumption" (the concentration of the product bacterial solution is extremely low, which does not change the fluid density and viscosity), only considering molecular diffusion and convective transport, and realizing the coordinated calculation with laminar flow and moving mesh through multi-physics coupling. Finally, it outputs evaluation indicators such as the spatiotemporal distribution of product concentration, centroid trajectory, and centroid velocity.

Model Assumptions and Core Objectives
(1) Model Assumptions
a. Dilute Substance Characteristics: The initial concentration of the bacterial solution is 1 mol/m³, which is much lower than the fluid saturation concentration. Its presence does not affect the density (1040 kg/m³) and dynamic viscosity (0.0014 Pa·s) of the intestinal contents, meeting the physical simplification conditions of "dilute substances".
b. Transport Mechanism: The transport of the bacterial solution in the intestine is dominated only by molecular diffusion and convection, and turbulent diffusion is ignored (since the Reynolds number of intestinal flow is Re ≈ 3, which belongs to the laminar flow category).
c. Initial Distribution: The bacterial solution is initially released in the form of an "elliptical bolus" at the jejunal axial position zc = 0.06125 m and radial position rc = 0 (near the axis). The semi-axes of the bolus in the axial (z) and radial (r) directions are az = 0.02 m and ar = 0.00692 m, respectively. The smooth function flc2hs is used to avoid numerical oscillations caused by hard edges (smoothing half-width sm = 0.00346 m).
d. Physics Field Coupling: The convective velocity field of the bacterial solution directly references the solution results of the laminar flow module (u_fluid, w_fluid), and the influence of mesh deformation on the concentration field is corrected by the displacement field (d(r, t), d(z, t)) of the moving mesh module to ensure the accuracy of concentration calculation in the dynamic intestinal environment.

(2) Core Objectives
a. Quantify the spatiotemporal concentration distribution of the bacterial solution under the synergy of small intestinal peristalsis, segmental motion, and gastric intermittent emptying mechanism, and evaluate the influence of intestinal motion on the mixing and propulsion of the product bacterial solution.
b. Calculate the axial centroid trajectory and axial dispersion coefficient of the bacterial solution, and verify the reproducibility of the model for the physiological function of "peristaltic propulsion - segmental mixing".
c. Provide the bacterial solution concentration boundary condition at the wall for the wall adsorption-desorption kinetics (boundary ODE module) to realize the full-process coupling of "transport - adsorption".

Governing Equations and Transport Mechanisms

The transport of the bacterial solution in the intestine follows the convection-diffusion equation. Considering the influence of moving mesh deformation on the concentration field (introducing the time derivative term of mesh displacement), its general form is as follows:

Where, Ji is the diffusion flux, c is the concentration of the bacterial solution, D is the molecular diffusion coefficient of the bacterial solution in the intestinal contents (m²/s), u is the velocity vector, and Rc is the reaction term. There is no chemical reaction in the main domain, and the adsorption-desorption reaction at the wall will be set separately through the boundary reaction, see Section 2.7 Boundary Ordinary Differential and Algebraic Equations.

Mesh Deformation Correction Term u

Due to the dynamic deformation of the mesh caused by the peristalsis and segmental motion of the intestinal wall, when COMSOL handles substance transport in the moving mesh, it will automatically introduce a mesh velocity term into the concentration equation to correct the calculation of the time derivative. The corrected concentration time change rate is:

ct is the intrinsic time derivative of concentration, cr and cz are the spatial gradients of concentration in the radial and axial directions, respectively. (∂r/∂t) and (∂z/∂t) are the radial and axial displacement rates of the mesh output by the moving mesh module. This correction ensures that the mesh deformation does not cause "spurious diffusion" in the concentration calculation.

Parameter Definition

See Section 2.1.1

Initial Concentration Distribution

It is assumed that the product passes through the stomach and duodenum under the protection of the outer coating and is released exactly at the entrance of the jejunum. The initial concentration of the bacterial solution is defined by an analytical function, adopting the form of "elliptical constraint + smooth step" to avoid numerical instability caused by sudden changes in the initial concentration. The expression at t = 0 s is as follows:

flc2hs is a built-in smooth Heaviside function in COMSOL. When the value in the parentheses is ≥ 0 (within the elliptical range), the function value approaches 1, and the concentration of the bacterial solution is C0 (C0 = 1 [mol/m³]); when the value is < 0 (outside the elliptical range), the function value smoothly transitions to 0, ensuring that there is no hard boundary in the initial concentration field.

Fig 10. Spatial distribution of product concentration in the 2D symmetric plane at t = 0, 17, 30 s.

Boundary Condition Settings

According to the physiological function of the intestine and the coupling requirements of the model physics field, the boundary conditions of the dilute substance transport module are divided into 3 categories, covering all geometric boundaries.

(1) No-Flux Boundary

A no-flux boundary is set on the solid walls without bacterial solution adsorption/release, such as the gastric wall, duodenal wall (non-adsorption region), and pyloric wall. The diffusion flux is 0, that is, the concentration gradient perpendicular to the wall is 0.

This means that the bacterial solution cannot diffuse through these walls and only moves with the fluid along the direction parallel to the wall or diffuses radially.

(2) Open Boundary

Open boundary conditions are set at the inlet and outlet of the overall digestive tract model, i.e., the end of the jejunum and both sides of the gastric antrum opening. The external concentration is set to 0.

To avoid reflux errors caused by sudden concentration changes at the boundary, the "upwind" numerical scheme is coupled, and the concentration value in the inflow direction is preferentially used for calculation to ensure flux continuity.

(3) Surface Reaction Boundary

A surface reaction boundary is set on the jejunal wall to simulate the intestinal mucosa region where the bacterial solution can be adsorbed. The flux of the bacterial solution is determined by the "adsorption-desorption kinetics". The dilute substance transport module only provides the bacterial solution concentration cw at the wall, and the flux calculation is dominated by the boundary ordinary differential and algebraic equation module.

The net flux of the bacterial solution at the wall is equal to the difference between the adsorption flux （Jads） and the desorption flux （Jdes）. The specific kinetic equation is detailed in Section 2.7 (Boundary ODE Module).

Multi-Physics Coupling

The dilute substance transport module needs to be directly coupled with the laminar flow and moving mesh modules, and the coupling with the boundary ordinary differential and algebraic equations will describe the bacterial solution adsorption-desorption process (see Section 2.7). The coupling logic is realized through the "Reactive Flow, Dilute Species" in the "Multiphysics" of COMSOL. The specific methods are as follows:

(1) Coupling with Laminar Flow

The "Fluid Flow (spf)" interface is adopted, and the velocity field (u_fluid, w_fluid) solved by the laminar flow module is directly substituted into the convection term u*(c∇) of the convection-diffusion equation.

(2) Coupling with Moving Mesh

d(z, t) of the moving mesh module is substituted into the concentration time derivative correction term:

Numerical Calculation
The "linear Lagrangian element" is used to match the velocity discretization of the laminar flow module to ensure the stability of the convection term calculation. Consistent with the overall transient settings, an adaptive time step is adopted, with a maximum step size of 0.1 s and a minimum step size of 1×10⁻⁴ s. The convergence criterion is "concentration absolute tolerance ≤ 0.001". The PARDISO direct solver is used to handle the sparse matrix of the concentration field to ensure the solution accuracy in the high-gradient region near the boundary layer.

Indicator Settings
The transport characteristics of the bacterial solution are analyzed from three dimensions: concentration spatiotemporal distribution, centroid trajectory, and axial dispersion coefficient to verify the physiological rationality of the model. See Section 3. Results and Analysis for the calculation result analysis.

Overview of the ODE/DAE Module
The boundary ODE/DAE module is the core module for quantifying the interaction between the bacterial solution and the intestinal wall. Its core function is to describe the dynamic equilibrium process of "adsorption-desorption" of the bacterial solution at the wall adsorption sites after reaching the intestinal wall in the dilute substance transport module, based on the Langmuir adsorption model and shear-enhanced desorption mechanism. This module only acts on the preset adsorption boundary (jejunal wall, see Section 2.6.5.1 Surface Reaction Boundary). Through coupling with the dilute substance transport module and laminar flow module, it realizes the quantitative transfer of "bacterial solution concentration → adsorption flux → wall adhesion amount", and finally outputs key indicators such as total wall adhesion amount, coverage area ratio, and net flux, providing a direct basis for evaluating the retention and distribution effect of the product bacterial solution in the intestine and a theoretical basis for the study of drug dosage and action sites.
Variables, Core Objectives, and Key Assumptions
(1) Variable Settings

Name	Expression	Description
Gamma	-	Wall adsorption amount
tauw_alt	spf.mu*spf.sr	Viscosity × shear rate (Pa)
koff	koff0*(1 + alpha*abs(tauw_alt))	Shear-enhanced desorption
theta	min(1,max(0,Gamma/Gmax))	Adsorption site occupancy rate
Jads	kon*c*(1 - theta)	Adsorption rate
Jdes	koff*Gamma	Desorption rate
M_Gamma	intop3( 2*pi*r*Gamma )	Total wall adhesion amount
CoverAreaFrac	intop3( 2*pi*r*flc2hs(Gamma/Gmax - 0.2, 1e-3) ) / intop3( 2*pi*r )	Coverage area ratio: the proportion of the area where Γ reaches the threshold (the threshold can be changed, e.g., 20%)
z_bar	intop3( 2*pi*r*z*Gamma ) / intop3( 2*pi*r*Gamma )	Coverage center position (where it is most concentrated)
J_net	intop3( 2*pi*r*(Jads - Jdes) )	Net flux (current intensity of deposition - desorption)
M0	intopD( 2*pi*r*an3(r,z) )[mol/m^3]	Initial total dosage (mol)
AttachYield	M_Gamma / M0	Adhesion yield (how much of the current dose remains)

(2) Core Objectives
Establish a kinetic model for the adsorption-desorption of the bacterial solution on the intestinal wall, and quantify the adsorption rate (Jads), desorption rate (Jdes), and net flux (Jnet); reveal the enhancing effect of intestinal wall shear stress on the desorption process, and clarify the influence of the shear sensitivity coefficient (α) on the adhesion stability; calculate the spatiotemporal distribution of the wall adsorption amount (Γ) in the "wall domain", and evaluate the coverage range and retention efficiency of the bacterial solution in the intestine.

(3) Key Assumptions
a. Monolayer Adsorption: The adsorption sites on the intestinal wall are limited and uniformly distributed. The adsorption of the bacterial solution follows the Langmuir model, with a maximum adsorption capacity of Gmax (no multi-layer adsorption).
b. Shear-Enhanced Desorption: The wall shear stress (τw) generated by the flow of intestinal contents will break the binding between the bacterial solution and the wall, leading to an increase in the desorption rate with the increase of shear stress. The quantitative relationship is related through the shear sensitivity coefficient α.
c. Non-Renewable Adsorption Sites: Within the model simulation duration (600 s), it is assumed that the adsorption sites on the intestinal wall do not regenerate or desorb, and only the cycle of "unoccupied sites → adsorption → desorption → empty sites" is considered.
d. Coupling Dependence: The adsorption rate depends on the bacterial solution concentration at the wall (cw) provided by the dilute substance transport module (see Section 2.6.5.3), the desorption rate depends on the wall shear rate (sr) provided by the laminar flow module, and mesh deformation does not affect the distribution of adsorption sites (the moving mesh module only corrects the geometric position and does not change the wall adsorption characteristics).

Governing Equations and Kinetic Mechanisms

(1) Core Equation of Adsorption-Desorption Kinetics
The adsorption-desorption process of the bacterial solution on the intestinal wall satisfies "dynamic equilibrium". The time change rate of the wall adsorption amount Γ (Gamma) is equal to the net flux (Jnet), which is the core governing equation of the boundary ODE:

Where, Jads is the adsorption flux, describing the migration rate of the bacterial solution from the fluid phase to the wall phase; Jdes is the desorption flux, describing the detachment rate of the bacterial solution from the wall phase to the fluid phase.

(2) Adsorption Flux: Langmuir Adsorption Model
The adsorption process is controlled by both the "bacterial solution concentration at the wall" and the "proportion of unoccupied adsorption sites".

Where, kon is the adsorption rate constant (m/s), reflecting the binding ability between the bacterial solution and the wall sites; cw is the concentration of the bacterial solution at the wall, provided in real-time by the dilute substance transport module, i.e., the calculation result of the bacterial solution concentration c at the jejunal wall; θ is the adsorption site occupancy rate; Gmax is the maximum adsorption capacity.

(3) Desorption Flux: Shear-Enhanced Mechanism
The desorption process is divided into two parts: "basal desorption" and "shear-enhanced desorption". Mechanical shear will enhance desorption (shear force promotes the detachment of the bacterial solution from the intestinal wall).

Where, koff0 is the effective desorption rate, which changes with shear stress; koff0 is the basal desorption rate, i.e., the desorption ability without shear; is the shear sensitivity coefficient, quantifying the enhancing effect of shear on desorption; τw is the wall shear stress (Pa), which is approximately calculated in the laminar flow module in the model.

Where, μ is the fluid dynamic viscosity, and sr is the shear rate.

Initial Conditions
At the initial moment of the model (t = 0 s), there is no bacterial solution adsorption on the intestinal wall, and the initial values of the adsorption amount and adsorption rate are set to 0, which is consistent with the physiological initial state of "the bacterial solution has not reached the wall".

Boundary Selection and Multi-Physics Coupling Mechanism

(1) Applicable Boundary
Only the jejunal wall is selected as the boundary. The jejunum is the core region for intestinal nutrient absorption and bacterial colonization. Selecting this segment of the wall to simulate adsorption is consistent with the product's target site requirements.
(2) Coupling Logic
It needs to be directly coupled with the dilute substance transport and laminar flow modules. The coupled data is realized through the "variable transfer" and "multi-physics interface" of COMSOL. The specific logic is as follows:

Coupled Module	Transferred Data	Transfer Direction	Physical Meaning	Implementation Method
Dilute Substance Transport (tds)	Wall bacterial solution concentration cw	Dilute Substance Transport → ODE/DAE	Provide the concentration driving term for the adsorption flux Jads	Reference the concentration calculation result c at the jejunal wall boundary of the dilute substance transport module
Dilute Substance Transport (tds)	Net flux Jnet	ODE/DAE → Dilute Substance Transport	Provide the wall flux boundary condition for the dilute substance transport	Substitute Jnet into the jejunal surface boundary of the dilute substance transport
Laminar Flow (spf)	Wall shear rate sr	Laminar Flow → ODE/DAE	Provide the shear driving term for the effective desorption rate koff	Reference the shear rate variable sr of the laminar flow module, and calculateτw​ in combination

(3) Verification of Coupling Equations
After coupling the dilute substance transport module, the complete format of the bacterial solution flux balance equation at the jejunal wall boundary is:

The left side is the wall diffusion flux of the dilute substance transport module (positive when pointing to the wall), and the right side is the net adsorption flux of the ODE/DAE module. The equality of the two ensures the mass conservation of the "fluid phase - wall phase", verifying the rationality of the coupling logic. n*(-D∇c)represents the diffusion of the bacterial solution through the jejunal wall, i.e., entering the "wall domain" from the "intestinal domain" or vice versa.

Numerical Calculation Settings
(1) Discretization: The shape function adopts the "discontinuous Lagrangian" quadratic element.
(2) Linear Solver: The PARDISO direct solver of the overall model is used to handle the sparse matrix equation of the adsorption amount Γ, avoiding the convergence risk of iterative solution.

Indicator Settings

Based on the transient simulation results of COMSOL Multiphysics 6.3 (600 s, output time step of 0.1 s, solution time: 13188 s), combined with the physiological function of the intestine, the analysis is carried out from three core dimensions: the hydrodynamic characteristics of the model, the transport behavior of dilute substances, and the wall adsorption kinetics. The influence of the synergistic effect of peristalsis and segmentation on product migration, mixing, and wall retention is quantified, the reproducibility of the model for physiological processes is verified, and data support is provided for the prediction of product action sites and dosage optimization.

Definition of Global Flow Rate and Monitoring Position
To test the ability of the model to reflect the basic characteristics of fluid transport in the digestive system, global flow rate variables such as inflow (jejunal inflow), outflow (jejunal outflow), infloww (gastric inflow into the pylorus), and outfloww (pyloric inflow into the duodenum) are defined in the intestinal model. These variables are the core indicators for quantifying fluid retention, propulsion efficiency, and inter-regional mass exchange in the intestine. These variables are defined through the "Global Equations" and "Integration Coupling" functions of COMSOL, directly related to the fluid flow (laminar flow module) and geometric boundaries (such as the jejunal inlet/end, pyloric interface). Their dynamic changes can verify the mass conservation of the model, reflect the regulatory effect of peristalsis and segmentation on the flow rate, and provide data support for the "fluid-driven mechanism" of product transport.

Global Variable Name	Expression	Monitoring Boundary	Physical Meaning
inflow (Jejunal Inflow)	intop1(2pir*w_fluid)	Jejunal Inlet	Instantaneous volume flow rate of fluid entering the jejunum from the duodenum/pylorus
outflow (Jejunal Outflow)	intop2(2pir*w_fluid)	Jejunal Terminal Outlet	Instantaneous volume flow rate of fluid flowing out of the model from the jejunum
infloww (Gastric Inflow into Pylorus)	intop5(2pir*w_fluid)	Stomach-Pylorus Interface	Instantaneous volume flow rate of fluid entering the pylorus from the gastric antrum
outfloww (Pyloric Inflow into Duodenum)	intop4(2pir*w_fluid)	Pylorus-Duodenum Interface	Instantaneous volume flow rate of fluid entering the duodenum from the pylorus
netflow (Jejunal Net Flow)	netflowt - (outflow + inflow)/2	Global Integration (Jejunal Domain)	Cumulative net retention of fluid in the jejunum (inflow - outflow)
netfloww (Pyloric Net Flow)	netflowwt - (outfloww + infloww)/2	Global Integration (Pyloric Domain)	Cumulative net retention of fluid in the pylorus (gastric inflow - duodenal outflow)

Flow Direction: w_fluid is the axial (z-direction) velocity component. A positive value indicates the positive direction along the z-axis (from stomach → pylorus → duodenum → jejunum). Therefore, a positive value of inflow/infloww represents fluid "entering" the target region, and a positive value of outflow/outfloww represents "flowing out".

Dynamic Change Law of Flow Rate
(1) Flow Rate Change in the Duodenal Domain

Further, the mass conservation of the model and the physiological rationality of the flow pattern are verified through the cumulative inflow and outflow curves at the jejunal site. The cumulative volume curve of the jejunum shows a typical "step + fine line" shape: the stepwise jump synchronized every 20 s reveals the intermittent emptying driven by the pyloric gating; the 6 s fine lines inside the steps result from the high-frequency modulation of the instantaneous flux by segmental motion. Reflux only occurs at the moment of valve switching with a limited amplitude, indicating that the valve control has a significant inhibitory effect on reflux. The slow fluctuation of the step height is consistent with the amplitude time course of the gastric antral peristaltic wave train (lifespan ≈ 61 s). The three curves satisfy the net volume relationship and rise approximately in parallel, indicating that the mass conservation is consistent with the boundary flux evaluation.
(2) Flow Rate Change in the Pyloric Domain

The cumulative efflux (blue), cumulative reflux (green), and cumulative net emptying (red) at the pylorus show a clear stepwise rise within 0–600 s. The step rhythm is completely in phase with the 20 s cycle "soft valve" gating in the model: the gating is applied to the pyloric segment with a length of L1.5 through the boundary displacement function, and the axial weight of sin2(ℼZ/L1.5) is used to limit the action range. Thus, the cumulative curve is nearly horizontal during the "contraction holding" period of the valve and triggers a stepwise jump during the "relaxation holding" period, reproducing the classic rhythm of intermittent gastric emptying. Fine undulations of approximately 6 s can be seen inside each step, which is consistent with the set segmentation frequency of the model (f=10 min-1), indicating that the duodenal/jejunal segmental motion causes high-frequency modulation of the instantaneous flux and near-wall shear rather than decisive propulsion. It can be obtained that at 600 s, the cumulative efflux volume is 2.3*10-4 m3, the cumulative reflux is 0.65*10-4 m3, and the cumulative net emptying is 1.45*10-4 m3. The reflux only accounts for 30% of the efflux, indicating that the gating and geometric contraction have a significant inhibitory effect on reflux; the estimated average net emptying rate is 2.4*10-7 m3s-1 approximately . The three curves maintain the relationship of "red ≈ blue − green" for a long time and rise approximately in parallel, indicating that the boundary integral operator used has good mass conservation and numerical consistency; at the same time, the slow envelope of the jump amplitude of each step fluctuates with time, which is consistent with the gastric antral peristaltic wave train lifespan (≈ 61 s) and the time course modulation of the amplitude envelope Aof, reflecting the synergistic effect of gastric antral peristalsis, pyloric gating, and segmental motion on multiple time scales.

Dynamic Change Law of Velocity Field

At t = 415.2 s (pyloric relaxation window), the local velocity field shows a typical "central jet + bilateral reflux vortex" structure: the jet is ejected from the narrowest part of the pylorus and forms a long and narrow high-velocity zone at the duodenal inlet, then decays rapidly along the radial direction; a pair of near-wall counter-rotating vortex cells are induced by the shear layer on both sides of the jet, and the streamlines are closed in the vortex cells, indicating that this region is mainly for mixing and recirculation and contributes little to net propulsion. The large shear gradient at the lip of the jet means that the wall shear stress reaches a relatively high value, which is consistent with the aforementioned division of labor of "segmentation enhances mixing on a short scale, and gating triggers jet transport during the window period", and explains the source of stepwise emptying and high-frequency fine lines in the cumulative curve from the perspective of the instantaneous field.

The flow field predicted by this model is highly consistent with the classic description of gastric motion when the contents are Newtonian fluids (0.0014 Pa·s, 1040 kg/m³). In particular, the flow field predicted in the pyloric antrum region shows two main flow patterns. When the gastric antral contraction wave moves toward the pylorus, the gastric contents are forced to reflux. When the gastric occlusion degree reaches 56%, a retrograde jet-like motion is formed. To quantify the development process of this jet motion, the maximum retrograde flow velocity was tracked within a 20-second flow cycle. As the gastric antral contraction wave approaches the pyloric sphincter, the increase in contraction wave velocity and occlusion degree causes the retrograde jet to accelerate continuously, reaching a maximum retrograde flow velocity of 8.8 cm/s at the narrowest part of the pyloric canal. Another main flow pattern identified in the pyloric antrum region is the annular vortex formed between two adjacent gastric antral contraction waves (Figure 14). Through strong mixing, friction, and grinding of the gastric contents, retrograde jet motion and vortex structures have long been considered the main drivers of gastric digestion of food^[1].

3.2 Transport Behavior of Dilute Substances

Spatiotemporal Distribution of Concentration

Fig 15. Shows the spatial distribution of product concentration in the 2D symmetric plane at t = 0, 150, 300, 450, 600 s.

Over time, under the combined action of peristaltic propulsion and segmental mixing, the concentration field gradually expands along the axial direction and diffuses radially. It shows monotonic dilution under the combined action of convection dominance, diffusion, and segmental mixing. At the same time, the attenuation of the bulk concentration is also continuously absorbed by the wall adsorption-desorption flux.

Centroid Position

This figure shows the time evolution of the axial centroid of the bolus. The curve shows a typical "step + fine line" characteristic: the large-scale steps are strictly in phase with the 20 s pyloric gating, indicating that dzc/dt increases significantly during the relaxation window and almost stagnates during the contraction holding period; the high-frequency fine teeth inside the steps are synchronized with the set segmentation frequency, reflecting the periodic modulation of segmental motion on near-wall mixing and instantaneous convection rather than decisive propulsion. The overall monotonic rise and no continuous retraction indicate that the net propulsion is dominated by gating, and reflux only occurs at the moment of opening and closing switching (consistent with the aforementioned "cumulative volume step"); the slight envelope fluctuation of the curve on a longer time scale is consistent with the modulation of the gastric antral peristaltic wave train lifespan and amplitude envelope Aof. In summary, the behavior of the centroid is highly consistent with the rhythm of volume efflux: gating determines the migration "jump distance", segmentation accelerates spreading and mixing, and adsorption/desorption only changes the bulk depletion and spatial distribution, but does not change the forward direction of the centroid.

Axial Dispersion Coefficient

The axial dispersion coefficient shows "a slowly increasing envelope + obvious periodic oscillation": within 0–600 s, its peak-to-peak amplitude and positive half-cycle average value continue to increase, indicating that the bolus advances into a region with stronger shear/stronger stirring (the central jet generated by pyloric windowing and downstream segmentation jointly enhance axial spreading); the large-scale oscillation in groups of 20 s is in phase with the pyloric gating: positive peaks appear in the relaxation window, and it is close to zero or negative valleys appear in the contraction holding period. The high-frequency fine lines (6 s) superimposed in the steps come from the modulation of segmental motion on instantaneous mixing. The negative value in the figure is not "physically negative diffusion", but "the instantaneous axial spreading rate is negative, which is consistent with the aforementioned "step + fine line" rhythm of centroid stepwise forward movement and cumulative volume.

Wall adsorption is the core link for the product bacterial solution to achieve retention and action in the intestine. Based on the Langmuir adsorption model and shear-enhanced desorption mechanism, key indicators such as total adhesion amount, coverage range, and net flux are used to quantify the dynamic behavior of the bacterial solution at the jejunal adsorption boundary, reveal the coupling law of "fluid shear - wall concentration - adsorption equilibrium", and verify the reproducibility of the model for the intestinal adsorption process. See Section 2.7.2.1 for evaluation indicators and calculations.

Total Wall Adhesion Amount

Within 0–600 s, the total wall adhesion amount MΓ increases monotonically with superimposed 20 s serrated oscillations; by 600 s, it has reached approximately 5.6×10−7mol, indicating that under the current dosage and parameter set, net adsorption is dominant for a long time, and the bulk phase is continuously "extracted" to the wall without obvious global saturation (Г/Gmax has not approached 1). The serrations are in phase with the pyloric gating: the relaxation window triggers a jet → the near-wall concentration c in the bulk phase increases sharply and the supply is sufficient, but at the same time, the shear increases to make koff larger; the shear decreases during the contraction holding period, the recirculation is enhanced, koff decreases, and Jads continues to work in the medium shear zone. The result of the alternation of the two phases is "scouring - redeposition" in each 20 s cycle, but the net increase is positive, resulting in the "sawtooth-like rise" shown in the figure. On a longer time scale, the slope rises slowly, which is consistent with the synergistic enhancement of near-wall supply and mixing by the gastric antral ACW wave train (61 s) + duodenal/jejunal segmentation: the stepwise forward movement of the bolus centroid periodically sends the high-concentration zone to the new downstream wall, and the product of the available area × near-wall concentration × contact time expands with time, so MΓ shows a nearly linear - weakly superlinear growth.

The coverage area ratio is "the proportion of the intestinal wall area that reaches the effective coverage threshold. Its formula is:

Where Ajej is the total area of the jejunal wall, Γ is the wall adhesion amount, Gmax is the maximum wall capacity, and​θthr is the set "effective coverage threshold" (usually 0.2, i.e., coverage ≥ 20% is regarded as "covered", that is, the adhesion amount reaching 20% of the maximum capacity is regarded as covered). Therefore, this indicator emphasizes the geometric coverage range rather than "how much quantity", which is suitable for answering where the drug/bacterial solution is distributed and how large the coverage area is, and is complementary to the total adhesion amount MГ​ ("how much quantity"). Figure 17 shows that the coverage area ratio increases monotonically within 0–600 s and presents a "step + fine line" shape: by 600 s, it has reached approximately 0.93, indicating that under the current administration and kinetic parameters, most areas of the jejunal wall have reached effective coverage. The steps of the curve are in phase with the 20 s pyloric gating: the jet pulse during the relaxation period quickly pushes the near-wall high-concentration zone downstream, triggering a large number of new walls to cross the threshold, and the coverage ratio jumps; the propulsion weakens during the contraction holding period, and the coverage ratio rises slowly/platforms. From the perspective of spatial distribution, the curve indicates that the model successfully reproduces the entire process of "pulsatile transport (gating determines the jump distance) - local mixing (segmentation determines spreading and near-wall supply) - wall capture".

Net Flux

The figure shows the global net flux Jnet, which can be understood as the instantaneous net rate of "adsorption - desorption" on the entire segment of the intestinal wall; it is the time derivative of the total wall adhesion amount MΓ​, so Jnet >0 indicates that the wall is increasing net. The curve in the figure shows three phases and is superimposed with periodic textures in phase with physiological driving forces:

(1) 0-60 s: Jnet jumps rapidly from zero, which comes from the initial bolus + the near-wall high-concentration jet supply generated by the first pyloric relaxation.

(2) 60-350 s: The 20 s sawtooth is in phase with the pyloric gating, and the 6 s fine teeth come from segmental mixing. The gating "batches" the high-concentration zone to the downstream, and the segmental motion enhances the near-wall supply. However, with the increase of the covered area and the continuous extraction of the bulk phase, the growth of Jads=konc(1-θ) is partially offset, thus entering a quasi-steady state limited by transport.

(3) 350-600 s: The bulk concentration further decreases, more walls approach the capacity limit, and the proportion of the desorption term koff=koff0(1+α|ﬢw|) in the high-shear window of the gating increases relatively. The combination of these factors leads to a decline in the net adsorption rate.

Adhesion Yield

The total initial administered dose (sum of the bulk phase and initial wall value) is denoted as M0​. The adhesion yield is defined as the normalized proportion of the cumulative retention on the wall, which measures "the proportion of the initial dose that has been captured by the wall by time t".

At t=600 s, AttachYield=0.14, indicating that approximately 14% of the initial dose has been transferred to the wall.

Given that distinct promoters may activate at different physiological conditions and that engineered modules often function sequentially, we designed simulations that mirror possible experimental subgroups. Six configurations were analyzed to assess metabolic burden:

(1) Adsorption module (2) Degradation module (3) Therapeutic module (4) Adsorption + Degradation + Therapeutic modules (5) Adsorption + Therapeutic modules (6) Degradation + Therapeutic modules

Flux Balance Analysis (FBA) is a standard computational framework in metabolic modeling. The general linear constraint equations are as follows:

(1) Steady-State Assumption The intracellular metabolite concentrations are assumed to remain constant under steady-state conditions, meaning that the production rate of each metabolite equals its consumption rate. Mathematically, this can be expressed as:

where S is the stoichiometric matrix of the metabolic network (rows represent metabolites, columns represent reactions), and v is the reaction flux vector.

(2) Flux Constraints Each reaction flux is bounded by upper and lower limits that define its reversibility and feasible rate range:

where li and ui​ represent the lower and upper bounds of reaction i, respectively. These constraints incorporate both thermodynamic directionality and experimentally determined limits on transport or enzyme capacity.

(3) Objective Function FBA typically assumes that cells have evolved to optimize certain physiological objectives, such as maximizing growth rate or biomass formation. The objective function is therefore defined as:

where c is a coefficient vector selecting the target flux .Solving this linear programming problem yields an optimal flux distribution v* and the corresponding optimal growth rate.

To present our results more intuitively, we visualized the simulation data as follows:

(1) Growth rate vs. expression level X-axis: Protein expression intensity (log scale) Y-axis: Growth rate (h⁻¹) Dashed line: Wild-type control

(2) ATP burden vs. expression level X-axis: Protein expression intensity (log scale) Y-axis: ATP consumption (mmol/gDW/h) A higher value indicates a greater energetic burden.

(3) GTP burden vs. expression level Same as the ATP burden plot, but representing GTP consumption.

A daily dietary supplement designed for people at high risk of microplastic exposure and health-conscious individuals (which can later be registered as a special medical purpose formula food).

Variables	Biological meanings
Ste2	Concentration of the receptor Ste2
α	Concentration of α-factor
Ste4-Ste18	Concentration of the βγ subunits dissociated from the G protein
α-Ste2	Concentration of the α-factor–Ste2 receptor complex
α-Ste2-G	Concentration of the complex formed by α–Ste2 and the G protein
Ste5-Ste4-Ste18	Concentration of the Ste5–Ste4–Ste18 complex
Ste5-Ste4-Ste18-Ste11	Concentration of the Ste5–Ste4–Ste18–Ste11 complex
Ste5-Ste4-Ste18-Ste11-Ste7	Concentration of the Ste5–Ste4–Ste18–Ste11–Ste7 complex
Fus3	Concentration of Fus3 (MAPK)
Ste12	Concentration of the transcription factor Ste12
αmRNA	Concentration of mRNA transcribed from the α-factor gene
enzyme	Concentration of the translated degradation enzyme

Parameters	Biological meaning
k₁	Dissociation rate constant of α–Ste2
k₂	Association rate constant between α-factor and Ste2
k₃	Dissociation rate constant of α–Ste2–G
k₄	Association rate constant between G protein and α–Ste2
k₅	Association rate constant between Ste5 and Ste4–Ste18
k₆	Dissociation rate constant of Ste5–Ste4–Ste18
k₇	Association rate constant between Ste5–Ste4–Ste18 and Ste11
k₈	Dissociation rate constant of Ste5–Ste4–Ste18–Ste11
k₉	Association rate constant between Ste5–Ste4–Ste18–Ste11 and Ste7
k₁₀	Dissociation rate constant of Ste5–Ste4–Ste18–Ste11–Ste7
k₁₁	Activation rate of Fus3
k₁₂	Activation rate of Ste12
k₁₃	Transcription rate of the α-factor gene
k₁₄	Degradation rate of α-mRNA
k₁₅	Translation rate of the α-factor gene
k₁₆	Degradation rate of α-factor
k₁₇	Transcription–translation rate of the degradation enzyme
k₁₈	Feedback regulation coefficient of Ste12
k₁₉	Inhibition coefficient of Fus3 under excessive α-factor

Based on literature data and known biological properties, the corresponding parameter values were determined and integrated into the simulation model The model primarily simulates the following four processes:

(1) Receptor activation stage
Chemical reaction:

ODE：

Simulation Results：
The concentrations of both α-factor and the α–Ste2 complex increased significantly, indicating that our system can not only secrete α-factor to participate in quorum sensing, but also effectively enhance the expression of downstream enzyme genes.

(2) MAPK Cascade Signal Amplification
Chemical reaction:

ODE：

Simulation Results：
The concentrations of intermediate signaling proteins in the MAPK cascade first increased and then decreased, consistent with the typical transient activation pattern observed in signaling intermediates. In particular, the elevated concentration of Fus3 suggests its promoting effect on the transcription of downstream target genes.

(3)α-Factor Secretion and Positive Feedback Formation
Chemical reaction:

ODE：

Simulation Results：
Simulation results showed an overall increase in the concentrations of the upstream transcription factor, α-factor mRNA, and α-factor itself. The rate of α-factor accumulation first increased and then slowed down. This is because, in the early stage, the positive feedback of α-factor synthesis and secretion accelerates its concentration growth. However, as the extracellular α-factor concentration becomes high, excessive Fus3 exerts a negative regulatory effect on Ste12 concentration^[3].

To represent this effect, a negative regulatory term was introduced into Equations (9) and (11) of the ODE system, reflecting Fus3-mediated inhibition of Ste12 and α-factor expression. This adjustment prevented an unrealistic continuous increase of α-factor concentration and made the simulation results more consistent with biological observations.

(1) he complete amino acid sequence (228 amino acids) of the human MXl1 protein was extracted from the UniProt-provided FASTA sequence.

(2) Motif scanning analysis

The SIN3-binding signature was Scanned using the regular expression pattern: [LIVFMWY][^P][^P][LIVFMWY][LIVFMWY].

Two potential SIN3-binding motifs were identified in the N-terminal region:

Motif 1: Position 9-13, sequence "MINVQ"

Motif 2: Position 90-94, sequence "LKVLI"

To systematically visualize and validate the interactions among core components of the transcriptional complex, a protein-protein interaction (PPI) network was constructed.

This network comprises five key proteins: MXI1, SIN3A, SIN3B, HDAC1, and HDAC2. Each node represents a protein, and its unique color and internal structure diagram facilitate visual distinction. The numerous connection edges (lines) between all nodes reveal a dense and highly interconnected network structure, indicating strong and multi-faceted interactions within the complex. Notably, the MXI1 node (red) has direct connections with both SIN3A (yellow) and SIN3B (green), providing computational evidence for our initial hypothesis and ChIP-seq findings: MXI1 acts as a key recruiter by directly binding to the scaffold proteins SIN3A and SIN3B. Additionally, Sin3 proteins exhibit direct interactions with histone deacetylases HDAC1 (light blue) and HDAC2 (blue), which is consistent with their known role in forming the core of the repressive complex. The presence of HDAC1 and HDAC2 in this network directly links our complex to the enzymatic activities responsible for chromatin condensation and transcriptional silencing.

Genome-wide mapping of the MXI1-Sin3a/b-HDAC1/2 complex revealed significant enrichment on gene regulatory elements. Annotation of the co-binding sites indicated that 92.3% of the co-binding peaks were located within promoter regions, with the vast majority (92.3%) within 1 kb upstream of the transcription start site (TSS). This significant promoter preference is highly consistent with the established role of Sin3-HDAC complexes as direct regulators of transcription initiation.

To investigate the functional state of chromatin at these binding sites, we evaluated the enrichment of the active enhancer mark histone H3 lysine 27 acetylation (H3K27ac). Notably, we found that the H3K27ac signal at MXI1-Sin3B-HDAC1 co-binding sites was 29.6 times higher than that at randomly selected genomic regions.

This difference was statistically highly significant (Wilcoxon rank sum test, p < 0.05). This finding demonstrates that a complex known for its repressive function is associated with an active chromatin mark - suggesting a complex model of gene regulation. It strongly implies that the MXI1-Sin3B-HDAC1 complex is not recruited to silent loci but is targeted to already active promoters to fine-tune or repress transcriptional output.

To validate the practical relevance of our molecular findings in a clinical context, we utilized cBioPortal to analyze the co-expression patterns of the components of the MXI1-Sin3-HDAC complex in a large number of human cancer samples from the TCGA database. Our analysis revealed a consistent and statistically significant positive correlation between MXI1 and its key partners. Specifically, the mRNA expression of MXI1 showed a strong positive correlation with SIN3A (Spearman's ρ = 0.39, p < 1.63e-16), and moderate positive correlations with HDAC1 (ρ=0.25, p=2.70e-7) and HDAC2 (ρ=0.26, p=8.46e-8).

In the protein screening, we conducted molecular docking between the expected HFBI protein to be screened and AGA1. The results of the docking provided support for the next step of the experiment. The specific results are as follows: