In this section, we aim to model the diffusion of lactose within the intestine, the transport of lactose by LacY, and the metabolic process mediated by lactase, with the goal of identifying key parameters that limit our metabolic efficiency.

Prior to the analysis, and to facilitate discussion by avoiding certain extreme special cases, we make the following assumptions:

1、The micro-environment at the bacterial cell surface is spatially discontinuous with the bulk intestinal lumen.

2、Each micro-environment is well-mixed, allowing us to neglect spatial gradients within it.

3、We assume the transport rate of LacY is linearly dependent on the proton gradient during the transport process.

4、For lactase, we assume that the adjustment of enzymatic activity in response to changes in substrate concentration is instantaneous, and that a basal degradation pathway for intracellular lactose exists.

5、To better analyze the key parameters of the engineered bacteria's metabolism, we employ a mean-field approximation to ignore lactose uptake rate variations caused by population effects. That is, the population behavior is described by an average density N, and the change in micro-environment concentration is proportional to both the population density and the per-capita uptake rate.

This equation quantifies the instantaneous rate of change of lactose concentration in the micro-environment at the bacterial cell surface. The change is governed by two primary terms:

The first term, adhering to Fick's First Law, represents the net influx of lactose from the bulk intestinal lumen into the micro-environment via passive diffusion.

The second term represents the total consumption rate of lactose by the entire bacterial population, driven by both basal transport and proton motive force-driven symport.

This equation defines the dynamic balance of lactose concentration within the cytoplasm. The net accumulation rate equals the transmembrane uptake rate minus the intracellular consumption rate. The uptake rate is the sum of basal transport and active transport mediated by the LacY protein, which is dependent on the proton gradient. The consumption rate includes the hydrolysis reaction catalyzed by lactase BgaB (following Michaelis-Menten kinetics) and degradation caused by non-enzymatic or bypass metabolic pathways.

This equation describes the dynamic balance of the transmembrane proton gradient that drives secondary active transport. Its rate of change is governed by three key processes: the first term represents the gradient generation rate driven by energy derived from lactose catabolism, which is proportional to the hydrolytic flux and drives proton pumping; the second term represents gradient decay due to passive processes such as membrane proton leakage; the third term represents the gradient consumption rate caused by the LacY protein consuming the proton motive force during lactose-H⁺ symport activity.

This equation describes the exponential growth kinetics of the bacterial population driven by energy derived from lactose metabolism. The specific growth rate of the population is proportional to the hydrolytic flux of BgaB, with the proportionality constant being the macro yield coefficient that converts substrate flux into biomass synthesis.

During our efforts to model this process, we noted that the specific degradation pathways of the BreR transcription factor and the molecular routes for the uptake and metabolism of deoxycholic acid by bifidobacteria are not yet fully elucidated. Furthermore, constrained by experimental timelines and necessary approval procedures, we are currently unable to investigate these issues in depth. However, to proceed with an initial attempt at this modeling approach, we hereby make the following assumption for this section: the intracellular concentration of the BreR transcription factor, expressed under the control of the strong promoter pgap, is constant, assuming no other influencing factors beyond the deoxycholic acid concentration.

For BreR (R) and deoxycholic acid (D), we have:

R + D → RD

This equation describes the kinetic process by which the BreR transcription factor binds to deoxycholic acid to form a complex. The association rate is proportional to the concentrations of free BreR [R] and deoxycholic acid [D], while the dissociation rate is proportional to the concentration of the complex [RD]. The model also accounts for the degradation of the RD complex.

This equation refines the expression kinetics of the MucBP protein. It incorporates the basal inhibition of the MucBP gene by free BreR [R] using an inhibitory Hill function, and the activation of the MucBP gene by the BreR-deoxycholic acid complex [RD] using an activatory Hill function. The parameter α represents the maximum induction multiple of DCA.

Considering the role of colonization capacity in intestinal microbiota growth and its contribution to lactose intolerance resolution in Lacbutler, we adjusted the population density changes of Bifidobacterium based on the adhesion system: We assumed that f(M) is a monotonically increasing function with a range of [f_min, 1], which describes the contribution of adhesion capacity to the actual growth rate.

