Optimal flow analysis

Introduction
In our iGEM project we engineered bacteria to break down toxic polycyclic aromatic hydrocarbons (PAHs) [1] to efficiently clean rivers. To ensure biocontainment, the bacteria shall be immobilized inside a device. Polluted water and native bacteria will pass this device, giving the GMOs the opportunity to degrade PAHs. To ensure effective cleanup, this device needs to be designed in a way, that polluted water comes in optimal contact with immobilized bacteria, while guaranteeing steady flow through of water.

Why create a model?

To better understand the behavior of a mechanism and to determine outcomes in advance, we need to implement a model. For simple, linear problems, such as mechanical equations or differential equations, we have tools that allow us to create a very good — and in some cases, perfect — analytical model of the process. Unfortunately, most real-world problems are not simple. They involve losses, friction, thermal dynamics, surface tension, and many other effects that are nonlinear and difficult to calculate by hand in a timely manner. For these types of problems, a numerical approximation or simulation using the Finite Element Method (FEM) [2] is best suited. Creating a model can reduce the time and cost of development and help detect unusual or unexpected behavior of the system without testing every possible condition.
Approach
For our project, we needed to create a model to analyze fluid flow around and through different 3D geometries. A better understanding of how these shapes would affect flow properties was crucial for optimizing the outcome. We aimed to design a structure that would convert the fast laminar flow of the Rhine, moving at 3–12 km/h [3], into a slower, more turbulent flow. This would maximize water’s contact with the inner structure of the vessel and give bacteria the opportunity to attack PAH molecules.
We came up with three different designs that we hoped would help achieve this. The first is a maze-like pattern (Figure 1), where our bacteria would be immobilized on the inner walls. The second design is cylindrical, featuring pillars arranged to create obstacles for the flow, breaking it up into turbulence (Figure 2). In this design, the bacteria would again be immobilized on the surface of the cylinders. The third design idea is, to use balls that could be simply placed in a container, creating a sponge-like pattern (Figure 3). The advantage of this design is that the balls, along with the bacteria immobilized on their surfaces, could be easily replaced. However, these were merely ideas, and we needed to develop a model to verify which design would perform best.

Figure 1: Maze structure.
Figure 2: Pillar structure.
Figure 3: Ball structure.
Simulating the models
After creating and preparing the 3D model we needed to simulate their behavior in a flow condition. Would they create the desired speeds but keep enough pressure to keep the flow going? We opted for ANSYS Fluent [4], an industry standard simulating software capable of simulating complex fluent mechanical problems. We created a fluent problem and started to mesh the geometries. This means creating a net of triangles and quadrans across the geometry, that have varying density and will be a single element when applying FEM solving algorithms. This was by far the trickiest part, and we had to go back and tweak the models because ANSYS was not able to mesh some of the geometries. Next, we had to set up our simulation properties and boundary-conditions: What material was used, which temperatures, which speeds and pressures are present at the inlet and outlet. We used the current data from the Federal Institute for Hydrology measured in the river Rhine at Mainz [5]. We opted for a steady state simulation, due to no moving parts in our models and a faster simulation time. The steady state results would show us how our geometry interacts with the flow after enough time has passed and a steady state has settled. This is what we were interested in. For all our models the residuals converged, meaning we would get a reliable result.
Evaluation of results
In our design, the most critical parameters are pressure and velocity. The velocity at the surface contact areas should be slow but still exhibit movement, with values near zero being ideal, to give the bacteria much time to attack the pollutants. Pressure should gradually decrease throughout the geometry to ensure a consistent flow.

Maze Structure

Figure 4 shows the velocity of water flowing through the maze. In the corners, the flow nearly stops, while maintaining moderate speed in the middle of the canal. As we move further along, we encounter more slow-flow regions. The inlet velocity (Inlet on the right side) increases from 6.2 m/s to around 12 m/s on the first corner and remains relatively steady across the geometry. Interestingly, the flow appears to accelerate again, particularly around the corner regions. This can be attributed to a decrease in flow area; in these cases, the very slow-flow regions act like walls, reducing the available flow area and consequently increasing the velocity. This explains the co-existence of both low and high-speed flow in the system.

