Overview

This modeling study is divided into two main parts. The first part aims to identify the optimal exogenous invasion protein from both phylogenetic and molecular perspectives, with the goal of enhancing the cell invasion efficiency of engineered bacteria. The second part focuses on constructing a partial differential equation (PDE) model to investigate bacterial diffusion, advection, and growth within host cells.

In the first part, Integrated FnBP Screening, we focus on fibronectin-binding proteins (FnBPs), which are widely present across various bacterial species. By analyzing the phylogenetic relationships of FnBP-containing strains, predicting protein complex structures with AlphaFold-Multimer, and performing molecular dynamics (MD) simulations using GROMACS, we systematically assessed the stability and binding free energy of the complexes. The results indicate that FnBPA from Staphylococcus aureus is the most promising candidate as an exogenous invasion protein. When expressed in engineered bacteria, it is expected to significantly enhance their ability to invade tumor cells.

In the second part, Intracellular Bacterial Diffusion, we developed a coupled diffusion–advection–growth PDE model and performed numerical simulations using COMSOL Multiphysics. These simulations capture the spatial migration and population dynamics of bacteria within tumor cells, providing insight into their diffusion, advection, and growth behaviors.

Integrated FnBP Screening

Objective

To enhance the invasion efficiency of Staphylococcus epidermidis and Staphylococcus xylosus into tumor cells and improve delivery efficiency, we aim to introduce surface membrane proteins that facilitate bacterial entry. From the perspectives of phylogenetic analysis and molecular analysis, we comprehensively screened for the most suitable exogenous invasion proteins.

Introduction to FnBPs

Fibronectin-binding proteins (FnBPs) are widely distributed among Staphylococcus and Streptococcus species. They are localized on the bacterial surface and can bind to components of the host extracellular matrix (ECM), including fibronectin, fibrinogen, and collagen. FnBPs mediate bacterial invasion into host cells by recognizing fibronectin (Fn) and leveraging the integrin signaling pathway.

FnBPs typically contain multiple tandem repeat domains. For example, S. aureus FnBPA possesses 11 FnBP repeat units (FnBPRs), each capable of binding fibronectin (Fn) (Meenan et al., 2007). These repeat units are intrinsically disordered regions (IDRs) in their unbound state, but upon binding to the N-terminal F1 module of Fn, they form a tandem β-zipper structure, resulting in high-affinity and stable interactions (Speziale & Pietrocola, 2020).

Upon binding to Fn, the Fn molecule undergoes conformational changes that expose cryptic integrin-binding sites, effectively serving as a “bridge” connecting bacteria to host cell integrin α5β1. This process induces integrin clustering, activates downstream FAK and Src kinase pathways, and drives cytoskeletal rearrangements, ultimately facilitating bacterial invasion via a “zipper-like” mechanism (Kember et al., 2022).

In this study, the interaction between Fn and integrins is kept constant and is not treated as a variable. Our modeling focuses on the FnBPs anchored on the bacterial surface, evaluating the stability and affinity of the complexes they form upon binding Fn.

Phylogenetic Analysis

Objective

Bacteria that are phylogenetically closer to the engineered strains are more likely to successfully express and functionally accommodate foreign proteins. In this study, we constructed a phylogenetic tree encompassing common FnBP-producing bacteria as well as the engineered bacterial strains. By comparing the evolutionary relationships among these bacteria, we aimed to infer the structural and functional conservation of FnBPs and thereby assess the potential expressibility and functional stability of exogenous proteins in the engineered hosts.

Procedure

Through a systematic literature review, we compiled a list of fibronectin-binding proteins (FnBPs) from various commonly studied bacterial species (Henderson et al., 2011). The 16S rRNA sequences of the corresponding bacteria were aligned, and a phylogenetic tree was constructed using the Maximum Likelihood (ML) method in MEGA12 to resolve evolutionary relationships between the source strains and the engineered bacteria. On this basis, we further constructed a genetic distance heatmap, providing an intuitive and quantitative visualization of the phylogenetic relationships among the bacteria.