among f_min: is the minimum colonization efficiency. When M = 0, i.e., there are no adhesion proteins present, f(M) =f_min. This means that the bacterial cells can only be present in a very low proportion (f_min) Grow by using the energy obtained from their metabolism.

K_M: Is the half-maximum efficiency constant. It reflects the sensitivity of MucBP protein concentration to colonization efficiency. When M = K_MAt this time, the difference between the maximum and minimum values of the planting efficiency f (M) is half. The original metabolic growth term $K_{H} \frac{V_{L,m}}{K_{c}+L_{m}}$ represents the potential growth capacity of the bacteria under ideal colonization conditions. After multiplying by f(M), the actual growth capacity is obtained.

As before, we assume that:

1、Considering that the corresponding mechanism of AI-2 has only been demonstrated to exist in Bifidobacterium, specific transport and decomposition mechanisms require further investigation. Additionally, we employed a strong promoter to activate it, thereby minimizing interference from other bacterial groups adopting the AI-2 strategy. This approach allows us to describe related processes using a constant parameter.

2、Considering that the promoter plsrA driving ccdB expression has been reported to have a threshold property for AI-2, we used a step function to approximate the description of its process

This equation quantifies the dynamic equilibrium of intracellular AI-2 signaling molecules. The rate of change is determined by four processes: sustained synthesis, where the LuxS enzyme maintains a constant rate kluxGeneration of AI-2, active uptake, i.e. the Mis kinetics describes the uptake of AI-2 from the environment, and active discharge, i.e. the first-order kinetics rate kexportTransport to the extracellular and intracellular degradation is given by the rate constant dinphysiolysis.

The equation describes the population accumulation dynamics of AI-2 signals in the environment. The first term represents the total amount of AI-2 discharged by all bacterial cells, establishing a quantitative relationship between population density and environmental signals. The contribution of individual bacteria is proportional to their intracellular AI-2 concentration and discharge rate.

The second item describes the reduction of environmental signals, which includes two mechanisms: physical dilution, caused by intestinal peristalsis and fluid flow, and chemical dilution, caused by the natural decomposition rate of AI-2 molecules in the environment.

The equation describes the expression regulation of the CcdB toxic protein in a suicide system. It employs a continuous threshold response model, using a Sigmoid function to simulate the switching characteristics of the plsrA promoter. The maximum expression rate is denoted by k_c, which characterizes the synthesis intensity of CcdB. Additionally, the degradation of the CcdB protein is described by the term d_c⋅C, representing the first-order kinetic rate of degradation.

And that gives us the final description of population density:

By considering the lactose content in the gut and assuming that the increase of lactose mainly comes from diet and other consumption items are described as a first-order kinetic function related to lactose concentration, we can derive the differential equation about the lactose content in the gut

The lactose input term utilizes methods from digestive dynamics, employing an exponential decay model for the initial maximum rate to simulate the dynamic process of lactose entering the intestinal tract from external sources. The engineered bacterial uptake term quantifies the total lactose intake rate by Bifidobacterium populations, while the natural consumption term represents lactose degradation caused by non-bacterial factors. Compare this formula with

as well as

By combining, the theoretical description of lactose content changes in the intestine can be described.

Given the current limitations in understanding critical components planned for our application, many parameters can only be estimated through empirical methods. However, all parameters can be validated through experimental or in vivo testing. Therefore, we consider these findings as references for anticipated outcomes in wet-lab environments and as attempts to mathematically model such components. The methodological value of this approach should far outweigh its practical data interpretation.

It is evident that lactase BgaB plays a crucial role in lactose decomposition. Given the significant industrial and domestic importance of lactase, we plan to conduct directed evolution simulations targeting BgaB to identify lactase variants with enhanced enzymatic activity.

Clearly, we need to first understand the mechanisms of action between BgaB and lactose. Therefore, we plan to use molecular dynamics simulations to elucidate their specific interaction mechanisms. The primary task involves determining the three-dimensional structures of both the substrate ——lactose and the BgaB protein.