Figure 4: Velocity of water flowing through maze, inlet is on the right side.

Figure 5 shows the pressure drop across the geometry. In incompressible flow, pressure typically increases as velocity decreases, and vice versa. In our case, we observe a pressure drop in each turn of the geometry, which suggests a flow from inlet to outlet, as fluid usually moves from high to low-pressure regions. Around the corners, there is a decrease in pressure, and at the outlet, we observe a negative pressure, indicating that the structure actively draws fluid back in.

Unfortunately, we see few regions where velocity and pressure have a steep gradient which suggests good mixing of the flow. It implies that the flow remains mostly laminar rather than turbulent. This limits the amount of water in contact with the walls and bacteria.
Figure 5: Pressure across maze design, inlet is on the right side.


Cylindrical Structure

Figure 6 shows the velocity of flow through the cylindrical structure, which looks very promising. There is a continuous, steady flow from inlet to outlet, ensuring that fresh water is transported in and old water is removed. Behind the cylinders, we observe areas of very low velocity, while the regions between the cylinders exhibit higher velocities. This separation of flow zones is desirable, as it creates slow flow areas where pollutants can be absorbed by the bacteria.

Figure 6: Velocity through cylindrical pillar structure, inlet is on the left side.

Overall, this design features more slow-flow regions compared to the maze, suggesting that we are on the right track. Additionally, in front of each cylinder, a stagnation point (or line, in 3D) forms where velocity is low and pressure is high. This indicates that the flow simulation is reliable. Figure 7 displays the pressure drop across the cylindrical design. The cylinders create bands of equal pressure, suggesting minimal mixing occurs within each band. From a high inlet pressure to a low outlet pressure, the continuous flow prevents the reverse flow issues seen in the maze design. However, the contact area in this configuration could be improved.

Figure 7: Pressure across cylindrical design, inlet is on the left side.


Ball Structure

Figure 8 presents a horizontal cut through the ball design. The flow pattern is similar to that of the cylindrical structure, but with significantly more contact between the surface and fluid, which is exactly what we intended. Around the balls, a layer of slower flow is observed, which is ideal for allowing bacteria to absorb substances more effectively. This design incorporates the advantages of the cylindrical model, but with increased surface area and improved slow-flow characteristics. Combined with the ease of production and replacement of the balls, this concept clearly outperforms the other designs. For these reasons, we proceeded with the ball design in our project.

Figure 8: Cut through ball design, velocity of flow-through is shown.


Room for improvement
The simulation helped us determine the best Design for our project. For even better and more realistic results a couple of things need to be done:
  • Simulate pollution and particles like dirt and wood in the water
  • Simulate the transient behavior
  • Simulate the uneven surface of our geometry when it is 3D printed
  • Simulate the wear and stress of the material to improve the design and increase its lifetime
All these aspects would help to improve our design further but are very time-consuming and not doable in the context of this project. Still further investigation might show interesting results.

Metabolic modeling

Preperation
Goal: Pseudomonas vancouverensis (DSM8368) [6] is known for its ability to degrade naphthalene and phenanthrene, two common polycyclic aromatic hydrocarbons (PAHs) [7]. However, it lacks the metabolic pathway necessary to break down pyrene, a more complex and toxic PAH. In this project, we focus on engineering P. vancouverensis to degrade pyrene by introducing the pyrene degradation pathway, and compare its performance to Pseudomonas putida KT2440 [8], which serves as a control strain. Using a computational approach helps us to understand the non-model organism and provides us valuable insight for the design and optimization of our PAH's degradation pathway.

Tools

We used COBRApy, a python package for Constraint-Based Reconstruction and Analysis (COBRA). This allowed us to simulate metabolic reactions and predict the fluxes under different environmental conditions [9].