Bacterium	Fibronectin-binding protein
Staphylococcus aureus	FnBPA
Streptococcus pyogenes	SfbI
Streptococcus zooepidemicus	FnZ
Streptococcus dysgalactiae	FnBB
Streptococcus equi 4047	FnEB
Streptococcus pneumoniae	Pfba
Listeria monocytogenes	Lmo1829
Staphylococcus epidermidis（engineered strain）
Staphylococcus xylosus（engineered strain）

Results

The results showed that Staphylococcus species clustered together with Listeria monocytogenes, while Streptococcus species formed a distinct branch. Within the Staphylococcus genus, Staphylococcus aureus was phylogenetically closest to Staphylococcus epidermidis, with almost no genetic distance between them. In addition, S. aureus, S. epidermidis, and S. xylosus formed a tightly associated sub-branch with relatively small genetic divergence. Therefore, this close phylogenetic relationship indicates that the fibronectin-binding protein (FnBP) of S. aureus shares high structural and functional similarity with that of S. xylosus, the chassis strain used for engineering. This suggests a high likelihood of successfully introducing key FnBPA motifs from S. aureus into the engineered strain.

Molecular Analysis

Objective

To identify FnBP candidate proteins with the strongest binding affinity and highest research and application potential, we designed an integrated workflow based on structure prediction → molecular dynamics (MD) simulation → free energy calculation. This approach allows not only rational prediction of FnBP–Fn complex binding modes but also quantitative evaluation of complex stability and affinity from both dynamic and thermodynamic perspectives, providing a solid theoretical basis for downstream candidate selection.

Workflow

(1) AlphaFold-Multimer Structure Prediction

In the initial stage, the amino acid sequences of candidate FnBPs were used as input for AlphaFold-Multimer to predict the structures of their complexes with fibronectin (Fn) (Evans et al., 2021). Predicted models were filtered according to the following criteria:

• Overall complex conformation is reasonable without obvious steric clashes;
• Protein–protein interface is clearly defined, with identifiable amino acid interactions;
• Predictions are consistent with existing experimental structures or literature reports.

Filtered candidate complexes were then selected as starting points for subsequent MD simulations to ensure structural rationality and biological credibility.

(2) Molecular Dynamics (MD) Simulation

Based on the predicted structures, GROMACS was employed to perform MD simulations, capturing the temporal evolution of FnBP–Fn interactions (Abraham et al., 2015).

Force Field: CHARMM36, which provides high accuracy in describing protein backbone and side-chain conformations, as well as protein–protein interactions.

System Setup: Complexes were solvated in a cubic water box with a minimum distance of 1.0 nm from the protein to box edges. TIP3P water molecules were used, and 0.15 M NaCl was added to mimic physiological ionic strength and maintain charge neutrality.

Simulation Conditions: Temperature was set to 310 K and pressure to 1 bar, approximating human physiological conditions.

Simulation Protocol:

Energy minimization → NVT equilibration (constant volume) → NPT equilibration (constant pressure) → Production run.

A 500 ns pre-run was performed to monitor RMSD trends and energy convergence, ensuring the system reached equilibrium.Subsequently, a 1000 ns production simulation was conducted to fully sample conformational changes and interaction patterns at the nanosecond to microsecond timescale.

MD simulations enabled observation of the FnBP–Fn complex dynamic stability and residue-level interface interactions under near-physiological conditions.

(3) MM/PBSA Free Energy Calculation

To quantitatively evaluate the binding affinity of different complexes from a thermodynamic perspective, we applied the MM/PBSA (Molecular Mechanics/Poisson–Boltzmann Surface Area) method to the stable RMSD segments obtained from the MD simulations (Genheden & Ryde, 2015).

Binding Free Energy Formula:

$$\Delta G = G_{\mathrm{complex}} - \bigl(G_{\mathrm{protein}} + G_{\mathrm{ligand}}\bigr)$$

