Description of Analog and Digital Work

Our modeling work follows the logical chain system of "light-flow-bacteria":

1. First, the Mie scattering-refractive index fluctuation coupling model was constructed to quantify the dynamic light intensity attenuation in the environment of Archimedes screw weak turbulence, so as to provide accurate dose for subsequent optical regulation;
2. The 7th room diffusion-mass transfer equation was used to describe the cross-barrier mass transfer process of engineering bacteria in fish, and the light-controlled bacterial quantity was transformed into the immune endpoint of chamber balance through model language;
3. Design the green/blue/red tricolor light gene circuit under NIR activation, and use the two-component kinase-response regulation network to realize the transient, reversible and saturation control of swimming mode of E. coli.

The three modules are seamlessly connected, covering three major scales: photon, fluid, and microbial species, to comprehensively explain the vaccine immunization results of our developed screw device. In model construction, we employed various numerical methods including Monte Carlo, Runge-Kutta, finite difference, Latin hypercube sampling, and Sobol global sensitivity analysis, combined with high-quality heat maps, response surface diagrams, and convergence analysis charts to enhance the model's intuitiveness and credibility. By introducing time-varying optical path length L(t), we helped understand how mechanical structures modulate light transmission paths; elucidated how optical signals regulate CheY phosphorylation states through CcaS/CcaR or YF1/FixJ systems, providing mechanistic insights into light-controlled biological systems. The system data for model construction were sourced from literature and experimental data, with parameters calibrated through literature references. This ultimately established a "light-bio-fluid" system modeling framework which is worth generalizing, offering researchers in related fields comprehensive methodological guidance, code, and data references.

Modeling and Analysis of Light Transmission Characteristics Based on Mie Scattering and Refractive Index Fluctuation Model

This model mainly constructs a mathematical model of light intensity attenuation for Edwardsiella vaccine solution in a weak turbulence system to analyze and predict the transmission behavior of light in the Edwardsiella vaccine solution. The model combines the Mie scattering theory (for particle scattering) and the refractive index fluctuation model (for turbulence scattering), and extends the model to a time-varying form for the dynamic environment of the Archimedes screw system. Numerical methods such as Monte Carlo and Runge-Kutta are used for solving, the light intensity attenuation characteristics are systematically analyzed, and the key influencing parameters are determined using Sobol global sensitivity analysis. The research provides important theoretical basis and design guidance for optimizing the optical detection and therapeutic technologies of vaccines.

1.Overview of Model Meaning

The study of light transmission characteristics in turbid media is an important topic in fields such as bio-optics, atmospheric physics, and ocean optics. Especially in the Edwardsiella vaccine solution system, the light transmission process is comprehensively affected by microbial particle scattering and turbulence effects, showing complex attenuation and polarization characteristics. The Mie scattering theory provides a strict mathematical basis for analyzing the scattering phenomenon when the particle size is comparable to the wavelength of light, while the refractive index fluctuation model well describes the random scattering effect caused by turbulence. This paper will combine these two theories to establish a complete light transmission model, solve it using numerical methods, and finally display the variation law of light transmission characteristics with various parameters through visualization methods. This study aims to construct a mathematical model that can accurately predict the transmission behavior of light in the Edwardsiella vaccine solution, providing a theoretical basis for optimizing vaccine optical detection methods and developing new optical therapeutic technologies. The model not only considers the scattering and absorption of microbial particles but also introduces the influence of turbulence effects, making the model more in line with practical application scenarios. Through systematic parameter research and visual analysis, this paper reveals the key physical mechanisms and laws in the light transmission process.

2 Model Establishment and Mathematical Expression

2.1 Mie Scattering Theory Model

The Mie scattering theory is a strict mathematical solution describing the scattering effect of uniform spherical particles on electromagnetic waves, whose core is to solve the boundary condition problem of Maxwell's equations in the spherical coordinate system. For microbial particles in the Edwardsiella vaccine solution, when their size is comparable to the wavelength of near-infrared light, the Mie theory can describe the light scattering phenomenon more accurately than the Rayleigh scattering theory.Scattering efficiency factor is a core parameter of the Mie theory, expressed as:

Among them, is the size parameter, is the particle radius, and is the wavelength of incident light.​and ​are the Mie coefficients, which are closely related to the complex refractive index of the particles.

​​The scattering phase function describes the angular distribution characteristics of light after scattering, and is approximately expressed using the Henvey-Greenstein phase function:

Among them, is the scattering angle, and g is the anisotropy factor, with a value range of [-1, 1], representing the degree of forward or backward preference of scattering.

For the Edwardsiella vaccine solution, the scattering coefficient ​can be obtained through integral calculation:

Among them, N(r) is the particle size distribution function, and is the scattering cross-section.

2.2 Refractive Index Fluctuation Model

2.2.1 Theoretical Basis and Wave Equation Derivation

The propagation of light waves in turbulent media follows the basic laws of electromagnetic waves and can be derived from Maxwell's equations. We first make the following assumptions:

1. The medium is a perfect insulator (conductivity is 0);
2. The magnetic permeability is 1;
3. The wavelength of the electromagnetic wave is much smaller than the inner scale of turbulence;
4. The time dependence of the electromagnetic field is determined by ;
5. 5.The statistical characteristics of the refractive index do not change with time.