Assumptions

  • Steady state assumption: By assuming a constant rate of production and consumption of the metabolites the calculation of the differential equation for the flux balance analysis can be simplified.
  • Dynamic pyrene uptake: The rate of pyrene uptake is modelled as a function of its concentration in the medium, this is done by applying Michaelis Menten kinetics. The approach is that the uptake of pyrene is dependent on the concentration. If concentration is low, the uptake increases linearly whereas at higher concentration the rate approaches a maximum value, because of the chassis' limited uptake capacity.
  • Unlimited Enzyme Capacity: Enzyme capacity limitations or regulatory mechanisms were not included in this model. This assumption simplifies the focus to the strain's metabolic potential.
  • Nutrient-Rich Medium: The medium contains enough glucose and pyrene, ensuring that neither carbon nor energy sources are limiting factors during the simulations.

Mathematical Background

The mathematic representation of a metabolic reaction for a genomic-scale metabolic model can be seen in Figure 9. It consists of a stochiometrix matrix which represents all the reaction in every column (n reactions) and the compounds in the row (m compounds). It is negative when a compound is consumend, positive when a compound is produced and zero if it is not involved in the reaction. The reaction fluxes are represented in a vector v which consists of all the fluxes. Applying and solving a kinetic model must contain the rate laws for each reaction which gives us a flux distribution lying at any point in the solution space. With the definition of some constrainment such as the steady state assumption, we assume a constant rate of production and consumption of compounds (Sv = 0). This limits the solution space. Additionally every reaction has an upper and lower flux constrain which adds an individual constrain for each reaction flux. After the calculation of different flux distribution is done, an optimization for a defined solution can be made. This is usually the biomass production rate [10].

Figure 9: The matrixes and vectors formulation of a genome-scale metabolic model. It consists of a stochiometric matrix S with n reactions and m compounds, a flux vector v and for the solution calculation the fomulation of the kinetics


For a dynamic flux balance analysis a complete steady state assumption is ignored, rather than that the variability of consumption or production rate is dynamic. In COBRApy this was achieved through a simple Michael Menten kinetics (Figure 10), where the assumption that the enzyme concentation or compound concentration is higher were made, such that the rate of the enzyme-substrate complex is zero [11]. The program COBRApy uses a static optimization approach (SOA) [12] meaning that snapshot were taken at different time points (here around every 23 min : 19h simulation with 50 steps if equal distributed).
Figure 10: Michaelis-menten kinetic with the metabolite concentration [A], the Michaelis-constant Km and the max. velocity of the reaction, vmax
1. Step: Testing the model
Goal: Load the Model, add the PAHs' degradation pathway, and simulate a successful pyrene breakdown.

Model loading

Because, to our knowledge, no metabolic model of P. vancouverensis exists yet, the first step is to load the COBRA metabolic model of P. putida using COBRApy. This model contains all known reactions involved in the strain’s metabolism, including energy production and biomass formation though lacking P. vancouverensis’ PAH-metabolism. Therefor we use it as basis for our simulation. The model is loaded from JSON-file: iJN1463.json.

Reaction Addition

Because the base model of iJN1463 does not represent the natural capabilities of P. vancouverensis, which includes complete naphthalene and phenanthrene degradation pathways, we added these reactions and this project’s pyrene degradation pathway. These pathways eventually produce catechol, which is further processed through the nitrobenzene degradation pathway already present in the model. This process ultimately leads to the production of pyruvate [13]. In this way the PAHs can be completely used as carbon sources.


Table 1
iJN1463 (before)
iJN1463 (after)
Genes
1462
1462
Metabolites
2153
2171
Reaction
2927
2946
Table 1: Number of genes, metabolites and reactions for the model iJN1463 before and after the addition of PAH degradation.

The test growth in Figure 11 shows that pyrene could be utilized by the model after the reaction was added, as the flux rate of the pathway's reactions were greater than 0.

Figure 11: Test grow for biomass optimization of our model iJN1463 with pyrene in the medium. The blue curve (left y-axis) represents the biomass concentration over time and the red curve (y-axis) shows the corresponding decline in pyrene concentration in the medium.