$$G = E_{\mathrm{MM}} + G_{\mathrm{sol}} - T\Delta S$$

Energy Decomposition:

Molecular mechanics energy \(E_{\mathrm{MM}}\): including van der Waals (\(E_{\mathrm{vdw}}\)) and electrostatic (\(E_{\mathrm{ele}}\)) interactions.

$$E_{\mathrm{MM}} = E_{\mathrm{vdw}} + E_{\mathrm{ele}}$$

Solvation energy \(G_{\mathrm{sol}}\): composed of polar Poisson–Boltzmann energy (\(E_{\mathrm{PB}}\)) and nonpolar solvation energy (\(E_{\mathrm{SA}}\)).

$$G_{\mathrm{sol}} = E_{\mathrm{PB}} + E_{\mathrm{SA}}$$

Entropic contribution: \(T\Delta S\), accounting for changes in conformational degrees of freedom upon binding.

By comparing the ΔG values of different complexes, we obtained the relative binding affinities of the candidate FnBPs to Fn. Further energy decomposition allowed identification of key amino acid residues contributing to binding, providing guidance for subsequent mutagenesis, experimental validation, and mechanistic studies.

Results

Stability analysis was performed on FnBP–Fn complexes from the 1000 ns MD simulations. To smooth transient fluctuations and highlight overall conformational trends, running average analyses were applied to the RMSD trajectories. This approach clearly visualized the stability, convergence, and long-term structural dynamics of the protein complexes.

Based on the RMSD trajectories smoothed using the running average, we selected the time segments where RMSD values were relatively stable and the complex docking conformations remained steady for MM/PBSA binding free energy calculations.

Complex	Time segment
Streptococcus pyogenes SfbI type1	0–800 ns
Streptococcus pyogenes SfbI type2	820–1000 ns
Streptococcus equi subsp. zooepidemicus FnZ type1	0–413.7 ns
Streptococcus equi subsp. zooepidemicus FnZ type1	779.1–1000 ns
Streptococcus dysgalactiae FnBB	650–750 ns
Streptococcus equi subsp. equi FnEB	410–522 ns
Staphylococcus aureus FnBPA type1	0–77.91 ns
Staphylococcus aureus FnBPA type2	526.31–937.94 ns

The results indicate that FnBPA binds most tightly to Fn, with \(\Delta G_{\mathrm{binding}}=-199.96 \pm 13.17~\mathrm{kcal\,mol^{-1}}\). The next strongest binders are SfbI and FnEB, with \(\Delta G_{\mathrm{binding}}=-187.84 \pm 11.33~\mathrm{kcal\,mol^{-1}}\) and \(\Delta G_{\mathrm{binding}}=-186.78 \pm 14.62~\mathrm{kcal\,mol^{-1}}\), respectively.

Given that the FnBB–Fn and FnZ–Fn complexes exhibited lower binding affinities, the MM/PBSA binding free energy results serve as a key criterion in the FnBP screening workflow. Consequently, subsequent molecular-level analyses focus on FnBPA, Sfbi, and FnEB.

To further dissect amino acid–level interactions, we generated a contact value heatmap for these complexes:

The analysis revealed that the contacts formed by FnBPA, FnEB, and Sfbi on Fn exhibit a regular distribution, with higher contact values observed in the repetitive sequence regions of Fn, consistent with the structural repeat units of Fn. FnI2–5 is composed of four tandem FnI domains, and upon binding to FnBPs, its disordered sequences can transition into an antiparallel β-sheet structure. Therefore, by performing contact value analysis of individual FnI domains with FnBPs and identifying their common features, we aim to uncover the binding patterns and underlying mechanisms.

The results indicate that FnBPs form significantly more contacts with FnI3, FnI4, and FnI5 than with FnI2, with interactions predominantly localized in the C-terminal regions of the domains, whereas the N-terminal regions exhibit relatively few contacts. This suggests that the N-terminal segments primarily function as linkers, while the C-terminal regions serve as functional sites more prone to binding with FnBPs.