Considering the random fluctuation characteristics of the refractive index, the refractive index is expressed as:

Among them, ​ is the average refractive index, and is the fluctuating component of the refractive index, which satisfies .

Starting from Maxwell's equations, the vector wave equation can be derived:

Among them, is the wave number in vacuum.

When the inner scale of turbulence ​is much larger than the wavelength ，the last term in the above equation is a high-order infinitesimal quantity and can be ignored. At the same time, the depolarization effect of the atmosphere can usually be ignored, so the vector equation can be simplified to a scalar equation:

2.2.2 Rytov Approximation and Perturbation Solution Method

For solving the above wave equation, the Rytov approximation (gentle perturbation method) is usually used. Let ，where is the complex phase, expressed as:

Here, is the logarithmic amplitude fluctuation, and is the phase fluctuation.

After substituting into the wave equation, the Riccati equation is obtained:

The perturbation method is used for solving, and is expanded as , where ​corresponds to the solution under the condition of no turbulence, and ​is the first-order perturbation term. For the case of plane wave incidence, .

The equation satisfied by the first-order perturbation term is:

The solution of this equation can be obtained using the Green's function method:

2.2.3 Refractive Index Fluctuation Statistics and Power Spectrum Model

The statistical characteristics of the refractive index fluctuation are usually described by its power spectrum function . Under the assumption of homogeneous and isotropic turbulence, there is a definite relationship between the power spectrum of refractive index fluctuation and the structure function of turbulence.

According to the Kolmogorov turbulence theory, the structure function of refractive index fluctuation is:

Among them, is the refractive index structure constant, which characterizes the turbulence intensity. The corresponding three-dimensional power spectrum function is:

In the actual atmospheric environment, it is necessary to consider the influence of the inner scale and outer scale of turbulence. The commonly used modified models include the Tatarski spectrum and the Von Karman spectrum:

Tatarski spectrum:

Among them, .

Von Karman spectrum:

Among them, .

2.2.4 Statistical Characteristics of Light Field Under Weak Turbulence Conditions

In the context of this study, the flow of the Edwardsiella vaccine solution in the Archimedes screw immersion device with a rotation speed of 3 r/min is regarded as a weak turbulence situation.

Under weak turbulence conditions (Rytov variance ), the first-order solution of the Rytov approximation can be used to derive the statistical characteristics of the light field. The variances of the logarithmic amplitude fluctuation χ and the phase fluctuation S can be expressed as follows:

For plane waves:

For spherical waves:

The relationship between the variance of light intensity fluctuation (scintillation) and the logarithmic amplitude fluctuation is:

2.2.5 Derivation of Turbulence Scattering Coefficient

Based on the above theory, we can derive the scattering coefficient caused by turbulence. The turbulence scattering coefficient is defined as the proportion of light intensity attenuated by turbulence scattering per unit distance. Considering the light intensity expression and combining the statistical characteristics of light intensity fluctuation under the Rytov approximation, the relationship between the turbulence scattering coefficient and the refractive index structure constant can be obtained:

This expression reflects the dependence of turbulence scattering on the wave number k, turbulence intensity , and transmission distance L. It should be noted that this formula is applicable to weak turbulence conditions and assumes turbulence in a uniform path. In the practical application of the Edwardsiella vaccine solution, it is also necessary to consider the influence of temperature gradient, concentration inhomogeneity, and flow characteristics inside the solution on the refractive index fluctuation. These factors together determine the effective turbulence intensity in the solution.

2.3 Light Transmission Equation

Comprehensively considering the absorption, particle scattering, and turbulence scattering effects, the transmission of light in the solution follows the extended form of the Beer-Lambert law:

Among them, α is the absorption coefficient, and s is the transmission path length. After considering the multiple scattering effect, the outgoing light intensity can be expressed as:

This expression uses the empirical correction factor to approximate the influence of multiple scattering.

2.4Extension of Time-Varying Light Transmission Model: Dynamic Light Intensity Attenuation in the Archimedes Screw System

Based on the parameters of the Archimedes screw system in the context of this study, the original model is extended, and the light transmission distance L is replaced by a time-varying function L(t) .

2.4.1. ​​Geometric and Kinematic Analysis of the Screw System

Core parameters:

Screw length:

Screw radius: R=0.5m

Pitch: P=0.5m

Rotation speed: 

Axial flow velocity:

​​It is approximately considered that light propagates along the axial direction of the screw, but affected by the spiral structure, the effective path length changes periodically with time. The actual path length L(t) is determined by the axial distance and spiral curvature together.

2.4.2 ​Mathematical Modeling of Time-Varying Transmission Distance L(t)​

The parametric equation of the spiral line of the screw (in the cylindrical coordinate system):

The differential of the actual path length of light in the spiral path:

After integration:

Calculate the constant part:

Therefore:

However, this model assumes that light always propagates along the spiral path. In practice, it is necessary to consider the mean free path of light in the fluid. A more reasonable simplification is:

Among them,  (average path), and (amplitude, caused by turbulence disturbance).

2.4.3. ​​Time-Varying Extension of the Light Intensity Attenuation Model

Original attenuation model:

After replacing L with L(t), the time-varying light intensity is:

​Update of the turbulence scattering coefficient :

3 Numerical Methods and Algorithm Design

3.1 Monte Carlo Method for Simulating Multiple Scattering

The Monte Carlo method simulates the transmission process of light in turbid media by tracking the random walk paths of a large number of photons. The fate of each photon is determined by probabilistic events such as absorption, scattering, and transmission distance.

The algorithm flow is as follows:

1. Photon initialization: Set the initial position, direction, and energy weight of the photon.
2. Step size calculation: Determine the photon free path according to the total attenuation coefficient , where is a random number uniformly distributed in (0, 1).
3. Absorption event: The photon weight is attenuated as .
4. Scattering event: Determine the new scattering direction by sampling the phase function.
5. Boundary processing: When the photon reaches the boundary of the solution, record the outgoing information.
6. Termination condition: Terminate the tracking when the photon weight is below the threshold or the photon escapes from the solution.
7. Result statistics: After repeating the tracking of a large number of photons, statistically analyze the outgoing light intensity, polarization state, and time broadening.

3.2 Fourth-Order Runge-Kutta Algorithm for Solving the Transmission Equation

For simple cases where multiple scattering is not considered, the fourth-order Runge-Kutta algorithm can be used to directly solve the light transmission equation:

The algorithm steps are as follows:

1. Set the step size h and the initial value
2. Calculate the iteration point
3. Calculate
4. Calculate
5. Calculate
6. Update
7. Repeat until s=L
8. Among them is the derivative function.

3.3 Parameter Setting and Grid Design

To accurately simulate the light transmission characteristics, it is necessary to reasonably set the numerical calculation parameters:

• Geometric parameters: The solution thickness , and logarithmic grid sampling is used.
• Light intensity parameters: The incident light intensity , and logarithmic sampling is used.
• Optical parameters: The absorption coefficient , and the scattering cross-section .
• Turbulence parameters: The refractive index structure constant .
• Particle parameters: The concentration , and the radius .

The grid generation adopts logarithmic equidistant division to ensure the stability and accuracy of numerical calculation:

4 Simulation Results and Visualization Analysis

Table of Key Parameter Simulation Results:

4.1 Analysis of Light Intensity Attenuation Characteristics

By numerically solving the light transmission equation, we obtained the variation law of the outgoing light intensity with the transmission distance and incident light intensity. Figure 1 shows a heat map of light intensity attenuation in a double logarithmic coordinate system, where the color represents the logarithm of the outgoing light intensity.

Figure 1: Heat Map of Outgoing Light Intensity Attenuation with Propagation Distance and Incident Light Intensity

The results show that with the increase of the transmission distance, the light intensity shows an exponential attenuation trend, which is consistent with the theoretical expectation.

The contour analysis shows that under a fixed transmission distance, the outgoing light intensity has a linear relationship with the incident light intensity. Under the condition of fixed incident light intensity, the outgoing light intensity decreases exponentially with the increase of the transmission distance. The specially marked 50 W/m² contour identifies the effective working boundary of the light transmission system. When the incident light intensity is higher than this intensity, the Edwardsiella vaccine will produce an effective immune response in the fish body.

4.2Analysis of Dynamic Light Intensity Attenuation Characteristics

Figure 2: Heat Map of Dynamic Light Intensity Attenuation

This heat map reveals the spatiotemporal dynamic characteristics of light transmission in the Edwardsiella vaccine solution in the Archimedes screw system. The horizontal axis (time) shows that the screw rotation has a period of 45 seconds. The vertical axis (incident light intensity, 0.1-100,000 W/m²) is expressed in a logarithmic coordinate system. The color mapping reflects the outgoing light intensity after attenuation by absorption, particle scattering, and turbulence scattering of the vaccine solution. The clear inclined stripe pattern in the figure indicates that to achieve a specific therapeutic effect (such as the 50 W/m² threshold), the required incident light intensity changes periodically with the position of the screw.

Figure 3: Time-Varying Curves of Transmission Distance and Light Intensity

The time-varying diagram of transmission distance intuitively shows the optical path length modulation caused by the geometry of the Archimedes screw. It can be seen that the path length fluctuation originates from the optical effect of the spiral geometry. The minimum value corresponds to the approximately straight path, and the maximum value corresponds to the fully developed spiral path. The time-varying diagram of light intensity clearly shows the dynamic response of the outgoing light intensity under different incident light intensity levels, which is synchronized with the change of transmission distance. The high incident light intensity (1000 W/m²) produces an effective therapeutic window of 25-100 W/m², while the low incident light intensity (1 W/m²) is almost completely attenuated. It is worth noting that the modulation depth of all curves remains consistent (about 75%), indicating that the path length change rate is the dominant factor, independent of the incident light intensity.

Figure 4: Coupling Heat Map of Turbulence Intensity and Light Intensity

