Our Model: AQUAINT
Introduction
A foundational understanding of how synthetic biology can be integrated into fluid dynamic environments is critical for advancing both theoretical and applied aspects of the field. Biological systems are inherently complex, and when introduced into environments governed by fluid flow—such as pipes, open water, or engineered surfaces—their behavior becomes even more difficult to predict. To address this challenge, a robust and modular mathematical framework is required to abstract and simulate the interactions between synthetic constructs and their surrounding fluid environments. Such a framework enables researchers to assess whether a construct can function reliably and sustainably under specific physical conditions.
The challenge lies in the vast diversity of synthetic biology applications and the multitude of possible deployment scenarios. Each use case—whether targeting microbial communities, therapeutic delivery, or environmental remediation—presents unique constraints. Additionally, incorporating fluid dynamics into these models is essential for realism, but doing so while preserving the ecological and community-level dynamics affected by genetic engineering introduces significant complexity. The interplay between engineered function and environmental transport must be carefully balanced to avoid oversimplification or loss of biological fidelity.
To begin addressing these challenges, we developed AQUAINT, a modular framework of novel mathematical models tailored to distinct fluid dynamic contexts. These include scenarios such as phage transport through plumbing systems, bacterial dispersal in flowing aquatic systems while preserving community interactions, and the deployment of biofilms onto various steel surfaces. Each model captures key aspects of fluid transport, biological persistence, and interaction dynamics, offering a foundation for predictive design and deployment. These models allow us to identify the key design principles for genetic engineering and effective treatment, guide the development of constructs, and evaluate how well a given chassis can function within its specific aquatic environment. With this framework, we developed novel mathematical models to describe the deployment of our chassis in aquatic environments, further extracting information of how the chassis may behave.
Additionally we incorporate the usage of our machine learning classifier AQUIRE to determine how biotic, abiotic, and mechanical factors of systems affect our chassis’ functionality in a given system.
Ultimately, integrating fluid dynamic principles into synthetic biology modeling of aquatic systems is not only essential for scientific rigor but also for improving the efficiency, reliability, and scalability of applications concerned with synthetic biology.
Open Moving Water
Chassis Transport
For transport of a synthetic biology chassis in running water, the transport can be described by a partial differential equation governing its time evolution as it is transported through advection and diffusion within a fluid environment. Given the computational complexity of computational fluid dynamics simulations using Navier-Stokes equations, a good initial test of feasibility begins with partial differential equations representing advection, diffusion, and reaction through the transport equation incorporated into the various dynamics.
In the above PDE, C represents the concentration of the chassis, and
where u is the averaged flow velocity, L is the characteristic length in the direction of the flow velocity, and D is the molecular diffusion coefficient of the chassis. When the Peclet number is extremely high, it indicates that the model can be approximated by ignoring diffusion, easing simulations.
The various other elements important to the system can be modelled similarly. For example, if S is a concentration dynamically related to C:
To initialize the system, we represented both the chassis and the target population using Gaussian spatial distributions. Each population was centered at distinct positions along the domain to reflect an initial separation between the engineered construct deployment and its intended target. The variance (σ) of each Gaussian was chosen to capture differences in spatial extent: the target population was assigned a broader distribution to represent its dominance in the environment, while the chassis population was initialized with a narrower distribution to reflect localized deployment. Here xc denotes the deployment location of the chassis in space and xt denotes the centered point in space of the target concentration which the chassis will affect:
This transport framework applies broadly to any synthetic chassis deployment where it is being transported in a moving water system.
By explicitly incorporating advection and diffusion into chassis design, synthetic biologists can better predict how constructs disperse, persist, and interact with surfaces or communities. Furthermore, they can identify engineering targets that improve performance in specific environments, such as adhereability of the construct. Through this framework, synthetic biologists can predictably optimize deployment strategies by matching chassis properties to flow regimes, geometries, and system contexts. Additionally, they can forecast chassis performance to determine whether their chassis can carry out its desired function. The advection-reaction-diffusion approach where u represents a streamline of the flow is a good starting point in the construction of models that incorporate fluid dynamics through computational simulation of Navier-Stokes equations for viscous flow.
Executive Summary
Mathematical modeling with parameter values representative of cyanophages identified the most promising engineering strategy for improving phage treatment—primarily the removal of the integrase to enforce a lytic cycle. However, the simulations also showed that even with substantially increased latent rates, cyanophage therapy would remain insufficient for neutralizing microcystins. Instead, the system tended toward coexistence between phage and Microcystis, allowing the harmful algal bloom and its problems to persist. These findings directed our wet‑lab efforts toward exploring an alternative chassis better suited for effective HAB suppression.
Background
Microcystis aeruginosa is a dominant bloom-forming cyanobacteria known for creating Harmful Algal Blooms (HABs) worldwide (Lim et al., 2023). These blooms disrupt aquatic ecosystems and release cyanotoxins, which have been implicated in adverse ecological impacts and linked to neurodegenerative diseases (Brenckman et al., 2025).
Cyanophage offers a novel solution to these problems through its unique ability to infect cyanobacteria and lyse the cell (Bhatt et al., 2023). However, cyanophage has various limitations, primarily in its lysogenic characteristics that make it difficult to successfully lyse a cell in a given environment (Guo et al., 2024). Currently, synthetic biology offers the ability to change these properties in a manner that could accelerate cyanophage’s capability of HAB treatment.
We looked to create a model of how this chassis would behave in flowing, aquatic environments to get a better understanding of the feasibility of treatment utilizing cyanophage. Models of cyanophage coupled with bacterial growth have been made, so we incorporated previous modeling frameworks to obtain a better understanding of the necessary design principles (Demory et al., 2020). We observed the following system as demonstrated with the following diagram:
Cyanobacteria grow logistically with rate α until reaching their carrying capacity KM. Upon infection of Cyanophage at rate ϕ, cells transition into an exposed and infected characterization. Granules are lysed by cyanophages according to a latent rate λ, releasing phage according to the burst size β and killing the cyanobacteria.
The cyanophage dynamics follow that it degrades over time according to a decay rate δ, is released and regrows according to a burst size β, and infects cyanobacterial cells according to an adsorption rate ϕ.
Based on these dynamics, we utilized the following system of ordinary differential equations:
Here, S denotes the susceptible cyanobacteria population, E the exposed population, I the infected population, and P the cyanophage population. The parameter lambda is proportional to the inverse of the latent period. While models often use a latent period, we instead reparameterize it as a latent rate (number of latent periods per hour) to avoid placing the parameter in the denominator.
Parameters
For the cyanobacteria, we look primarily at M. aeruginosa and thus incorporate respective parameter values:
Initial Conditions
For the initial conditions, we analyzed how a system with MOI of 1 with cyanobacteria and cyanophage concentrations of 107 CFU/mL and 107 PFU/mL, respectively.
Simulations
Initial simulations using Euler’s method (explicit) demonstrated the following cyanobacteria-cyanophage dynamics:
Given the initial parameters, the cyanophage is unable to suppress the bacteria, instead a coexistence occurs favoring M. aeruginosa sustainment, signifying the need for genetic engineering cyanophage to become capable of cyanobacteria treatment.We thus performed a parameter study to determine how lysis and adsorption rates influence algicidal treatment:
The graph highlights the critical role of a high latent rate. Natural cyanophages typically exhibit latent rates around 0.1, which is insufficient for effective treatment of M. aeruginosa. Therapeutic efficacy appears to require rates closer to 1. While adsorption rate also matters, increasing it by two orders of magnitude likely exceeds current genetic engineering capabilities—unless the phage can inhibit cyanobacterial defense gene pathways.
Our simulations indicated that for cyanophages to serve as a viable treatment for M. aeruginosa, extensive engineering would be needed to enhance their lytic activity. Nonetheless, the models revealed a parameter space where cyanophage deployment could successfully lyse M. aeruginosa.
Design Principles
Our simulations highlight key engineering strategies for enhancing cyanophage-based treatment of Microcystis. Increasing the adsorption rate would facilitate more effective suppression; however, achieving such a high rate through genetic engineering is likely infeasible. One potential approach involves inhibiting defense genes in M. aeruginosa, assuming adsorption can be modeled as a function of their expression:
Realizing this strategy would require extensive research to characterize differential gene expression in M. aeruginosa and to quantitatively link it to adsorption dynamics. In contrast, increasing the latent rate is more tractable. Removing the integrase from the cyanophage could accelerate lysis, shifting the latent period from days to hours. This would improve treatment efficacy, though it may reduce burst size due to shortened replication time.
Construct Deployment
After characterizing the cyanophage-bacteria dynamics and necessary modifications of cyanophage for M. aeruginosa treatment and HAB suppression, we looked to expand it to a larger scale model to determine whether or not cyanophage could sufficiently suppress M. aeruginosa in a flowing environment. We looked specifically at engineered lytic cyanophage based on our previous simulations to see if it could provide a plausible solution to M. aeruginosa under flowing water conditions.
Harmful Algal Blooms caused by M. aeruginosa form in large bodies of water, for example Lake Erie (Watson et al., 2016). The computational cost and complexities associated with making a full model, comprehensive of all vital factors and the fluid dynamics is extremely difficult, and more so difficult to interpret, so we looked to reduce the system’s complexity to get a better observation of how the chassis will behave when deployed.
Fluid Dynamic Assumptions
Due to the system being in a turbulent flow regime, the system becomes very complex and difficult to model. So, we reduce the system from modeling a full lake with varying complex geometries to a 1D system where the water transports the bacteria along the axis the fluid is flowing along.
For expanding this model to one on a macroscopic level, we make the following assumptions: The M. aeruginosa concentration is constant along its depth. This assumption holds generally with cyanobacteria where it is approximately constant in the epilimnion of a lake (Dokulil & Teubner, 2000); The M. aeruginosa concentration is constant along the width. This assumption does not generally hold due to the varying wave movements in a body of water, however it allows us to reduce the system to 1 dimension and help understand the feasibility of deployment. The cyanophage concentration is deployed at a constant concentration along the width of the body of water. The speed is averaged across the channel width, allowing a constant velocity profile along the body of water. Molecular diffusion is insignificant compared to the displacement induced by the water. This can be seen with a calculation of the Peclet number, where for these bodies of water it is much greater than 1, indicating advection dominates transport.
Chassis Transport Equation
Under these assumptions, the transport of the chassis can be done using the transport equations, reducing the model to a 1D Advection-Reaction system of PDEs:
Parameters
We modified the parameters such as the growth rate to be more indicative under real world conditions, where the carrying capacity of M. aeruginosa is less. We adjusted the latent rate to be reflective of a lytic cyanophage and due to the real world environment, the ability and probability for cyanophage to adsorb to M. aeruginosa is going to be less than what is observed in lab conditions. We also assumed that the burst size is significantly less due to the faster latent rate. We estimated this impact by decreasing the adsorption rate by an order of magnitude. We also scaled the parameters to be for the spatial case of meters:
Initial Conditions & Boundary Conditions
For the initial conditions of the simulation, we want to observe how the HAB would operate with and without the deployment of lytic cyanophage. So, we conducted a simulation with cyanophage and without and watched how the HAB developed over time. To initialize the M. aeruginosa population, we used a gaussian bump centered 100 meters away from the population of cyanophage, similarly initialized with a gaussian bump at 100 meters away from the origin. We additionally set σ for the M. aeruginosa to be bigger to initialize the HAB as taking up significantly more space than the deployment of the lytic cyanophage:
The attached populations were initialized to be 0 everywhere to properly simulate the event of chassis deployment and attachment.
For the boundary conditions of the system, we initialized an inlet of velocity but not a continual deployment of the chassis or M. aeruginosa. We are assuming that the deployment is upstream from the HAB location and will thus be transported toward the HAB.
We also implemented an explicit Euler scheme to update each grid point at every time step. To ensure numerical stability, the simulations satisfied the Courant–Friedrichs–Lewy (CFL) condition, expressed as:
Simulations of Chassis Deployment
The time slices above demonstrate the capability of cyanophage to successfully suppress a HAB under ideal conditions where it is able to thrive. There is an initial HAB suppression, however the M. aeruginosa is able to coexist with the phage and grow spatially. There additionally is a wave-like propagation of the M. aeruginosa profile that goes unsuppressed by the chassis, indicating that the HAB will still cause negative effects, though at a greatly decreased magnitude.
However, due to the nature of phage, ultimately the cyanobacteria will still release microcystins that are unaffected by the phage deployment. Additionally, the HAB is not completely suppressed and is nonetheless able to sustain itself, forming a front that the phage can not get to and treat. Phage is also limited in the ability to get a substantial amount of it to solve a system, getting enough phage to treat a HAB would be the depth of the HAB and length of the channel. Even a channel 10 meters wide would imply a need for at least 60,000 liters of phage, which given the growth rate of cyanophage is near impossible. However, a cyanolytic bacteria such as Acinetobacter, could nonetheless be grown both much faster and require a lower concentration, making it a much more feasible treatment solution.
Conclusion
In the conditions where lytic cyanophage is a chassis fit for the aquatic environment it is deployed into, it could offer many benefits to HAB suppression. However, the cyanophage is unable to fully suppress the M. aeruginosa population or limit the production and exposure of microcystins significantly. In fact, though cyanophage suppresses a majority of the HAB, the simulations illustrate an inability of the cyanophage to contain and suppress the entirety of the HAB, implying that the detrimental effects of the HAB are nonetheless sustained.
This informed and guided our wet-lab to look into another chassis, Acinetobacter, as a possible treatment method for HABs, as it has well documented cyanolytic potential and furthermore an ability to prevent microcystins from aquatic environments (Li et al., 2016).
The transport model can be modified to fit any range of bacterial dynamics to get a preliminary understanding of how a chassis will behave in moving water environments where advection dominates the transport.
Future Directions and Limitations
Though this model demonstrates the feasibility of various chassis in water and can give a predictive approach to deployment, further enhancements can be made through the usage of computational fluid dynamics and solvers that would allow for a more accurate analysis of chassis deployment in specific fluid dynamic environments.
Executive Summary
Mathematical modeling suggested that Acinetobacter strains are better suited for HAB suppression, requiring far lower concentrations than cyanophages to achieve comparable effects. Guided by this insight, our wet‑lab conducted microcosm experiments to observe the interaction dynamics between Acinetobacter and M. aeruginosa. These experiments provided critical data on whether Acinetobacter could sustain itself in moving, dynamic water systems, which in turn informed the next stage of our modeling. Using these results, we extended the framework to simulate aquatic environments in silico, allowing us to further test the potential of Acinetobacter for effective M. aeruginosa and HAB suppression.
Background
Microcystis aeruginosa is a dominant bloom-forming cyanobacteria known for creating Harmful Algal Blooms (HABs) worldwide (Lim et al., 2023). These blooms disrupt aquatic ecosystems and release cyanotoxins, which have been implicated in adverse ecological impacts and linked to neurodegenerative diseases (Brenckman et al., 2025).
Acinetobacter has been reported to exhibit cyanolytic activity, with the ability to lyse cyanobacteria such as M. aeruginosa, making them promising candidates for the biocontrol of M. aeruginosa and furthermore HAB suppression (Li et al., 2016). The interaction typically begins with Acinetobacter adhering to cyanobacterial cells, followed by the recruitment of additional bacterial cells to form aggregates on the membrane. Over time, Acinetobacter produces and releases enzymes that deform and ultimately lyse the cyanobacterial cells. This process is also associated with inhibition of photosynthesis, giving Acinetobacter an additional light‑suppressing effect on M. aeruginosa (Qin et al., 2016).
We developed a novel system of ordinary differential equations to describe the Acinetobacter-M. aeruginosa population dynamics, as illustrated in the following diagram:
M. aeruginosa grows logistically with rate α until reaching their carrying capacity KM. Upon attachment of Acinetobacter at rate ϕ, cells transition into a granule state with no growth due to shutdown of photosynthetic pathways. Granules are lysed by Acinetobacter at rate λ, releasing attached cells and killing the cyanobacteria.
Acinetobacter itself grows at rate β until reaching its carrying capacity KA. The cyanolytic interaction is governed by the attachment rate ϕ, after which attached Acinetobacter continue reproducing at rate β. Accordingly, we classify Acinetobacter into two states: unattached (Au) and attached/granule-associated (Ag). Based on these dynamics, we formulated the following system of ordinary differential equations:
Here, Mu represents the unattached Microcystis aeruginosa population, and Mg the granulated Microcystis aeruginosa population. Au denotes the unattached Acinetobacter population, which grows competitively with other bacteria, decreases through attachment, and is replenished upon lysis of granules. Ag denotes the attached (granule‑associated) Acinetobacter population.
Parameters
We estimated the lysis rate based on the typical time required for Acinetobacter to lyse the Microcystis aeruginosa cell membrane (Osman et al., 2017). The Microcystis aeruginosa growth rate was scaled to reflect laboratory conditions, while the attachment rate was approximated from cell motility, cell radii, and an attachment probability influenced by the extracellular polymeric substances (EPS) of the Acinetobacter outer membrane, based on the data presented by Li et al.:
Initial Conditions
For the initial conditions, we adopted values reported by Li et al., which conducted cocultures of Microcystis aeruginosa and Acinetobacter to simulate their interaction dynamics. The initial population values for the system are as follows: The unattached Microcystis aeruginosa starts at 1.5 × 106, while the granule/attached Microcystis aeruginosa begins at 0 CFU/mL. For Acinetobacter, the granule/attached population is initially 0 CFU/mL, and the unattached population starts at 1 × 103 CFU/mL.
Simulations
Initial simulations using Euler’s method verified consistency with reported dynamics of M. aeruginosa populations in coculture with Acinetobacter (Li et al., 2016):
The agreement of experiment and model justifies the chosen parameters (see table above) thus we performed a parameter study to determine how lysis and attachment rates influence algicidal treatment.
The graph above indicates the importance of Acinetobacter attachment for proper treatment of theM. aeruginosa. Lysis rate plays an important role, but the attachment rate of the Acinetobacter is the determining factor of M. aeruginosa suppression.
The simulations suggested to our wet lab that Acinetobacter could be a promising treatment for systems containing M. aeruginosa under the right conditions. However, Acinetobacter and M. aeruginosa behave differently in real environments compared to cocultures. Thus, it is essential to account for the influence of water system fluid dynamics to better predict their interactions and extend the model to real‑world applications.
Design Principles
Our simulations identify key engineering strategies for optimizing M. aeruginosa treatment with Acinetobacter. Increasing the attachment rate is critical, as it depends on the adherence probability to cyanobacterial surfaces and is linked to genes regulating extracellular polysaccharide (EPS) production. Overexpressing these genes is therefore central to success in aquatic systems. With further knowledge to connect qualitative data of protein production, specifically in the K locus, we could create a quantitative form of relating gene expression of Acinetobacter to the attachability of Acinetobacter (Geisinger & Isber, 2015).
The lysis rate, governed by the genetic pathway for enzyme production that degrades cyanobacterial cells, also contributes to efficacy. However, attachment rate and EPS production exert a far greater influence on system treatability than lysis alone.
Construct Development
After characterizing the bacteria dynamics of Acinetobacter and M. aeruginosa, we looked to expand it to a larger scale model to determine whether or not Acinetobacter could sufficiently suppress M. aeruginosa in a flowing environment.
Harmful Algal Blooms caused by M. aeruginosa form in large bodies of water, for example Lake Erie. The computational cost and complexities associated with making a full model, comprehensive of all vital factors and the fluid dynamics is extremely difficult, and more so difficult to interpret, so we looked to reduce the system’s complexity to get a better observation of how the chassis will behave when deployed.
Fluid Dynamic Assumptions
Due to the system being in a turbulent flow regime, the system becomes very complex and difficult to model. So, we reduce the system from modeling a full lake with varying complex geometries to a 1D system where the water transports the bacteria along the axis the fluid is flowing along.
For expanding this model to one on a macroscopic level, we make the following assumptions: the M. aeruginosa concentration is constant along its depth. This assumption holds generally with cyanobacteria where it is approximately constantly concentrated in the epilimnion of a lake (Dokulil & Teubner, 2000). The M. aeruginosa concentration is constant along the width. This assumption does not generally hold due to the varying wave movements in a body of water, however it allows us to reduce the system to 1 dimension and help understand the feasibility of deployment. The Acinetobactor concentration is deployed at a constant concentration along the width of the body of water. The speed is averaged across the channel width, allowing a constant velocity profile along the body of water. Molecular diffusion is insignificant compared to the displacement induced by the water. This can be seen with a calculation of the Peclet number, where for these bodies of water it is much greater than 1, indicating advection dominates transport.
Chassis Transport Equation
Under these assumptions, the transport of the chassis can be done using the transport equations, reducing the model to a 1D advection-reaction system of PDEs:
Parameters
We modified the parameters such as the growth rate to be more indicative under real world conditions, where the carrying capacity of M. aeruginosa and Acinetobacter is less. Additionally, due to the real world environment, the ability and probability for Acinetobacter to adhere to M. aeruginosa is magnitudes less due to the other bacteria present in the lake, disrupting attachment. We estimated this impact by decreasing the probability of attachment and thus attachment rate by an order of magnitude. We also scaled the parameters to be for the spatial case of meters:
Initial Conditions & Boundary Conditions
For the initial conditions of the simulation, we want to observe how the HAB would operate with and without the deployment of our chassis, Acinetobacter. So, we conducted a simulation with Acinetobacter and without and watched how the HAB developed over time.
To initialize the M. aeruginosa population, we used a gaussian bump centered 100 meters away from the population of Acinetobacter, similarly initialized with a gaussian bump at 100 meters away from the origin. We additionally set σ for the M. aeruginosa to be bigger to initialize the HAB as taking up significantly more space than the deployment of Acinetobacter:
The attached populations were initialized to be 0 everywhere to properly simulate the event of chassis deployment and attachment.
For the boundary conditions of the system, we initialized an inlet of velocity but not a continual deployment of the chassis or M. aeruginosa. We are assuming that the deployment is upstream from the HAB location and will thus be transported toward the HAB.
We also implemented an explicit Euler scheme to update each grid point at every time step. To ensure numerical stability, the simulations satisfied the Courant–Friedrichs–Lewy (CFL) condition.
Simulations of Chassis Deployment
The time slices above demonstrate the ability of Acinetobacter to successfully suppress a HAB under conditions where it is able to thrive.
With the knowledge that attachment rate will dominate the sustainability of the chassis, we created the following heatmap to observe how the attachment rate changed the final concentration profile of the M. aeruginosa:
From this, it is seen that for complete treatment of the M. aeruginosa, an attachment rate around 10-14 m/h is necessary for HAB suppression. This adsorption rate is possible naturally, however, due to the many factors that come about in water, the bacteria's success becomes much more volatile and the need for genetic modification of the exopolysaccharide production becomes necessary for proper implementation in ideal conditions.
Conclusions
In the conditions where Acinetobacter is a chassis fit for the aquatic environment it is deployed into, it could offer many benefits to HAB suppression. Furthermore, genetic engineering of the Acinetobacter to overexpress genes related to EPS production is necessary for proper treatment.
The transport model can be modified to fit any range of bacterial dynamics to get a preliminary understanding of how a chassis will behave in moving water environments where advection dominates the transport.
Future Directions & Limitations
Though this model demonstrates the feasibility of various chassis in water and can give a predictive approach to deployment, further enhancements can be made through the usage of computational fluid dynamics and solvers that would allow for a more accurate analysis of chassis deployment in specific fluid dynamic environments.
Pipes & Plumbing
Chassis Transport Equation
The transport of a synthetic biology chassis through a plumbing system, where it reacts with a surrounding substance, can be described by a partial differential equation governing its temporal evolution as it disperses or diffuses within the fluid. For example, the chassis concentration can be modeled as a function of the radial coordinate under the assumptions of azimuthal and translational symmetry to understand its ability to be correctly transported within a system. This framework captures both passive transport (diffusion, advection) and active interactions (adsorption, reaction, or colonization).
In many engineered environments, eddy diffusion dominates molecular diffusion. We approximate this using standard chemical engineering relations, where the dimensionless constant α characterizes the flow regime, u is the averaged bulk flow speed, and R is the pipe radius (e.g., α = 0.05 in pipe flow; Venkatram & Weil, 2021).
At the bulk–surface interface (e.g., biofilm, tissue, or material surface), reaction rates govern how the chassis interacts with its environment. This can be incorporated through the following mixed boundary condition,
The relevant timespan is determined by the interaction window between the flowing chassis and the target surface. This can be calculated from the fluid velocity u and the characteristic length of the surface region, L, given by the length of the reactive surface.
The system can be initialized at various chassis concentrations, assumed spatially uniform, and scaled to the geometry of the system.
This transport framework applies broadly to any synthetic chassis deployment into a plumbing environment or pipe system (household plumbing, tubes in medical constructs, etc.). By explicitly incorporating fluid dynamics into chassis design, synthetic biologists can better predict how constructs disperse, persist, and interact with surfaces or communities. Furthermore, they can identify engineering targets that improve performance in specific environments, such as adhereability of the construct. Through this framework, synthetic biologists can predictably optimize deployment strategies by matching chassis properties to flow regimes, geometries, and ecological contexts.
Executive Summary
The math modeling demonstrated the feasibility of phage in suppressing biofilms when deployed at sufficiently high titers, showing that under controlled conditions they can be an effective treatment strategy, a promising result for our wet-lab to test. However, insights from the wet‑lab highlighted resistance evolution as a significant limitation and the difficulty not accounted for in diffusion and adsorption of phage in transport, feeding back into the modeling to emphasize long‑term sustainability concerns. The final analysis revealed that while phage can function effectively in pipe systems, challenges arise at lower titers where stochastic effects with adsorption and resistance reduce reliability. The key takeaway is that phage are indeed viable in controlled environments, but their success ultimately depends on strategies to mitigate resistance and enhance adsorption.
Background
Phage–bacteria dynamics have long been studied, and many models describe these interactions (Demory et al., 2020). Phages are viruses that infect bacterial cells, replicate within them, and are released when the host cell bursts. Each phage has a specific host range, meaning it can only infect compatible bacteria.
A common framework divides bacterial populations into three groups: susceptible, exposed, and infected (Palacios & Segal, 2025). This structure parallels the SEIR model, where phage replication within host cells and release upon lysis drive changes in both bacterial and phage populations. Because phage abundance directly influences infection rates, the phage population is also tracked closely. The system’s dynamics can be captured by the following diagram:
The population dynamics are shaped by interactions between phage and bacteria. We first examine phage behavior within this system.
Phages infect host cells indiscriminately, so their loss through infection depends on the total number of host cells rather than just the susceptible population. This infection rate is defined by the adsorption parameter ϕ. Free phages also degrade at a rate δ. Propagation occurs when an infected cell lyses, releasing on average β new phages per event.
Bacterial dynamics depend on phage presence. In the absence of phage, bacteria remain susceptible and grow toward their carrying capacity at a maximal rate α. When phages are introduced, susceptible bacteria transition to exposed and then to infected states. While this intermediate “exposed” stage may seem unintuitive, it provides a more realistic, gradual representation of infection and replication; without it, the system evolves unrealistically fast. The transitions from exposed to infected, and from infected to lysis and death, are governed by the latent rate λ, which reflects the time required for phage-induced lysis.
Based off this information, we observed how the following system of ordinary differential equations behaves:
Here, S denotes the susceptible bacteria population, E the exposed population, I the infected population, and P the phage population. The parameter lambda is proportional to the inverse of the latent period. While models often use a latent period, we instead reparameterize it as a latent rate (number of latent periods per hour) to avoid placing the parameter in the denominator.
Parameters
We drew parameter values from the literature for mycobacteria and mycobacteriophage, ensuring consistency in system dimensionality.Concentrations are expressed with spatial units of 10 microns, effectively normalizing bacterial populations between 0 and 1. Using dimensional analysis, we converted literature values to this scale for accurate calculation and analysis.
Initial Conditions
For the initial conditions, the bacterial population is set at 1 to represent maximal concentration, while the exposed and infected populations are initialized at 0 to capture post‑infection dynamics. The initial phage population is 1 plaque-forming unit ( PFU) per pL.
The initial conditions for the system are as follows: the concentration of susceptible bacteria (S) is 1 CFU per picoliter, while both exposed (E) and infected (I) bacteria begin at 0 CFU per picoliter. The phage population (P) starts at 1 PFU per picoliter.
Simulations
The resulting system of nonlinear, cross‑coupled ODEs is analyzed through a parameter study, focusing on terms modifiable by environment or genetic engineering. For numerical simulation, we used Euler’s method and enforced positivity. We analyzed the system’s initial behavior across biologically realistic parameter ranges, beginning with the adsorption rate, which reflects phage infectivity. Simulations revealed oscillations between phage and bacterial populations: phages initially suppress bacterial growth, then degrade, allowing bacteria to rebound exponentially, before being suppressed again. These oscillations were driven primarily by the adsorption rate, so we conducted simulations to identify the parameter values that determine whether the system oscillates indefinitely or converges to a steady state.
The above plot illustrates the range of adsorption rates that determine whether the system reaches a steady state—indicated by near‑zero oscillation amplitude—or undergoes sustained oscillations with large amplitudes. Results show that phages with adsorption rates near 10-1 produce indefinite oscillations, while lower rates drive the system toward steady state. We therefore examined the steady states achieved across adsorption rates that favor convergence.
The above graph shows that phage–bacteria coexistence occurs at adsorption rates between 10-3 and 10-1. Within the narrower range of 10-2 to 10-1, the system tends toward successful phage treatment, where bacterial populations are effectively suppressed.
This serves as a precursor for evaluating whether a phage can adsorb to and lyse host cells sufficiently before application in real systems. This analysis highlights how adsorption rate influences treatment outcomes and provides insight into the parameters that govern phage efficacy against biofilms.
We additionally did a sweep over the latent rate and burst size to get a better understanding of the treatment of phage. To do this, we iterated over the range of latent rates and burst sizes, fixing adsorption rate to be 0.01 pL/h:
As expected, phages with longer latent periods require larger burst sizes to outcompete bacteria. The results demonstrate the necessity of lytic biased phages and demonstrate that burst size becomes a predominant factor once the latent rate reaches a certain point: even with extremely high burst sizes, lysogenic phages fail to suppress bacterial populations. Using literature‑based parameters for mycobacteriophage concentration, we then tested whether our phage solution would be effective.
After collecting all the necessary values from the literature, we conducted the following simulations under an appropriate time scale, we chose a time span of approximately one month.
The results initially suggested that phage, at sufficiently high concentrations, could suppress bacterial populations. However, over longer time scales the bacteria recovered and restabilized, indicating that an additional mechanism for bacterial removal may be necessary for sustained treatment. At these scales, other factors also become significant—for example, in pipe systems, constant water flow could wash phage into the effluent, enabling bacterial regrowth.
Design Principles
Our simulations highlight key engineering strategies for optimizing phage therapy. Increasing the adsorption rate—potentially by inhibiting bacterial defense protein expression—enhances infection efficiency. Though data supporting this claim is too limited for quantitative analysis, qualitatively this equation relates the differential gene expression of defense proteins to the adsorption rate of phage:
More critically, removing the integrase from a phage to bias a lytic lifecycle is essential, as even high burst sizes cannot compensate for lysogeny. While increasing burst size further improves efficacy, integrase removal remains the dominant factor for successful bacterial suppression.
We additionally recognize the fact that the matrix is inhibiting phage infection and the development of enzymatic phage could allow phage infection to excel. According to this, we developed the following set of equations describing enzymatic phage infection where Z is the enzyme concentration, M is the extracellular matrix sugar concentration and A is the attached-enzyme-matrix concentration:
γ represents the production rate of the matrix degrading enzyme, μ is the matrix production rate and k1 is the forward reaction rate of the enzyme attaching to the enzyme, k-1 is the respective reverse reaction rate, and k2 is the reaction rate of the attached-enzyme-matrix leading to the destruction of the matrix, release of the enzyme, and production of an inert product with respect to the reaction.
However, this is only known on a qualitative basis for M. smegmatis biofilms, though a quantitative basis is possible with further research, especially with differential gene expression and how the expression of genes responsible for the production of the extracellular matrix are up or downregulated. This provides a framework for enzyme-phage modeling with the potential to solve biofilms.
Construct Development
After characterizing the phage–bacteria dynamics for M. smegmatis and mycobacteriophage, we next examined phage transport to the M. smegmatis biofilm through pipes.
We first calculated the Reynold’s number of the system to better understand the flow regime the transport would occur in. In practical plumbing systems, such as the ones found in everyday households, the fluid speed tends to be at least 0.5 m/s; our phage concentration will approximately be characteristically similar to water and thus we use parameters that respectively represent water. In plumbing systems such as household pipes, fluid speed (u) is typically at least 0.5 m/s. Since the phage solution behaves similarly to water, we use water-based parameters for our calculations. Using a pipe diameter of 0.2 cm, we calculated the Reynolds number:
With this in mind, we now look to model how the phage would be transported in this turbulent regime, however incorporating all of the factors that come about in turbulent flow is very complicated; so, we looked at the most important factors that come about in a turbulent flow regime.
Due to the system being in a turbulent flow regime, the system becomes very complex and difficult to model. So, we reduce the system from modeling the full flow within a pipe to monitoring how one spatial instance of the phage concentration operates when in the presence of the biofilm’s surface. We can model this as a reduced system through the following assumptions: The phage concentration is well mixed and the concentration is approximately constant; phage adsorption is uncompetitive because when one phage adsorbs to the surface of the biofilm, that does not impact the ability of another phage to adsorb to the surface of the biofilm; the biofilm height and bacterial density is constant along its length, we let this be an overapproximation of the amount of bacteria, thus giving a “worst case” scenario for the phage deployment; phage infection and reproduction is ignored, thus we only look at how the initial amount of phage is able to behave and adsorb. This is because the time scale of the system is significantly less than the time scale required for phage infection and cell lysis.
With these assumptions in mind, we can adopt a model to describe the system through the usage of a partial differential equation that represents how the phage is being transported and reacting to the surface of the biofilm. Below is a diagram portraying the geometry of the system:
This calculation indicates that the flow regime of the system will likely be turbulent rather than laminar flow. Using cylindrical (axi-symmetric) coordinates and assuming that the phage concentration profile is uniform in the transverse (z) direction, the dependent variable of the model is the radial concentration profile of phage, denoted P(r,t).
Fluid Dynamic Assumptions
Due to the system being in a turbulent flow regime, the system becomes very complex and difficult to model. So, we reduce the system from modeling the full flow within a pipe to monitoring how one spatial instance of the phage concentration operates when in the presence of the biofilm’s surface. We can model this as a reduced system through the following framework previously described.
Chassis Transport Equation
Under these assumptions, phage transport can be analyzed with a partial differential equation describing its time evolution as it diffuses along the pipe’s radial axis, assuming azimuthal symmetry and a transitional symmetry:
Parameters & Initial Condition
For the diffusion constant, we assume eddy diffusion dominates molecular diffusion. We approximate this process using a standard chemical engineering relation, where α is a dimensionless constant characterizing the flow regime, in pipe flow, α = 0.05 (Venkatram & Weil, 2021).
At the bulk–biofilm interface, there will naturally be a reaction rate determining how the phage infects the biofilm. This will be governed by the surface concentration of bacteria and adsorption rate of phage:
Due to the geometry of the system, the time span of the simulation is when the phage concentration is interacting with the biofilm surface. This will be calculated through the following formula, relating the water speed and length of the biofilm to the timespan of simulation:
We initialized the system at various phage concentration values that are constant in space according to the following equation:
Due to the geometry of the system, the phage concentration in the simulation is attained by the following calculation where c is the phage concentration per cubic meter and R is the pipe radius:
We then use the following parameters that describe mycobacteriophage and mycobacteria:
Boundary Conditions
At the center of the pipe, there should be no flux given the azimuthal symmetry so we utilize the following Neumann boundary condition:
At the interface between the biofilm and bulk phage concentration, we implemented a Robin mixed boundary condition that incorporates the adsorption rate previously calculated according to the following equation:
Numerical Method
The PDE with a Robin boundary condition was solved using an explicit finite‑difference scheme: forward Euler in time and central differences in space with N uniform grid points. Stability was ensured by enforcing the CFL condition, maintaining the required inequality for time‑step size relative to spatial discretization.
Simulations of Chassis Deployment
The plot shown above is the phage concentration (PFU per unit length) over time with an initial concentration corresponding to a concentration of 1012 PFU/mL. Limited adsorption to the biofilm indicates a reaction‑limited system, though reaction is occurring.
The graph indicates an approximately linear increase in attached phage concentration over time, with the product of concentration and volume remaining constant. This is best seen through the log plot where a concentration that is ten times more concentrated behaves ten times better than the diluted concentration. In other words, flushing ten times with a 10-1 dilution of my original concentration will behave the same as if I did not dilute my concentration. This suggests that treatment efficacy depends on the total phage present rather than concentration alone.
These simulations provide a generalizable framework for phage transport in pipes, adaptable to different diameters and flow speeds, and highlight how phage adsorption to biofilms can be understood under fluid dynamic conditions. While simplified by assumptions, the model captures key features of phage diffusion and reaction, supporting its applicability. Overall, the results show promising adsorption dynamics in flow environments and continue into our next component of monitoring how the phage reacts within the biofilm after transport.
Phage-Biofilm Infection
The extracellular matrix of the biofilm creates spatial dependence in phage infection, as phages penetrate the biofilm but have reduced access to bacteria at greater depths. To capture this, the model is expanded from ordinary differential equations to partial differential equations that vary in both space and time.
Since infection occurs after the concentration flush, the PDEs exclude fluid advection and instead focus on phage diffusion and reaction within the biofilm. Because biofilm height is much smaller than the pipe radius, it suffices to model penetration of phage into biofilm as a 1D Cartesian reaction–diffusion system, leading to the following set of partial differential equations for dynamics of phage and bacteria:
Diffusion Coefficients
Diffusion coefficients for M. smegmatis (Db) and phage (Dp) were chosen to reflect hindered movement in the biofilm matrix. We assume a reduction of two orders of magnitude compared to free media (e.g., water), yielding the following values scaled to our spatial and temporal dimensions while assuming phage diffuse an order of magnitude faster than bacteria:
Initial Conditions & Boundary Conditions
The system was initialized with bacteria near carrying capacity, while exposed and infected populations were set to zero prior to phage infection.
To avoid numerical instabilities from initializing phage at a single point, we approximated the initial condition with a Gaussian bump at the origin. Its magnitude was set to the average phage concentration adsorbed after the pipe transport simulation, calculated by averaging the total adsorbed phage over both distance and circumference of the biofilm surface.
The system uses Neumann boundary conditions enforcing zero flux at both the biofilm surface and its maximal depth at the pipe wall. These conditions ensure that populations cannot leak through the pipe material or extend into regions without matrix (i.e., into the air):
We additionally used a Forward-Time-Centered-Space numerical method to simulate the system of partial differential equations and ensured numerical stability using the CFL condition mentioned above.
Simulations
The heatmaps show mycobacteria-mycobacteriophage interactions over time, highlighting phage diffusion and reaction within the biofilm that ultimately leads to infection of the bacterial population.
The heatmaps show the temporal dynamics of exposed and infected mycobacteria lead to a propagating pulse along the bacteria–phage interface.
The above simulation shows that the SEIR population dynamics and diffusion of phage is sufficient to suppress the mycobacteria throughout the depth of the biofilm. The velocity of the simulated wave front is approximately 2 microns an hour.
Conclusion
The numerical studies presented above suggest that phage can potentially treat bacterial biofilms in pipes, provided bacteria do not develop resistance. Identifying the genes responsible for resistance could enable genetic engineering strategies to inhibit their activation, leading to more robust phage‑based treatments. This gave our wet-lab a promising result, however dilute phage concentrations could cause adsorption to become more stochastic in nature and thus the model would not be able to properly represent these occurrences. However, at concentrations of phage at a high enough titer, the model becomes more deterministic in nature and more representative of reality. Additionally, the model does not fully incorporate the adaptability of the mycobacteria to become resistant to the phage.
Future Direction & Limitations
The model could be adapted to different hydrodynamic conditions of pipe flow, and used as a starting point for simulating alternative treatment methods or synthetic biology chassis, capturing their transport and dynamics to assess treatment efficacy. Further modeling using Navier-Stokes equations to represent fluid dynamics would be an important future step to get the best understanding of how the various chassis would be transported and react in plumbing systems. Additionally, further research into the exact quantitative implications of gene expression on phage infection ability is needed for model implementation.
Material Attachment
For deployment of a synthetic biological material, we reduce the complexities of the system so that the material is able to be simply modeled by shear and normal detachment through an assumption of homogeneity.
Due to the fluid dynamics that interact with the material, we use simple equations of stress and strain and contact angles, θ, to understand the variable deployment positions that favor sustainability. This composes our shear, τ, and normal, σn stress equations into the following dependent on the aquatic environment, where L is the material thickness, ρ is the fluid density, μ is the viscoelasticity, and u is the speed of the fluid:
We can then look to see the various contact angles and properties of water where deployment is possible or impossible. Additionally, we can determine the critical shear stresses and critical normal stress of a given material to find out what genetic design principles need to be applied to improve various elements of the construct. Here, ks and kn are detachment coefficients and n and m characterize the level of detachment after reaching the threshold detachment value:
Additionally, when a construct is deployed into an environment, various other factors can affect the surface structure, such as ion degrading a surface or other solutes reacting with the surface. This can be incorporated into the equations using a Michaelis-Menten type of reaction following this scheme, where S is the reactive concentration in the aquatic environment and E is the interface of the material:
With these kinetics in mind, a preliminary model of deployment can be made to discover the various aquatic conditions and environments where deployment is feasible or infeasible. With this framework of model, genetic design principles can be made to understand the environments where a reactive substance or the fluid properties cause the construct to not sustain its deployment.
Executive Summary
The mathematical modeling indicated key design principles for biofilm deployment, indicating the necessity to raise the critical shear and normal stresses to ensure viable deployment. Our wet-lab then conducted various genetic engineering on B. Subtilis in genetic pathways that would increase the feasibility of deployment. We overall discovered that B. Subtilis could be deployed in oceanic environments so long as it is able to self-sustain itself when deployed.
Background
Bacillus subtilis biofilms have been shown to inhibit microbiologically influenced corrosion (MIC) in marine environments by forming a stable biomineralized film that protects steel surfaces (Guo et al., 2019). With this in mind, we looked to see how B. Subtilis biofilms could be modified and whether it can sustain itself in marine environments. Biofilm growth and detachment are pivotal components of understanding the deployability of biofilms into aquatic environments to prevent corrosion (Picioreanu et al., 2000). It is critical to understand the range of environments in which a deployed biofilm could sustain itself, as opposed to detaching or being gradually degraded by the surrounding fluid. As a starting point, we developed a minimal model of the biofilm height at a particular location.
An ordinary differential equation describes the logistic growth dynamics of the biofilm thickness over time.
Where the parameter μg determines the time scale for biofilm thickening, and Lmax is the maximal thickness of the biofilm given by 4 x 10-5/h and 0.08 cm, respectfully.
Simulations
The graph on the left was from data collected by (Wang et al., 2015) the graph on the right is our initial simulations for B. Subtilis’ thickness based on a logistic growth. After this initial validation, we incorporated the various parameters that would be implied when deployed in a fluid dynamic environment.
Fluid Dynamic Assumptions
Due to the high content of exopolysaccharides and sugars that compose of the biofilm, the material will be reacting with the bulk environment primarily depending on the salinity of the water. This is incorporated into the ODE model for biofilm height using as a loss term that is proportional to biofilm height a Michaelis-Menten type of reaction between the biofilm interface and the water:
Here, the salinity of the water is treated as a constant due to the rapid influx and outflux of incident waves upon the biofilm surface. kb is the salt damage constant and Ks is the half saturation constant for salt damage.
Additionally, the B. Subtilis must be able to sustain themselves through the stresses and strains the bulk flow would induce onto the biofilm. We analyze in particular how the normal and shear stresses would cause detachment of the biofilm. To do this, we introduce a parameter θ that represents the average angle of incidence upon the biofilm surface. This angle determines the component of normal stress and shear stress. From that, we introduce a shear detachment term due to the shear stress with the respective term, where τ is the shear stress and τcrit is the critical shear stress value, where once τ is greater than τcrit shear detachment occurs:
Additionally, we introduce a normal detachment term for the normal stress where σn denoted the normal stress induced on the biofilm and Pcrit is the critical normal stress that would cause detachment of the biofilm:
With these additional terms in mind, we expanded the previous model for logistic growth into the following differential equation:
Parameters
We used the following parameters to properly simulate the environments of deployment. Since the biofilm will be deployed in environments where corrosion is affecting the environment which heavily occurs in salt water environments, we let the conditions pertaining to the bulk flow reflect this through Sconst and a modification of the fluid density. With that in mind, we collected the following parameters to conduct simulations with ranges being realistic:
Construct Deployment
We then conducted various parameter studies to see how the biofilm would behave in various aquatic systems. First we determined how the flow velocity and contact angle determines the feasibility of deployment by doing a parameter sweep at various flow velocities and contact angles.
The above graph demonstrates the various circumstances that the biofilm would be effective at adhering and sustaining itself. At contact angles that favor normal stress, the biofilm will have a less sustainable profile. With this, we determined next to look at how the salinity and growth rate of the biofilm determined the feasibility of sustained growth with a similar approach by creating a heatmap of a parameter sweep over various growth rates and salinities.
The above graph demonstrates the various circumstances that the biofilm would be effective at adhering and sustaining itself. It is easily seen that over various realistic salinity values, the resulting sustainability is invariant. The predominant factor is thus seen to be the growth rate of the biofilm as to whether it will sustain itself or not.
The above graph portrays the values of critical stress and normal stress necessary for sustained deployment. Due to the heaviside nature of the detachment terms, it is no surprise that there seems to be a threshold value for both critical shear and critical normal stress needed for deployment and a near discontinuous nature. These are the predominant values that must be modified through genetic engineering.
Conclusion
Through these simulations we were able to determine the various states and situations where biofilm deployment is a feasible option and determined key parameters to modify to ensure deployment success. Based on these simulations, the salinity of the water isn’t an extremely important factor; additionally, the growth rate of the biofilm, especially being nonzero, is necessary for a sustained approach.
It is clear that the most important factors, aside from the deployment technique, are the ability of the B. Subtilis to grow at a net positive rate and that the critical shear stress reaches a threshold value to ensure detachment occurs at an infrequent rate. With this in mind, it thus makes sense to upregulate pathways that lead to a higher growth rate and EPS excretion while additionally upregulating genes that would cause a higher structural integrity for managing the shear stresses. The genetic pathways that would help the sustainability of the biofilm were thus genes that upregulated biofilm related genes, ones that ease the shear stress through a more hydrophobic surface, or genes related to the production of amyloid fibers to help critical detachment parameters.
Future Direction & Limitations
Though this model portrays characteristics that would help understand the genetic engineering design principles and the necessary angles and abiotic factors of deployment, further structural analysis that incorporates computational fluid dynamics and the usage of stress tensors would lead to a better understanding of the behavior of deployment for the system.
This differential equation works well for homogenous structures as each point can be treated more or less the same, however in nature there will be discrepancies in the homogeneity of the biofilm leading to structural weak points that must additionally be accounted for and treated differently depending on that.
Categories of Aquatic Systems
AQUIRE-ING MACHINE LEARNING INTO FEASIBILITY ANALYSIS
With a foundational understanding of how the various deployment techniques and situations can be treated, a commonly overlooked aspect of synthetic biology is the natural compatibility of a bacteria to sustain itself in the environment into which it is deployed.
Therefore, a predictive model is needed to characterize a construct’s capacity for self-sustainability, enabling its design to be tailored to specific environments and ensuring effective deployment where it can persist. For our model to be all encompassing of the ability and predictive aspect of synthetic biology constructs in model, we couple our database model AQUIRE to incorporate the biotic and abiotic factors that an ecosystem would imply on a chassis’ success.
We do this through an adjustment of the constructs carrying capacity (K) and growth rate (α) based on the fitness parameter assigned by AQUIRE based on the following mapping where γ is the fitness parameter of the chassis:
To demonstrate the importance of this parameter for predictive approaches of chassis deployment, we first look at the differences of behavior of Acinetobacter under differing γ values.
At gamma values below 0.5, the HAB remains above the threshold required for full suppression. Although the chassis functionally reduces HAB levels, key deterrents such as microcystins are unsolved and persist into further waterways. This indicates that synthetic biologists if willing to deploy in a less fit environment must compromise or engineer further their construct to be specifically effective for a given environment.
We additionally applied this transformation through a gamma term onto the B. Subtilis deployment model to better understand the environments where it can sustain itself properly.
We additionally see that the ability of the B. Subtilis biofilm depends greatly on its level of fitness in an environment. If its fitness isn't high enough, it will not sustain itself and will thus be dysfunctional within its environment.
While the precise link between chassis fitness and environmental performance remains unclear, we incorporate this scalar into our model to have a predictive understanding of the deployability of a chassis.
The feasibility of any synthetic biology chassis depends not only on its engineered traits but also on its ecological fitness within the deployment environment. AQUIRE’s machine‑learning–derived fitness parameter allowed us to quantify this compatibility and incorporate it into predictive simulations, making the modeling both more realistic and more actionable. This combined framework ensures that future chassis design can be tailored not just for function, but for persistence and effectiveness in their targeted real‑world systems.