However, a high number of contacts does not necessarily correlate with strong binding. Therefore, combined with MM/PBSA analysis, we further identified amino acid residues that contribute most significantly to the binding free energy.

By integrating the single FnI contact value heatmap with amino acid/MM-PBSA analysis, it was found that ARG46, ARG91, ARG128, ARG136, ARG134, ARG173 and ARG179 in Fn make significant contributions to the binding free energy. In FnBPA, Sfbi, and FnEB, the residues contributing most to the binding free energy are mainly ASP, ARG, LYS, and GLU, with ARG32 in Sfbi being particularly prominent. The following figure shows selected hotspot docking images. It can be observed that ARG173 in Fn forms bonds with ASP14 in FnBPA at distances of 1.6 Å and 2.5 Å, while GLU12 in FnBPA interacts with ARG179 in Fn at distances of 1.8 Å and 1.7 Å. Based on the contact values between Fn and FnBP, as indicated by the contact heatmap and the intrinsic structural characteristics of Fn reported in the literature, we can infer that these residues constitute key interaction hotspots characterized by charged hydrophilic side chains, which mediate hydrogen bonds and extended electrostatic interactions essential for stabilizing the Fn–FnBP complex.

Integrative Analysis of Phylogenetic and Molecular Results

At the phylogenetic level, the evolutionary tree shows that among bacteria harboring common FnBPs, the engineered strains Staphylococcus epidermidis and Staphylococcus xylosus are genetically closest to Staphylococcus aureus, all belonging to the same genus. Listeria monocytogenes is slightly more distant.

At the molecular level, FnBPA from S. aureus exhibits the tightest binding to Fn, with a ΔG_binding of −199.96 ± 13.17 kcal/mol, although the overall complex stability is moderate. In contrast, the Sfbi–Fn complex from Streptococcus pyogenes demonstrates higher stability, with smaller RMSD fluctuations and a ΔG_binding of −187.84 ± 11.33 kcal/mol.

Integrated analysis indicates that FnBPA is both phylogenetically closest to the target engineered strains and exhibits the lowest binding free energy in MM/PBSA calculations. Although its RMSD-based structural stability is not the highest, considering both evolutionary proximity and binding affinity, FnBPA emerges as the most suitable candidate for introduction into engineered bacteria as an exogenous invasion protein.

Reference

---

Intracellular Bacterial Diffusion

Objective

This study aims to simulate the motility, growth, and metabolic processes of bacteria within tumor tissue using a mathematical modeling approach. The focus is on analyzing bacterial migration patterns and population distribution in complex microenvironments. By integrating fluid dynamics with biochemical reaction mechanisms, the model provides insight into bacterial colonization within tumors and their therapeutic potential, offering a theoretical reference for future experimental work.

Mathematical Modeling

Assumptions

1. Continuum Assumption
   Bacteria and nutrients are treated as continuous concentration fields rather than discrete individuals.

2. Simplified Motion
   Bacterial motion is reduced to passive diffusion and advection:

   • Diffusion: Isotropic random diffusion driven by concentration gradients (Fick’s law).
   • Advection: Bacteria are treated as passive tracers transported by pre-defined fluid flows.

3. Simplified Growth Kinetics
   Growth depends only on a single substrate and bacterial density:

   • Nutrient dependence: Growth is limited by a single substrate (glucose) and neglects other nutrients, metabolic waste, pH, and oxygen.
   • Density inhibition: Population growth is constrained by local bacterial density.

4. Simplified Tumor Microenvironment
   The environment acts as a passive, static background. Tumor tissue complexity is simplified. Host immune responses (macrophages, neutrophils) are ignored.

5. Simplified Coupling
   Interactions are unidirectional. Feedback effects, such as bacterial proliferation altering tissue porosity or fluid properties, are neglected.

Model Parameters