Based on the mathematical model constructed above, the above figure represents the relationship between the outgoing light intensity, the fluid turbulence degree, and the propagation distance when the fixed incident light intensity is = 1000 W/m². It can be seen from the above figure that there is a nonlinear coupling effect between turbulence intensity and light intensity: there is a strong nonlinear coupling relationship between the turbulence intensity and the transmission distance L. When the transmission exceeds a certain distance (about L = 0.1 m), even weak turbulence will lead to significant light intensity attenuation. In addition, there is a critical threshold phenomenon: there is a critical turbulence intensity threshold (about ). Beyond this threshold, turbulence scattering becomes the dominant mechanism of light attenuation. In the context of this study, when Edwardsiella flows in the Archimedes screw with weak turbulence, if exceeds , the dominant mechanism of light attenuation is turbulence scattering; otherwise, the dominant mechanisms are still absorption and particle scattering.

5.Sensitivity Analysis

To evaluate the influence degree of each parameter in the light transmission model on the output results, the model adopts the Sobol global sensitivity analysis method. This method can comprehensively evaluate the first-order effects, high-order interaction effects, and total effects of parameters, and is suitable for the sensitivity analysis of the nonlinear mathematical model of the weak turbulence fluid system in the Archimedes screw immersion device under the background of this study.

5.1 Range of Analysis Parameters

Sample setting strategy: Sobol sequence sampling is used, and a total of 8192 sample points are generated to ensure uniform coverage of the parameter space.

5.2 Sensitivity Analysis Results

5.2.1 First-Order and Total Effect Analysis

Figure 5: Comparison Diagram of First-Order Sobol Index and Total Effect Index

Figure 1 shows the comparison results of the first-order Sobol index and the total effect index of each parameter. It can be seen from the analysis results that:

The attenuation coefficient α has the highest sensitivity, with a first-order index of 0.4133 and a total effect index of 0.7613, indicating that α not only has a significant impact on light intensity transmission itself but also has a strong interaction with other parameters.

The sensitivity of the particle scattering coefficient is the second highest, with a first-order index of 0.2078 and a total effect index of 0.5023, indicating that it also plays an important role in the model.

The influence of the turbulence scattering coefficient is relatively small, with a first-order index of only 0.0267 and a total effect index of 0.1104, indicating that its direct impact on light transmission is limited.

5.2.2 Parameter Interaction Effect Analysis

Figure 6: Heat Map of Second-Order Interaction Effects

This heat map of second-order interaction effects reveals the coupling effect between parameters: there is a significant interaction effect (0.128) between α and , which indicates that the attenuation process and the particle scattering process are coupled with each other and jointly affect the light transmission characteristics. The interaction effects between other parameters are relatively weak, all below 0.05, indicating that the influence of these parameters on the model output is mainly through independent effects.

5.3 Discussion of Results

It can be seen from the result diagram that in the light transmission model, the attenuation coefficient α is the most critical influencing factor, with a total effect index of 0.7613, which dominates. Although the first-order effect of the particle scattering coefficient is only 0.2078, its total effect is increased to 0.5023 through the interaction with other parameters, showing a strong influence potential. In contrast, the influence of the turbulence scattering coefficient is the most limited, with both direct and interaction effects being relatively small, indicating that there is an obvious synergistic effect between α and .

Conclusion

This model establishes a comprehensive model that can accurately describe the transmission of light in dynamic and turbid fluids, covering multiple physical processes from micro-particle scattering to macro-turbulence effects. From the simulation results, the absorption coefficient α of the solution should be prioritized for control and accurate measurement, as it has the greatest impact on the final light intensity. For process control, the turbulence intensity of the solution should be controlled below the critical threshold, which can prevent turbulence scattering from becoming an unnecessary attenuation source and improve the light transmission efficiency.

The code of index fluctuation Model are available in our GitLab Software Tool repository. The requirements are as follows:

Fish Vaccine Diffusion Mass Transfer Model

1.Overview

In aquaculture, vaccination serves as a vital measure to prevent disease transmission. However, traditional vaccination methods suffer from inefficiency and challenges in precise control. To address these issues, we developed a diffusion-transport model for vaccine administration, aiming to investigate the transport process of vaccines within fish bodies through mathematical modeling and numerical simulations.

The model takes into account the following key steps:

1. Diffusion of water body to mucous layer: The vaccine diffuses from the water body to the mucous layer on the surface of the fish.
2. Mass transfer between mucus layer and gill mucosa: The vaccine is transmitted to the gill mucosa through the mucus layer.
3. Internal diffusion of the gill mucosa: the process of vaccine diffusion in the gill mucosa.
4. Transport in the body circulation: The vaccine enters the blood from the gill mucosa and is distributed in the body circulation.
5. Mass transfer from blood to intestinal mucosa: The vaccine is transferred from blood to intestinal mucosa.
6. Internal intestinal mucosal spread: The vaccine spreads in the intestinal mucosa.
7. Colon colonization: The vaccine is ultimately colonized in the intestinal lumen.

In this model, the finite difference method was used for numerical solution, the implicit format was used for long-term simulation analysis, and the Sobol index was used for parameter sensitivity analysis to evaluate the influence of each parameter on vaccine effect, so as to provide a scientific basis for the optimization of vaccination strategy.

Figure 1. Flow chart of bacterial infection

Figure 2. Schematic diagram of bacterial infection

2. Main models and interpretations

Step 1: Water diffusion into the mucus layer