2. Iteration: Applying and comparing with control organism
Goal: Compare the results of biomass formation and pathway flux in P. vancouverensis with and without the pyrene degradation pathway, and with P. putida as a control which can not metabolize PAHs.

Figure 12: Biomass formation of P. vancouverensis DSM8368 on the left and P. putida KT2440 on the right with and without pyrene degradation plasmid. The lined curve represents the unmodified strain and the solid lines represents the transformed strain.

As shown in Figure 12, P. vancouverensis and P. putida with the additional pathway shows higher biomass production compared to the unmodified. This increase in biomass is attributed to the effective utilization of pyrene as a carbon source by the introduced pathway. The modified strain can metabolize pyrene, leading to enhanced growth and thus demonstrating the successful incorporation of the degradation pathway.
3. Addition: Defined medium for experimental reproducibility
Continuing with one strain for degradation optimization: P. vancouverensis will be further analyzed. To ensure reproducibility later in the experiments M9 medium will be used as it is well defined in contrast to LB medium. Since fluxes have the unit mmol/gDW*h (concentration per gram dry weight of cell and hour) for the exchange bound the approximation will be made that the medium supply is 1 mol of our substance x every 1h to 1g of bacteria as suggested in [14]. Water/OH- uptake was set to -100 to 100 mmol/gDW*h and oxygen to -100 - 20 mmol/gDW*h as suggested in the paper where the breakdown of plastic waste was improved [15]. The M9 growth media will be coded, with the recipe from Thermofisher that we ordered (A1374401) [16].

Table 2
compound
mass in 1L in g
Molar mass [g/mol]
Mmol/h*gDW
Sodium phosphate (dibasic) heptahydrate
12.8
268.07
47.7
Monopotassium phosphate
3
136.09
22
Sodium chloride
0.5
58.44
8.6
Ammonium chloride
1
53.49
18.7
D-Glucose
4
180.16
22.2
Magnesium sulphate
0.241
120.37
2
Calcium chloride
0.111
110.98
0.1
Table 2: M9 compound according to Thermofisher for 1L

Inspired by the approach from the Aalto iGEM Team in 2023 the ions of the salts will be added separately to the medium [17]. This allows us to adjust the bounds of the exchange reactions, ensuring proper simulation of nutrient availability. Trace elements were added because no growth in M9 was calculated. By analyzing the reduced cost, we identified the most limiting reaction for our biomass/growth rate. For instance, increasing iron supply from 0.007 to 0.009 mmol/L resulted in a significant increase in growth rate from 0.48 1/h to 0.613 1/h. This strategy helped us to optimize the nutrient supply.


Figure 13: 3D plot showing the relationship between the growth rate (1/h in the Z-axis and heatmap) and the uptake rates of pyrene and glucose for the transformed P. vancouverensis DSM8368 in M9 medium with trace elements. The X-axis represents the glucose uptake rate and The Y-axis shows the pyrene uptake rate.

In Figure 13 it seems in M9 Medium the pyrene uptake rate and thus the concentration has no impact on the growth rate. This is unexpected, as it contradicts the findings above. The solution to this problem is still being worked on.
Future improvements
For future models an implementation of enzyme cost would provide more realistic outcomes by accounting for energy and resource demands of enzyme production and activity. Furthermore, incorporating gene knockout experiments into the model would enhance our understanding of the genetic determinants of pyrene degradation. By systematically removing key genes and observing the resulting metabolic shifts we could try for improvements of the pyrene enzyme fluxes such that only the most efficient and essential enzymes are produced. In combination with the enzyme cost implementation this could give rise to a better and more realistic design.

References