Symbol	Biological Meaning	Value
\(D_b\)	Random motility coefficient of bacteria	\(1\times10^{-13}\ \mathrm{m^2\,s^{-1}}\)
\(\vec{u}\)	Convective transport velocity due to interstitial fluid flow	\(1\times10^{-8}\ \mathrm{m\,s^{-1}}\)
\(\mu_{\max}\)	Maximum bacterial growth rate under nutrient saturation	\(1\times10^{-4}\ \mathrm{s^{-1}}\)
\(K_m=K_s\)	Half-saturation constant for nutrient-dependent growth	\(1\times10^{-2}\ \mathrm{mol\,m^{-3}}\)
\(K_B\)	Bacterial carrying capacity	\(1\times10^{-2}\ \mathrm{mol\,m^{-3}}\)
\(Y=Y_1\)	Yield coefficient (substrate to biomass)	\(0.5\)
\(C_0\)	Initial substrate concentration within tumor cells	\(2\ \mathrm{mol\,m^{-3}}\)
\(B_0\)	Initial bacterial concentration at cell interface	\(3.17\times10^{-11}\ \mathrm{mol\,m^{-3}}\)
\(D_c\)	Substrate diffusion coefficient within tumor cells	\(5\times10^{-10}\ \mathrm{m^2\,s^{-1}}\)

The Mass Conservation Equation

Considering that the local concentration of bacteria within tumor tissue dynamically changes over time, accompanied by transport processes such as diffusion and convection, and that bacterial growth or death generates or consumes corresponding source terms, the model needs a mathematical framework capable of simultaneously describing temporal evolution, spatial migration, and reaction processes. The mass conservation equation precisely fulfills this requirement: through the local rate-of-change term \( \frac{\partial \mathrm{b}}{\partial t} \), the flux divergence term \( \nabla\!\cdot\!\mathbf{J}_{\mathrm{b}} \), and the source term \( R_{\mathrm{b}} \), it comprehensively characterizes bacterial diffusion, migration, and growth/death within the tissue. The same formalism is applied to the substrate, with the terms as described above.

The governing equations are:

$$\frac{\partial b}{\partial t}+\nabla\cdot\mathbf{J}_{b}=R_{b}$$

$$\frac{\partial c}{\partial t}+\nabla\cdot\mathbf{J}_{c}=R_{c}$$

Flux Term（Diffusion and Convection）

Considering that bacteria and substrate undergo random migration within the tumor cell, and their diffusion is influenced by local pore structures and viscosity, diffusion term \( -D_{b}\nabla b\) is introduced into the model to characterize the non-directional diffusion of bacteria and substrate driven by concentration gradients in space.

Since the interstitial fluid within tumor cell typically exhibits slow flow, this flow can passively carry bacteria, resulting in directional migration. Therefore, a convection term \( \vec{u}\,b\)is introduced to simulate the passive movement of bacteria caused by fluid percolation. A similar term in the substrate equation,\( \vec{u}\,c\)is used to describe the transport of nutrients by the flowing fluid.

The Flux (\(\mathbf{J}\)) is decomposed into diffusive and advective components:

Diffusion (Fick’s law):

$$\mathbf{J}_{\mathrm{diff},b}=-D_{b}\nabla b$$

$$\mathbf{J}_{\mathrm{diff},c}=-D_{c}\nabla c$$

Advection:

$$\mathbf{J}_{\mathrm{adv},b}=\vec{u}\,b$$

$$\mathbf{J}_{\mathrm{adv},c}=\vec{u}\,c$$

Source Terms (Growth and Substrate Consumption):

Considering that the bacterial growth rate is closely related to substrate concentration—particularly under the hypoxic and nutrient-deficient conditions typical of tumors— the Monod Kinetic Form \( \mu_{\max}\cdot\frac{c}{K_{m}+c} \) is adopted to characterize the nutrient-dependent growth behavior. This formulation naturally captures the decline in growth rate when the substrate becomes scarce.