When examining the immersion process of fish in vaccine solution, we observe a thin mucilage layer on their bodies. Consequently, the vaccine spreads into this mucilage layer first as the water flow propagates. Given the complex composition of the mucilage, we simplify its properties by assuming it matches the water composition. Additionally, since the water flow volume is significantly larger than the fish's body volume, the disturbance to water remains minimal during immersion. This results in weak convective effects, allowing us to disregard the convective term when applying Fick's Second Law. The simplified equation can be expressed as:

Among them, cwi refers to the vaccine concentration in the mucus layer (CFU/m³), Dw refers to the diffusion coefficient of the vaccine in the main body of the water flow (m²/s), δ refers to the thickness of the mucus layer (m), and t refers to the time (s).

Step 2: Mucosal boundary and mucosal boundary mass transfer

After the vaccine spreads to the boundary of the mucosal layer, the classical two-film theory is considered to simulate the vaccine transfer between the mucosal layer and the gill mucosa. Using the interfacial mass transfer formula, it can be obtained that:

Among them, Nw→m refers to the rate of vaccine transmission from mucous side to mucosal side (CFU/s·m2), kw→m refers to the total mass transfer coefficient of bacteria in mucous side and mucosal side (m/s), cwi refers to the concentration of mucous layer (CFU/m3), and cm refers to the concentration of mucosal layer (CFU/m3).

Step 3: Internal mucosal self-diffusion

The chemical composition and structure of mucous membranes in organisms are highly complex. To better describe vaccine diffusion within mucosal tissues, we propose the following assumption: The mucosa acts as a homogeneous medium where vaccines can diffuse freely. Considering its biological effects, we incorporate first-order elimination reactions on the right side of the diffusion equation, yielding:

Among them, cm refers to the vaccine concentration on the mucosal side (CFU/m3), Dm refers to the diffusion coefficient of the vaccine on the mucosal side (m2/s), km is the elimination coefficient inside the mucosa (S-1), and Lm refers to the mucosal thickness (m).

Step 4: Systemic circulation compartmentalization balance

When the vaccine spreads within the mucosal lining, it gradually enters the bloodstream. At this stage, the vaccine traverses the submucosal microvessels. Given their extremely fine and thin nature, their mass transfer effects can be disregarded. Consequently, the model is further simplified to direct contact between the mucosa and blood, resulting in mucosal-blood mutual mass transfer. The equations can be expressed as:

Among them, Nm→b refers to the rate of vaccine transmission from mucosal side to blood side (CFU/s·m2), km→b refers to the total mass transfer coefficient between mucosal side and blood side (m/s), cm refers to mucosal layer concentration (CFU/m3), and cb refers to blood layer concentration (CFU/m3).

In fish, the blood volume is relatively small compared to the mucosal surface area. This allows for instantaneous homogenization and equilibrium of vaccine distribution in the bloodstream without prolonged diffusion mixing. Based on this mechanism, we can further simplify the model: The blood concentration of vaccine components originates from direct diffusion through the skin and gill mucosa, while diffusion loss occurs at the intestinal mucosa and natural clearance within the bloodstream. The simplified equation can be expressed as:

Among them, Vb refers to the blood volume of fish (m3), cb refers to the vaccine concentration in blood (CFU/m3), kb refers to the clearance coefficient in blood (S-1), Jb→g refers to the flux of vaccine from blood to intestinal mucosa (CFU/s), and Jm→b refers to the flux of vaccine from mucosa to blood (CFU/s).

Step 5: Blood-intestinal mucosal mass transfer

When the vaccine is in equilibrium with the blood, the vaccine will undergo blood-intestinal mucosal interaction mass transfer, and it can be obtained:

Among them, Nb→g vaccine transmission rate from blood side to intestinal mucosal side (CFU/s·m2), kb→g refers to the total mass transfer coefficient of vaccine from blood side to intestinal mucosal side (m/s), cb refers to the concentration of blood side, and cg intestinal mucosal side concentration (CFU/m3).

Using (4) and (6), the expressions for flux Jm→b and Jb→g are obtained:

Here, cm(Lm, t) represents the concentration (CFU/m³) in the mucosal lining at time t, cb(t) denotes the vaccine concentration in blood at time t (CFU/m³), Am→b indicates the contact area between mucosa and blood (m²), Ab→g refers to the contact area between blood and intestinal mucosa (m²), and cg(0, t) stands for the vaccine concentration (CFU/m³) on the outer surface of intestinal mucosa at time t.

Step 6: Intestinal mucosal spread

In the mucosa of fish intestines, we also simplify it as a homogeneous medium in which the vaccine diffuses freely. Similar to the gill mucosa, this area also has biological effects. By introducing a first-order reaction elimination term, we can obtain:

Among them, CG refers to mucosal side concentration (CFU/m3), Dg refers to intestinal mucosal diffusion coefficient of vaccine (m2/s), kg refers to intestinal mucosal clearance coefficient (S-1), and Lg refers to intestinal mucosal thickness (m).

Step 7: Gut Balance

After the vaccine passes through the intestinal mucosal barrier and reaches the intestinal lumen, it can play a role. Here we also assume that the intestinal lumen can reach instantaneous equilibrium, so the following expression can be obtained by using the vaccine flux:

Among them, the expression of flux Jg→in is as follows:

Among them, kin intestinal lumen clearance coefficient (S-1), kg→in total mass transfer coefficient of the vaccine on the intestinal mucosal side and intestinal lumen side (m/s), Jg→in refers to the flux of the vaccine from the blood into the intestinal mucosa (CFU/s·m2), Ag→in refers to the contact area between the intestinal mucosa and intestinal lumen (m2), Vin refers to the intestinal lumen volume of the fish (m3).

In conclusion, we get the system of equations

initial condition ：

Since the fish body has not been exposed to the vaccine immersion solution at the initial inoculation, the vaccine concentration is zero everywhere in the fish body. In addition, since the water volume is much larger than the fish body volume, the vaccine concentration in the water remains stable for a long time, so the initial conditions (13) ~ (18) can be obtained.

boundary condition ：

First of all, based on dimensional analysis, the relationship between mass transfer rate and concentration change can be obtained:

Using Equation (19), we can obtain the mucosal-mucosal boundary (Equation (21)), mucosal-blood boundary (Equation (22)), blood-intestinal mucosal boundary (Equation (23)), and intestinal mucosal-lumen boundary (Equation (24)).

3. Numerical simulation and stability analysis

In order to solve the model's partial differential equation, we consider using the finite difference method for solving. Based on the characteristics of the model, we decide to use the implicit difference for numerical solution.

Take the internal diffusion of the mucosa as an example (see Equation (3) for details):

First, we discretize cm and use the symbol umn, where the superscript refers to time and the subscript refers to displacement.

The second step is to replace the differential form with a difference. In time, we adopt the Euler front format and get:

Next, the space is trapezoidal, so we get:

Since (24) is equal to (25), we can get:

Simplifying , we have

Furthermore, we use the Von Neumann stability analysis to derive its growth factor G(k), which helps us to classify the formats.

Introducing Fourier mode errors

Substituting (28) into (27), we get:

Divide both sides by, and simplify to get:

For the growth factor , if the amplitude ratio between two consecutive oscillations is not more than 1 to maintain the stability of numerical simulation, it can be deduced that:

Using Euler's formula to simplify equation (30), we can get:

Using trigonometric functions, we can further simplify it as follows:

It can be deduced from Equation (31)

Similarly, we can treat Equation (1) and Equation (9) to get:

In this case, Equation (35) is always true, and Equation (34) and Equation (36) are conditionally true.

After completing the construction and numerical simulation of the vaccine diffusion model, we conducted in-depth analysis of how various parameters influence the model's output results. To clarify which parameters play a dominant role in vaccine diffusion processes, we performed parameter sensitivity analysis. This research provides crucial insights for subsequent experimental design and model optimization.

Based on this, we selected the area under the blood administration curve (AUC) as the evaluation index, that is:

Furthermore, we use Sobol sensitivity analysis, that is, using different output parameters to obtain the function value for variance decomposition to evaluate the significance of each parameter, and introduce the main effect sensitivity analysis index (Si) and total effect sensitivity index STi.

among ：

So we define the output variable matrix Y as:

Introducing the parameter input matrix A (N×11), B (N×11) can be obtained

By successively replacing the same columns of matrix A with columns of matrix B, 11 mixed matrices can be obtained, which are denoted as (i represents the replacement of the i-th column of A by the i-th column of matrix B).

For example:

Furthermore, we bring the parameter matrix into the model for calculation, and we get:

In the same way, 11 mixed matrices can be calculated:

To calculate Si and STi, the calculation formula is introduced as follows:

Among them, it is the variance of the expected value of the output Y with other variables X~i under fixed Xi. It refers to the average influence of the input variable Xi on the output variable Y considering the uncertainty of all other input variables.

4. Operation parameter table and result analysis

Figure 3. Time variation of vaccine concentration in fish

Figure 3(A) illustrates the vaccine diffusion process from the aqueous phase to the mucosal layer as described in Model (1). Simulation results show that vaccine concentration initially rises rapidly, then declines sharply due to mass transfer between the mucosal layer and gill mucosa (Model (2)). This trend aligns with the dynamic behavior predicted by the model. Figure 3(B) demonstrates that blood concentration first surges before gradually decreasing, consistent with the mucosal-blood interfacial mass transfer and clearance effects described in Models (4) and (5). The peak concentration may relate to transport rates through submucosal microvessels, while subsequent decline likely stems from blood clearance coefficient (kb). Figure 3(C) reveals that mass transfer from blood to intestinal mucosa and internal diffusion drive concentration elevation under Models (6) and (9). The observed concentration peaks and subsequent decline match model predictions, illustrating the distribution and elimination processes of vaccines in intestinal mucosa. Figure 3(D) shows concentration trends consistent with Models (10) to (12). After crossing the intestinal mucosal barrier, vaccines achieve transient equilibrium in the intestinal lumen, followed by gradual concentration reduction due to intestinal clearance coefficient (kin), reflecting their dynamic behavior within the intestinal cavity.

Figure 4. Sobol index convergence diagram