We obtained the molecular structure information of lactose (PubChem CID: 5913) from PubChem, a widely used chemical information database, and provided the molecular structure data. The data was converted into PDB format using Avogadro to facilitate future simulations. Additionally, for better visualization, we recommend using PyMOL, a molecular visualization tool, to examine the obtained molecular structures.

As for lactase BgaB, we did not find its structural data in the existing database. Therefore, we found the FASTA sequence of the lactase we planned to use on NCBI and used AlphaFold3, which is one of the most accurate protein structure prediction algorithms at present, to generate its 3D structure.

It can be noted that the structure we generated has ipTM = -pTM = 0.97, so we believe that it has high credibility.

One way we need to understand how BgaB interacts with lactose is through molecular dynamics simulations, and we will do that next.

When performing molecular dynamics simulations, a force field parameter file is essential to define atomic types and various interactions. We selected CHARMM36, available from the MacKerell Labs website, as a widely-used force field for modeling protein-small molecule interactions that has been rigorously tested through numerous experiments. Alternatively, other force fields provided by CGenFF could be considered. After preparing the force field parameter file, we utilized

To convert the BgaB structure into a GROMACS - readable .gro file, you can select the force - field parameters and water model during the process.

When generating CHARMM36 parameters for BgaB, we encountered issues. The BgaB structure from AlphaFold3 had significant differences in atomic type definitions compared to CHARMM36 parameters. Inconsistent atomic type definitions are unacceptable in simulations.

After discussion, we decided to use CHARMM36 force - field parameters uniformly. We manually modified the atomic definitions in the BgaB’s pdb file and created the force - field parameter file. After verification, it was used in subsequent steps.

For small molecules like lactose, file format conversion is also required to enable GROMACS calculations. Since the CHARMM36 force field cannot recognize lactose molecules, we need a method to generate parameters for lactose within this framework. Here, we opt to use CGenFF for parameter generation (incidentally, quantum mechanical methods could also be employed to derive consistent solutions for original force field parameters of this small molecule, though this might require some computational effort). Therefore, obtaining .mol2 files compatible with CGenFF becomes essential.

Since the CHARMM force field uses a full-atom approach, hydrogen atoms are explicitly represented. However, our lactose crystal structure lacks hydrogen coordinates, requiring the use of Avogadro software to add hydrogen atoms and generate .mol2 files. We observed that some programs fail to list bond indices in ascending order when creating .mol2 files, which may lead to issues when constructing accurate topology files with matching coordinates. We recommend using Justins sort_mol2_bonds.pl script to resolve this issue. Consider the following:

Next, we upload and generate the .str file about lactose on CGenFFs official website. It should be mentioned here that the obtained .str file contains CGenFFs score for charge distribution and dihedral angle parameters, which is generally recommended to check that it is less than 10, and manual assignment is recommended for those greater than 50

However, the generated .str file cannot be used directly by GROMACS, so we use the auxiliary python script provided by CHARMM36 to help us convert it into a format recognizable by GROMACS.

It should be noted here that this script requires the NetworkX 1.11 version package and is compatible with Network 2.Version 0 is not compatible. You can download it from the NetworkX website or use

After this command, you will get a .mol2 file with the same name as the original .mol2 file. Use this command to convert it into a .gro file

Next, we need to merge our two research objects ( the .gro file of BgaB and the lactose molecules) into a single gro file for subsequent operations. Given that GROMACS built-in commands make this task challenging, we recommend performing the manual merging process instead.

Create a new .gro file (here we call it complex.gro)

Copy Protein.gro (here BgaB.gro) to complex.gro.