[1] A. T. Lawal, "Polycyclic aromatic hydrocarbons. A review," Cogent Environmental Science, vol. 3, no. 1, p. 1339841, 2017, doi: 10.1080/23311843.2017.1339841#d1e147.
[2] K.-J. Bathe, "Finite Element Method," in Wiley Encyclopedia of Computer Science and Engineering, B. W. Wah, Ed.: Wiley, 2007, pp. 1–12.
[3] Jörg Krause, Baden im Rhein kann lebensgefährlich sein. Mindestens 378 Menschen in 2023 ertrunken. [Online]. Available: https://polizei.nrw/Lebensgefahr%20im%20Rhein (accessed: Sep. 25 2024).
[4] J. E. Matsson, An Introduction to Ansys Fluent. Mission (Kans.): SDC Publications, op. 2023.
[5] Landesamt für Umwelt Rheinland-Pfalz (LfU), Gütemessstellen im Rheingebiet: Mainz-Wiesbaden, Rhein. [Online]. Available: https://undine.bafg.de/rhein/guetemessstellen/rhein_mst_mainz_wiesbaden.html (accessed: Sep. 1 2024).
[6] W. W. Mohn, A. E. Wilson, P. Bicho, and E. R. Moore, "Physiological and phylogenetic diversity of bacteria growing on resin acids," Systematic and applied microbiology, vol. 22, no. 1, pp. 68–78, 1999, doi: 10.1016/S0723-2020(99)80029-0.
[7] Y. Yang, R. F. Chen, and M. P. Shiaris, "Metabolism of naphthalene, fluorene, and phenanthrene: preliminary characterization of a cloned gene cluster from Pseudomonas putida NCIB 9816," Journal of bacteriology, vol. 176, no. 8, pp. 2158–2164, 1994, doi: 10.1128/jb.176.8.2158-2164.1994.
[8] M. Bagdasarian et al., "Specific-purpose plasmid cloning vectors. II. Broad host range, high copy number, RSF1010-derived vectors, and a host-vector system for gene cloning in Pseudomonas," Gene, vol. 16, 1-3, pp. 237–247, 1981, doi: 10.1016/0378-1119(81)90080-9.
[9] J. Schellenberger et al., "Quantitative prediction of cellular metabolism with constraint-based models: the COBRA Toolbox v2.0," Nat Protoc, vol. 6, no. 9, pp. 1290–1307, 2011, doi: 10.1038/nprot.2011.308.
[10]J. D. Orth, I. Thiele, and B. Ø. Palsson, "What is flux balance analysis?," Nature biotechnology, vol. 28, no. 3, pp. 245–248, 2010, doi: 10.1038/nbt.1614.
[11]L. Michaelis, M. L. Menten, K. A. Johnson, and R. S. Goody, "The original Michaelis constant: translation of the 1913 Michaelis-Menten paper," Biochemistry, vol. 50, no. 39, pp. 8264–8269, 2011, doi: 10.1021/bi201284u.
[12] 15. Dynamic Flux Balance Analysis (dFBA) in COBRApy — cobra 0.27.0 documentation [Online]. Available: https://cobrapy.readthedocs.io/en/latest/dfba.html (accessed: Oct. 1 2024).
[13] Nitrobenzene Degradation Pathway. [Online]. Available: http://eawag-bbd.ethz.ch/nb/nb_map.html (accessed: Sep. 14 2024).
[14] 12. Growth media — cobra 0.27.0 documentation. [Online]. Available: https://cobrapy.readthedocs.io/en/latest/media.html (accessed: Sep. 18 2024).
[15] Leah A Lewis, Matthew A Perisin, and Alexander V Tobias, "Metabolic modeling of Pseudomonas putida to understand and improve the breakdown of plastic waste," 2020.
[16] M9 Minimal Salts (2x). [Online]. Available: https://www.thermofisher.com/order/catalog/product/A1374401?SID=srch-srp-A1374401 (accessed: Sep. 20 2024).
[17] metabolic_modeling/modeling.ipynb · main · 2023 Competition / Software Tools / Aalto-Helsinki · GitLab. [Online]. Available: https://gitlab.igem.org/2023/software-tools/aalto-helsinki/-/blob/main/metabolic_modeling/modeling.ipynb?ref_type=heads (accessed: Sep. 18 2024).