Analysis of the graph reveals that k_bg (blood-intestinal mucosa mass transfer coefficient), Dg (intestinal mucosa mass transfer coefficient), k_g_del (intestinal mucosa elimination coefficient), k_in_del (intestinal lumen elimination coefficient), and k_g_in (intestinal mucosa-lumen mass transfer coefficient) show negligible impact on AUC (area under the curve of blood vaccine concentration) within minor fluctuation ranges. Six parameters exert significant influence on AUC: Dw (main water flow mass transfer coefficient), k_wm (mucus layer-mucosal mass transfer coefficient), k_m_del (mucosal elimination coefficient), Dm (mucosal mass transfer coefficient), k_mb (mucosal-blood mass transfer coefficient), and k_b_del (blood elimination coefficient). Notably, AUC demonstrates highest sensitivity to changes in k_mb (mucosal-blood mass transfer coefficient) and k_b_del (blood elimination coefficient). Additionally, the Sobol index converges after approximately 900 simulation runs, indicating that the optimal number of simulations should be N> 900.

Figure 5. Sensitivity analysis of vaccination model

The figure reveals that AUC demonstrates the highest sensitivity to blood elimination coefficient, followed by parameters such as the mucosal-blood mass transfer coefficient, Dw (mass transfer coefficient in the fluid phase), k_wm (mucosal-luminal mass transfer coefficient), k_m_del (mucosal elimination coefficient), Dm (mucosal mass transfer coefficient), k_mb (mucosal-blood mass transfer coefficient), and k_b_del (blood elimination coefficient). Other parameters show slightly lower sensitivity. Notably, k_bg (blood-intestinal mucosa mass transfer coefficient), Dg (intestinal mucosa mass transfer coefficient), k_g_del (intestinal mucosa elimination coefficient), k_in_del (lumen elimination coefficient), and k_g_in (mucosal-lumen mass transfer coefficient) remain unaffected by minor fluctuations in AUC (area under the curve of vaccine concentration in blood).

Dynamical Analysis of Multi-Light-Regulated Biological Systems

1. Proteomic analysis of proteins involved in directional movement of Escherichia coli under green and blue light regulation

Modeling idea:

Step 1: Input of optical signals: Firstly, CcaS (histidine kinase) is activated. It should be noted that red light and green light are antagonistic to each other. In addition, we believe that during the reaction, there is a certain rate at which CcaS protein can be converted into a non-activated state, thus obtaining the following equation:

The second step is the phosphorylation of CcaR. At this time, we introduce the literature and relevant formulas to obtain phosphorylation and dephosphorylation:

among ：

Step 3 phosphorylated CcaR promotes cheA transcription activation

Step 4: CheA catalyzes the conversion of CheY to CheY_P. Similarly, we also use the formula in this paper:

among ：

Figure 1. Phosphorylation and dephosphorylation transformation diagram

Based on the above model, we carry out numerical simulation to assist experimental research. The simulated values of each parameter are as follows

First, we conducted numerical simulations using the aforementioned parameters to obtain Figure 1.2. The results showed that CcaS and CcaS* reached steady states within an extremely short time (less than 41 seconds), consistent with our literature review findings that light-sensitive proteins undergo conformational changes rapidly under illumination to perform their functions. Other key components (CheA and CheY_P) gradually stabilized over approximately 500 seconds of exposure, aligning with our expectations. Western blot analysis revealed increased levels of these critical components after prolonged illumination, while repeated wet-lab experiments confirmed minimal pre-and post-exposure variations in their concentrations (<5%), which perfectly matches our simulation predictions.

Figure 1.2 Concentration variation curve of each component over time ()

Furthermore, we attempted to analyze the specific effects of light intensity on protein CcaS and CcaS* using the aforementioned parameter values. It should be noted that we have not found direct literature specifying their interaction mechanisms. Our approach was based on a genetic design concept, where we hypothesized antagonism between green and red light (already reflected in previous equations) for numerical simulation. Through computer simulations capturing protein concentration data after 1500 seconds of single-light exposure (Figure 1.3), it is evident (as shown in the upper two panels of Figure 1.3) that minimal response occurred at extremely low light intensities (maintaining original concentrations). At 25lx light intensity (applicable to both green and red light), the photosensitive protein concentration remained constant. This indicates that exceeding 25lx at 1500 seconds achieves equilibrium in the photosensitive protein, representing the saturation light intensity.

To evaluate the appropriateness of our light capture timing parameters, we conducted steady-state duration monitoring across multiple illumination conditions. As shown in the two subplots of Figure 1.3, green light at 50lx achieves rapid protein activation within <40 seconds, while red light demonstrates a significantly delayed response, reaching approximately 400 seconds at 100 lx. It should be noted that our simulations utilize estimated parameters which may deviate substantially from real-world conditions. Nevertheless, this study establishes a framework and methodology for systematic analysis.

Figure 1.4 Response surface analysis of green and red light on S/S* concentration

Through comprehensive analysis of green and red light, we developed Figure 1.4. The figure clearly demonstrates a distinct plateau phase in their interaction, indicating that even with sufficient green light intensity, saturation can still be achieved under red light interference (acquisition duration: 1500 seconds). We plan to validate this simulation result through experimental verification in subsequent studies.