Copy the coordinates from Compound.gro (this is the coordinate part in the (Lac.gro) file to the complex.gro file, and paste it between the long vectors of the protein nucleus box.

According to the total number of atoms in the Protein.groand Compound.grocomplex.grofile, the number of atoms is usually on the second line

In the previous gmx pdb2gmx command, we set up a topol.topfile for the protein, and we only need to

Add it after the "Protein" molecular type section in the topol.top file #include " compound.itp".

Given that the ligand molecule contains new dihedral angle parameters in the compound.prm file, it is necessary to insert a parameter at the top of the topol.top file.

#include " compound.prm"

Finally, we need to modify the [molecules]. We need to add the compound molecules to the topology file

Note that the inclusion of the #include statement is very strict and must be added before the [moleculetype] because, according to the calculation logic, all atomic types must be known before the bonding parameters.

Therefore, parameters must be included before constructing any molecule. They must also appear after the #include declaration that contains the main force field.

Next we use:

To complex.gro create a box for a regular dodecahedron, and set the distance between the box and the molecule

This command is used to help generate the newbox.gro to generate a solvent environment according to the TIP3P water model.

Heres the deal: If youre paying attention, you’ll notice the protein carries a -17 charge in the previous pdb2gmx output. For those who missed it in the pdb2gmx output, check the [atoms] section of the topol.itp file. Youll spot "qtot-17" right at the end of the last atom line. Lets tackle this issue head-on.

We first obtain a .tpr file using Grompp from the mdp file. For computational efficiency, we recommend using the energy-minimized mdp file (you can find details about its generation in the next section). Run:

Then pass the obtained .tpr file to genion:

The ion names defined by pname and-nname are standard names in GROMACS 4.5. Finally, your [molecules] section should look like this:

Now that the system is set up, we will use grompp to generate input information. The input file and mdp file need to be designed by yourself.Here we present our mdp file as a reference:

Then call mdrun for energy minimization (EM):

First of all, in order to constrain the ligand molecule, it is necessary to generate the position constraint topology file of the ligand.

Create an index group for the compound that contains all atoms except hydrogen: Then, execute the genrestr module and select the index group obtained earlier

The obtained information is then incorporated into the topology file. Depending on environmental conditions, there are several different methods for adding these components. Here we recommend separately constraining proteins and ligands, using different #ifdef statements to control the inclusion of the ligand constraint file. Here we provide an example:

In this case, to constrain both the protein and the ligand simultaneously, you need to specify define = -DPOSRES-DPOSRES_LIG in the mdp file.

To avoid the problem that the temperature coupling function lacks sufficient stability in controlling the kinetic fluctuations of a small atomic system, we consider the protein and complex, ions and solvents as two wholes for temperature coupling. Therefore, we need to create an index group

Then perform NVT balance, the .mdp file is here.

After NVT balance, perform NPT balance ,the corresponding .mdp file is here:

We performed a MD simulation here. Here we did a 5-ns MD simulation, and the corresponding. mdp file is here.

After searching the relevant literature, we found that the active catalytic site of BgaB was at Glu303 and Glu148, so we adopted

To calculate the distance between lactoses hydroxyl group and the carbonyl oxygen atom of Glu303, our obtained average data is 0.30±0.04 nm. After analyzing the amino acid residues around the catalytic site, we concluded that hydrogen bonding is the primary mechanism driving the catalysis. Through referencing relevant literature and conducting visualization analysis using pymol, we confirmed our findings:

It is evident that the Glu148 and Glu303 sites play primary catalytic roles, with additional contributions from Trp311, His354, Glu351, Arg109, and Asn147. Our molecular dynamics simulations revealed each residues average energy contribution, with key residues highlighted in red circles. These critical sites will be the focus of our subsequent analysis.

Now that we have a better understanding of the catalytic mechanism of BgaB for lactose, we aim to identify amino acid sites with high evolutionary variability that enhance the enzymes catalytic activity. Using Hotspot Wizard, we conducted an analysis to identify additional potential sites. Our research has newly identified cysteine (Cys) at position 113 and asparagine (Asn) at position 147 as key contributors.

After identifying mutation-prone sites, we observed an enormous number of potential combinations. To streamline our approach, we first conducted preliminary stability screening to eliminate combinations that negatively impacted BgaB stability. For stability analysis, we utilized PyroSetta, a tool that can only handle single-point mutations. You may reference our reference scripts, which significantly reduce workload while offering strong adaptability for other tasks.

After conducting a comprehensive analysis of the structure and referencing relevant literature, we identified BgaB as an industrial-grade lactase requiring high thermal stability. Our research revealed its capacity to produce galactooligosaccharides with significant probiotic benefits, though at relatively low yields. To enhance these advantages, we plan to engineer mutations at four key sites: Y272A, E351R, F341T, and R109W. For site 272, we propose converting the benzene ring and hydroxyl group to methyl groups to reduce the catalytic cavitys compatibility with galactose hydrolysis products, thereby decreasing galactose competition at the active site. At site 351, we aim to stabilize lactose binding by modifying its electrical properties. Regarding site 341, we observed that the benzene ring may exhibit stronger stabilizing effects on smaller molecules, suggesting replacement with a smaller amino acid could diminish its galactose-supporting capacity. Finally, for site 109, we consider introducing an indole ring to boost hydrophobic interactions and strengthen hydrogen bonding efficiency.

Similarly, we aim to conduct molecular dynamics simulations on our modified BgaB molecule again and analyze its properties in detail to validate the modification effects. Here, we will not repeat the molecular dynamics simulation process. Using the molecular structure obtained from AlphaFold3, we observe that the generated structure shows an ipTM = -pTM = 0.97 value, indicating high reliability.

We start with

And use xmgrace to visualize the results, but you can also use python, you can refer to this

We noticed that the fluctuation was between 4.8 and 5.0 nm. This indicates that the structure of the enzyme remained relatively stable during the simulation, without large conformational changes, and was relatively stable.

We further analyze the variation of the radius of gyration of the modified BgaB

The rotation radius can be seen to range from 2.6 to 2.7 nm, which reflects the flexibility and dynamic changes of the protein during the simulation. Smaller fluctuations indicate that the core structure of the protein is stable.

Finally, we tried to quantify this binding strength as an alternative way to evaluate our enzyme modification

We determined the binding energy to be-231.06 kJ/mol. By comparing this with the original BgaB-lactose system (binding energy of-212.85 kJ/mol), our modification shows significant improvement in binding stability and accuracy. However, the enhancement remains modest and unlikely to cause excessive product release, suggesting potential enzyme activity enhancement. Further validation across additional parameters is required to confirm these findings.

To evaluate our modification effects, we employed Auto Dock Vina for molecular docking. Since the software cannot directly process PDB files, we utilized AutoDockTools to convert lactoses PDB file into a PDBQT format. Given that proteins exist in complex conformations, removing water molecules from the original file significantly reduced computational complexity. Additionally, adding hydrogen atoms—a crucial factor in determining protein bioactivity—was deemed essential. This approach ultimately yielded the final PDBQT file.

When performing Vina docking, we need to create a config.txtfile manually to define input and output files, determine the binding center coordinates, and specify the docking box dimensions. You can either use pymol to calculate the binding center coordinates or rely on data from the PDB file for this purpose. The docking box size must be larger than your target molecules dimensions. After obtaining the parameters:

We are showing here the affinity scores for each combination of BgaB, lactose and galactose, Glu148 and Glu303 before and after the modification to comprehensively evaluate the effect of our modification

1、After modification, BGA B and lactose are at the Glu148 site

2、The original BgaB and lactose are at the Glu148 site

3、After modification, BGA B and lactose are at the Glu303 site

4、The pre-modified BgaB and lactose at the Glu303 site

Clearly, both the Glu303 and Glu148 site modifications in the modified BgaB strain demonstrated improved activity scores. Lower affinity generally correlates with enhanced catalytic efficiency, indicating our modifications have achieved notable progress. However, we aim to investigate its binding affinity with galactose, as this interaction may not only influence catalytic performance but also interact with galactooligosaccharides—a probiotic component that plays a significant role in gut health.

1、The pre-modified BgaB and galactose at the Glu303 site

3、The pre-modified BgaB and galactose at the Glu148 site

4、The modified BgaB and galactose are at the Glu148 site

It can be observed that the affinity of galactose at the Glu303 site increases, while the Glu148 site remains essentially unchanged. This suggests that galactose has a more difficult time binding to the lactose active site, which may make oligogalactose formation easier. This indicates that our Lacbutler exhibits stronger probiotic effects than conventional Bifidobacterium beyond its basic functions.

Next, we will try to measure the modified BgaB in the wet laboratory and detect its relevant data as the final verification.