Within tumor tissue, spatial constraints lead to both physical and resource-based inhibition when bacterial density becomes high. To reflect this realistic limitation, a Logistic Inhibition Term \( \left(1-\frac{b}{K_{B}}\right)\) is introduced into the model, ensuring that as bacterial concentration approaches the carrying capacity \( {K_{B}}\), the growth rate gradually slows and stabilizes.

Since the substrate not only provides energy for bacterial growth but is also continuously consumed during the process, a consumption term is introduced into the Substrate Equation: \( -\frac{1}{Y}\,\mu_{\max}\frac{c}{K_{s}+c}\,b\) where Y is the yield coefficient, representing the quantitative relationship between substrate consumption and bacterial biomass production. This term enables dynamic coupling between bacterial growth and substrate utilization.

$$R_{b}=\mu_{\max}\cdot\frac{c}{K_{m}+c}\cdot\left(1-\frac{b}{K_{B}}\right)\cdot b$$

$$R_{c}=-\frac{1}{Y}\,\mu(c)\,b=-\frac{1}{Y}\,\mu_{\max}\frac{c}{K_s+c}\,b$$

Final Coupled Equations

Bacteria:

$$\frac{\partial b}{\partial t}-\nabla\!\cdot\!(D_{b}\nabla b)+\vec{u}\!\cdot\!\nabla b =\mu_{\max}\cdot\frac{c}{K_{m}+c}\cdot\left(1-\frac{b}{K_{B}}\right)\cdot b$$

Substrate:

$$\frac{\partial c}{\partial t}-\nabla\!\cdot\!(D_{c}\nabla c)+\vec{u}\!\cdot\!\nabla c =-\frac{1}{Y}\,\mu_{\max}\frac{c}{K_{s}+c}\,b$$

This framework converts complex biological processes into a computationally tractable model with clear biophysical meaning.

COMSOL Simulation

To visualize bacterial diffusion and concentration distributions within tumor tissue, the model was implemented in COMSOL Multiphysics, which solves coupled multiphysics PDEs while accounting for realistic tumor geometry and parameter heterogeneity.

Setup

Transport of diluted species physics.

Simulation duration: 6 hours, consistent with co-culture experiments.

Tumor cells: 3D spheres of radius 0.005 cm.

Initial bacterial concentration at cell boundaries: 3.17e-11 mol/m³.

Results

Radial profiles:

It can be observed that, as the simulation time progresses, bacteria with an initial concentration of 3.17 × 10⁻¹¹ mol/m³ located at a radial distance of 0.005 cm diffuse toward the central region (radius = 0). With increasing simulation time, the bacterial concentration near the sphere center gradually increases. Once a certain cellular concentration is reached, the values stabilize. Overall, the trend is consistent with expectations.

2D cross-sectional distribution:

The heatmap on the yz-plane through the spheroid center shows clear bacterial growth and convective diffusion within the tumor region from 0 to 5 hours, evidenced by a steady concentration increase that aligns with the radial profile. The consistent patterns between the 5- and 6-hour maps suggest concentration stabilization has been reached.

Optimization

The experimental results show that at 6 hours, the intracellular bacterial concentration is still increasing, which is inconsistent with the COMSOL simulation results where the bacterial concentration has already begun to plateau. However, as the co-culture time extends, the experiments indicate that the bacterial concentration eventually stabilizes. Meanwhile, the stable concentration predicted by the COMSOL model is slightly higher than the experimental value. We speculate that this discrepancy may stem from the current model’s simplification—that bacterial growth is solely determined by a single substrate and the bacterial density—leading to an earlier predicted stabilization and a slightly higher steady-state concentration. To optimize the model, we plan in future work to introduce more complex growth-limiting mechanisms, including multi-substrate-dependent growth equations, and to account for environmental factors such as pH and oxygen that regulate bacterial growth and diffusion, thereby improving the model’s fidelity to actual biological processes.

Reference