Finally, to identify which parameters most significantly influence the model's output, we employed Sobol's global sensitivity analysis (for detailed methodology, see "Model Development for Vaccine Free Diffusion into Fish"). The core approach involves breaking down the variance in outputs into smaller components: individual parameter contributions (Sᵢ) and their combined interaction effects with other parameters (STᵢ). Both Sᵢ and STᵢ range between 0 and 1, where higher values indicate greater parameter significance.

The operation steps are as follows:

1. Select 6 potentially sensitive parameters: Vmax_g, Vmax_r, Vmax_A, β_A, V1_Y, V1star_Y, each of which varies within the range of 20% of the baseline value.
2. 1000 parameter combinations are generated by Latin hypercube sampling to ensure uniform coverage.
3. Run the model for each group of parameters for 2000s, and take the steady-state concentration of phosphorylated CheY (CheY_P) as the output.
4. Calculate the first-order index Sᵢ and the total order index STᵢ according to the Saltelli method; the total order index can find the "hidden" interaction effect.
5. In order to reduce the random error, the sample size is gradually increased from 200 to 1000, and stopped when STᵢ changes by less than 1%, and finally it is determined that 1000 samples can meet the accuracy.

By comparing the size of STᵢ, the "high sensitivity" parameters can be quickly screened out to provide a basis for subsequent experimental determination or model simplification.

Figure 1.5 Sobol parameter sensitivity analysis

We found that the concentration of CheY_P had the most significant effect, which provided a direction for our subsequent experiments.

2.Protein component analysis of E. coli rectal motility under blue light regulation

Running logic:

Under blue light: YF1 is inactivated, FixJ can not be phosphorylated → inhibit pFixK2 promoter → decrease of cI protein → increase of CheZ expression → dephosphorylation of CheZ CheY-P → decrease of CheY-P → enhance of linear movement of cells.

Catalog of parameters

First, when exposed to blue light: The YF1 structure undergoes structural alteration and inactivation, losing its ability to phosphorylate the FixJ protein. At this stage, we believe that the rate of YF1 protein inactivation is influenced by two factors: the intensity of blue light and the concentration of YF1. It's important to note that for light-sensitive proteins, the response to light is extremely rapid, allowing us to formulate the following equation:

Next, for FixJ phosphorylation and dephosphorylation, we constructed our equations using steady-state enzyme kinetics formulas reported in the literature:

among ：

Next, we describe the relationship between phosphorylated FixJ and the repressor protein cI (using Hill equations to describe it):

Then the interaction between cI and chez:

Finally, for the transition between cheY and cheY_P:

among ：

Similarly, we carry out numerical simulation based on the above model to assist experimental research. The simulated values of each parameter are as follows:

Table 2.1 Parameter value estimation table

The analytical approach for this pathway is fundamentally consistent with the previous study. First, we conducted numerical simulations using the aforementioned parameters to obtain Figure 2.1. Our findings reveal that YF1 and YF1* reached steady states within an extremely short duration (<12s), which aligns with the magnitude of results from literature reviews—photosensitive proteins undergo conformational changes under light stimulation within brief periods, thereby altering their functions. Other key components (CheZ, CheY_P) achieved steady states over approximately 300 seconds after blue light activation, consistent with our expectations. Western blot analysis demonstrated increased levels of these critical components following prolonged illumination. Repeated wet-lab experiments confirmed minimal temporal variations in CheZ and CheY_P concentrations (less than 5%), which corroborates our simulation results.

Figure 2.1 Concentration curves of each component over time (dark for the first 500s, after 500s)

Figure 2.2 Light intensity response analysis

Furthermore, we attempted to analyze the specific effects of light intensity on proteins YF1 and YF1* using the aforementioned parameter values. It should be noted that we have not found direct literature specifying their interaction mechanisms. Through computer simulations capturing protein concentration data after 2000 seconds of blue light exposure (Figure 2.2), it is evident that both YF1 and YF1* exhibit extreme sensitivity to blue light. At an irradiance of IB=15 lx, they reach half-saturation light intensity. To evaluate the rationality of our illumination duration settings, we conducted steady-state time monitoring across multiple lighting conditions. As shown in the right panel of Figure 2.2, low light intensities (approximately 15 lx) can rapidly activate the proteins within <6 seconds. While it should be emphasized that our simulations utilize estimated parameters which may deviate significantly from real-world conditions, we have nevertheless established a framework and methodology for analytical purposes.

Finally, to identify the parameters with the most significant impact on the model's output, we conducted parameter sensitivity analysis. Seven key parameters (Vmax_B, IB_50, k_YF1, V1_F, V2_Y, Vmax_Fp_cI, Vmax_cI_Z) were selected for ±10% relative perturbations while keeping other parameters constant. After each disturbance, the ODE system was re-integrated during the blue light activation phase (500-1000 s), and the steady-state CheY-P concentration Yp at t = 1000s was extracted.

The sensitivity index is defined as follows:

among :

: CheY-P steady-state concentration under reference category;

: CheY-P steady-state concentration after increasing or decreasing the parameters by 10%;

Figure 1.5 Parameter sensitivity analysis (1)

Figure 1.6 Parameter sensitivity parameter ranking (2)

According to the simulation results, we found that the maximum reaction rate of CheY phosphorylase and the maximum response rate of blue light were the most sensitive to the biological system, which provided a preferred target for subsequent experimental verification.