Model | OhioState - iGEM 2022

OhioState iGEM 2021’s team made a model that depicted how the population of bacteria changed after the addition of bacteriophages, using a paper written by Dr. Lenski¹. The paper contained a system of ODEs that could be used to describe the relationship between bacteria and bacteriophages. However, due to complications with two parameters within the model, the model was forced to be created using the steady state approximation. The two parameters S’ and P’, the concentrations of free phage and uninfected bacteria at time t-τ¹, were the parameters causing the complications.

While OhioState iGEM 2021’s model and findings were valid, they only applied to situations where the system was at or approximately at steady state. Due to this limitation, the utility of the model was largely reduced since it could only be applied in specific situations.To widen the scope of the model, new methods of solving the system of ODEs were tested to create a model that could be used in both steady and unsteady state systems.

After re-analyzing OhioState iGEM 2021’s steady state model, different possible solutions were determined to address the S’ and P’ parameters. The team then emailed Dr. Lenski to see what methods could be used to address the parameters. Dr. Lenski provided two possible routes the team could go to addressing the parameters:

Using a smaller time step in the simulation
Utilizing the Runge-Kutta method to create the model

In OhioState iGEM 2021’s model, an incredibly small time step was used. When testing to see if decreasing the time step further addressed the complications, the model displayed the same issues as OhioState iGEM 2021’s simulation and took longer to run. Therefore, the team moved to utilize the Runge-Kutta method to develop the model. The Runge method integrates ODEs by using a trial step at the midpoint of the half step to cancel out the error². The method uses the additional parameters (k) as correction factors, improving the model’s output. The model was developed using the RK4 method, meaning 4 parameters were used to determine the slope between two points.

Figure 1. RK4 Fundamental Equations

For the first iteration, the model was created with the assumption that S’=S and P’=P. The parameters S and P are the concentrations of free phage and uninfected bacteria at time t. This was done to see whether the Runge-Kutta method could eliminate the complications without significant alterations to OhioState iGEM 2021’s code. Using OhioState iGEM 2021’s model, setting S’=S and P’=P led to concentrations of free phage and uninfected bacteria either not changing with time or going to positive/negative infinity. The RK4 developed model produced a more reasonable model with the same assumption:

Figure 2. Plots Generated during first iteration, S'=S and P'=P

Although the model created using the RK4 method produced more reasonable results than the OhioState iGEM 2021 model when using the S’=S and P=P assumption, the model did not properly display the changes in bacteriophage and bacteria populations over time. There was one main potential issue that was determined when analyzing the produced plot. The first being the slope at the beginning of the simulation for both uninfected bacteria. When analyzing the zoomed in plot, the concentrations of both the free phage and uninfected bacteria are both decreasing. The plot implies that phage is attaching to bacteria, which makes sense as to why the concentration of uninfected bacteria is decreasing. However, the rate of which they’re decreasing doesn't match, which implies that the population of uninfected bacteria is decreasing due to other factors. Overall, the general plot displayed the predator/prey relationship between the concentrations; however, we weren’t sure whether that initial decrease in both concentrations was correct. Therefore, we ran another iteration with a slight alteration of how S’ and P’ were defined.

For the final iteration, a simulation was created that would assign S’ and P’ with their associated value based on the time t-τ. τ represents the period after infection in which beta phages are produced, while t represents the actual time. τ was set to 0.5 hours, and it was tracked as the simulation progressed. When t<=τ, S’ and P’ were set to 0 within the simulation. When t>τ, S’ = S(t-τ) and P’ = P(t-τ). Running the simulation produced the model shown in Figure 3.

Figure 3. Plots Generated during final iteration, S'=S(t-τ) and P'=P(t-τ)

The plots generated from the updated values still portrayed the same predator-prey relationship as the previous iteration. However, here the difference can be seen in the zoomed in plot. It displays that as the concentration of phage increases, the concentration of bacteria also decreases. This seems more logical than the previous iteration’s assumption that the concentration for both would decrease. We can see that bacteria’s concentration begins to increase before the phage concentration flattens out. This could be due to the bacteria beginning to split themselves as the phages are still lysing the infected bacteria. Once the concentration of uninfected bacteria reaches its apex, the concentration of free phage begins to drop as it attaches itself to uninfected phage. A similar cycle will occur as time passes. This iteration is seen as an improvement of the previous iteration, the plots seem a bit more logical than the previous iteration’s plots. The population model was constructed in MATLAB. A commented pdf file containing the assumptions, paramaters, ODEs, and code can be downloaded here.

References

(1) - Lenski, Richard E. "Dynamics of interactions between bacteria and virulent bacteriophage." Advances in microbial ecology (1988): 1-44, https://lenski.mmg.msu.edu/lenski/pdf/1988,%20Adv%20Microb%20Eco,%20Lenski.pdf

(2) - Weisstein, Eric W. "Runge-Kutta Method." From MathWorld--A Wolfram Web Resource. https://mathworld.wolfram.com/Runge-KuttaMethod.html

To successfully integrate bacteriophages into the medical field as a measure against sepsis, a method of mass producing the phages is necessary. After researching common processes, the most suitable method for bacteriophage production was determined to be via bioreactors. A bioreactor is a unit used to grow organisms under controlled conditions¹. The team determined early in the process that the process would need other units to develop the phages. Given the nature of bacteriophages, specifically the fact that they are unable to reproduce by themselves, a minimum of one unit would be needed to grow the bacteria and another that would introduce the phage to the bacteria.

For reference, the team used Escherichia coli (E. coli) as the bacteria that phage would infect. Most strains of E. coli are harmless², which would, in the event of some sort of unit failure, reduce danger for operators. Although most are harmless, thorough research would need to be done into the chosen strand to ensure the safety of individuals operating the system.

Before committing to building the bioreactor, the team analyzed past team’s pages to determine whether the build would be viable with the given time. After looking through different pages it was determined that building it would not be possible given the team’s resources. Therefore, the team decided to model an aspect of the bioreactors. The team decided to generate a dynamic model for the bioreactors. The goal of the dynamic model was that it would be the first step towards developing a control system that could be used in industry. A process control system maintains a unit at its ideal conditions³. For example, if the temperature in the bioreactor increases, the thermal jacket (used to control the temperature) of the unit can be adjusted to remove heat from the bioreactor. With a developed control system, the bioreactors would be able to operate at optimal conditions while protecting individuals operating the systems.

Before developing the dynamic model, the team needed to determine a specific objective that needed to be achieved. The team decided to create a dynamic model for a single bioreactor, with the modeled bioreactor being the bioreactor that would introduce phage to the bacteria. The unit that would develop the E. coli will not be included, the model will assume that the incoming stream of bacteria will be from an upstream unit. Our goal is to create dynamic models that control the leaving stream’s flow rate and temperature. The following assumptions were made when creating the model:

The bioreactor is well mixed
Pressure in the bioreactor is maintained at a constant pressure P₁
Bioreactor pH is maintained at a constant pH
Heat capacities are independent of temperature
Negligible kinetic and potential energy, internal energy dominates the system
Boundary work and work due to the stirrer can be neglected
The internal energy is approximately equal to enthalpy
Bacteria and phage interation occuring in an dominating inert solvent, assuming constant density

After defining the assumptions, a schematic was drawn as a reference to the calculations:

Figure 1. Bioreactor Schematic

Within the bioreactor, the interaction between the phages and bacteria are treated as a reaction. Unfortunately, the kinetics could not be tested in the lab, so a baseline assumption will be made regarding the reaction. We'll assume that one phage and one bacteria will form 1 phage. We'll also assume that this reaction is elementary, since the phage and bacteria interacting is common (when treated as bimolecular reaction). Assuming an elementary reaction, the reaction rate will be a function of stoichiometric coefficients.

Figure 2. Left image is the reaction equation. Right image displays the reaction rate.

To start, a general mole balance is created for the batch system. This will be the baseline for the following calculations.

Figure 3. General Mole Balance

In the general mole balance equation:

V represents the volume of the bioreactor
Q_feed represents the volumentric flow rate of the inlet streams
Q represents the volumetric flow rate of the outlet streams
C_phage represents the concentration of phage in the stream or bioreactor
C_{E. coli} represents the concentration of E. coli in the bioreactor
R_phage represents the production rate of phage within the bioreactor.

At this point, there are simplifications that can be made of the system. Since the bioreactor is running a batch process, the volumetric flow rates are equal to 0, as nothing is entering or leaving the system as the phage is reproducing. This simplification produces the following equation:

Figure 4. Simplified Mole Balance

The simplified mole balance displays that the change in phage concentration in the bioreactor is dependent on the production rate of phage. However, a further simplification can be made with assumption 8, the constant density assumption. The validity of the assumption could be argued against, since phage is being produced from bacteria, the concentration of phage would increase, and the concentration of bacteria would decrease. Therefore, there would most likely be a change in density within the bioreactor, with the extent of change being dependent on the changes of both components. However, the assumption being made inherently is that the change in density of the reaction is negligible, since there is an inert solvent that dominates and sets the density of the solution. Since, the constant density assumption is being used, the volume within the bioreactor would be constant. The V term within the derivative can be removed and canceled out with the volume term on the right side of the equation, giving the following equation:

Figure 5. Simplified Mole Balance with the constant density assumption

The R_phage term is equal to the reaction rate of the species, which in this case is phage. Since the overall production of phage is 1, one phage begin consumed to make two, R_phage is equal to the reaction rate.

Figure 6. Mole Balance with Production Rate Substitution

We can integrate this equation to find an equation that determines C_phage at any time.

Figure 7. Phage Concentration equation

Here we have our final phage concentration equation, and we can see that it is also dependent on the concentration of E. coli. Figuring out how long the batch process should run will be important to operate the bioreactor. To understand the optimum time, we also need to understand how the concentration of E. coli changes over time. The derivation is similar to the phage concentration derivation. The only difference is the definition of the production rate for E. coli. In this case, we can see that the overall production rate of E. coli is 1 according to the chemical equation. This means that the production rate for E. coli is equal to the following equation:

Figure 8. Production Rate for E. coli

After integrating the a final equation for E. coli is found:

Figure 9. E. coli concentration equation

Using the two derived equations a model can be created to see how these concentrations change over time. That model would allow us to select an ideal run time for the bioreactor, with the goal of maximizing phage production. More models can be created from here to control other factors such as temperature and pressure. With all aspects modeled, a bioreactor can be made to test whether these models are valid for use, or if adjustments need to be made to reach the desired phage output.

References

(1) - “What Is a Bioreactor?” Bioreactors.net https://www.bioreactors.net/what-is-a-bioreactor

(2) - Sepsis Alliance. Sepsis and Intestinal E. coli Infections. 2022. https://www.sepsis.org/sepsisand/intestinal-e-coli-infections/

(2) - A Short Introduction to Process Dynamics and Control. http://www.users.abo.fi/khaggblo/PDC/PDC%20intro%20-%20longer.pdf

Our team wanted to provide a visual depiction of the interaction between bacteria and bacteriophages, specifically the benefits of using a phage cocktail.
To model the benefits of a phage cocktail, a cellular automata method was created (one really famous example of cellular automata is Conway’s Game of Life). For this method, there are discrete cells that can have different states. These states are determined by rules and these rules are user determined. Another important aspect of cellular automata is the neighborhood. This neighborhood describes the area around each cell that that cell can impact.

Figure 1. Depiction of Two Common Neighborhoods for Cellular Automata, von Neumann and Moore, and How the Radius Parameter Affects Their Shape

The model simulated a very simplified human body by treating the system as a 2D plane. Many other assumptions were made to help simplify the model and allow it to focus on the difference between using 1 or 2 phage types. To better understand our model, the final model assumptions are listed below:

Unlimited resources available to both the phage and bacteria
All bacteria divide equally
The initial phage dosed into the system surround their corresponding bacteria right away (type 1 match up and type 2 match up)
Phage type 1 can only infect bacteria type 1 and phage type 2 can only infect bacteria type 2
Bacteria type 1 uses less resources to replicate and so will always beat out bacteria type 2 if there is room to grow
All types of phage can surround all types of bacteria at the same time. There is no space limiting dimension, meaning that bacteria 1 can be surrounded by phage 1 and 2 at the same time and the same for bacteria 2
When a bacteria is surrounded by any phage, it is under pressure and so cannot divide
Phages do not perish or disappear so a type 1 phage surrounding a type 2 bacteria will stay there forever and there will always be a surrounded bacteria
When a phage lyses the bacteria only phage of the same type that lysed the bacteria (and same type as bacteria based on assumption 4) are released an any of the same phage type on that site are fully replaced with the progeny phage. For example, if there were ten phage type 1 surrounding bacteria type 1 and they infect, when it is lysed, the bacteria releases progeny # of phage type 1 (the amount (ten) that were previously there do NOT matter and don't get added to the total)

On the 2D plane, there were eleven possible states each cell could have. To create a visual output, each state was linked to a certain color and both state and color descriptions are described below.

State 0: dead bacterial cell or bare space - white

State 1: healthy bacteria of type 1 - green
State 2: healthy bacteria of type 2 - orange

State 3: surrounded bacteria type 1 by phage type 1 - cyan
State 4: surrounded bacteria type 1 by phage type 2 - purple
State 5: surrounded bacteria type 2 by phage type 1 - yellow
State 6: surrounded bacteria type 2 by phage type 2 - pink

State 7: infected bacteria type 1 - blue
State 8: infected bacteria type 2 - red

State 9: surrounded bacteria type 1 by both phage types - dark green
State 10: surrounded bacteria type 2 by both phage types - dark orange

Figure 2. Example Output Showing Most of the Possible States

These states are determined by our model’s rules. These rules seek to parallel actual interactions between bacteria and phage.

Rule 1: A healthy bacteria has p = #phage * search chance of being surrounded for each phage type. A portion of these phage (based on search) are then put into the same location of as the "found" bacteria
Rule 2: Dead/empty cells have p = growth chance of becoming an alive bacteria if surrounded by exactly 1 or 2 healthy bacteria. These must be the same type (the surrounding type will dictate what type of bacteria is grown and type 1 will beat out type 2). If bacteria are surrounded they are under pressure and will not replicate
Rule 3: Surrounded bacteria have p = #phage * lethality chance of being infected by one phage(only the number of phage of the same type are counted into this)
Rule 4: After 1 generation, infected bacteria die and release a fraction of progeny # of phages to the cell (phage released are of the same type as the bacteria dying and all of the other phage type are assumed inactive, thus after lysing the only phage in that cell will be a fraction of progeny # of phage of the same type as the bacteria that died)

As seen in the above rules, there are some parameters that change probabilities and the ultimate outcome of the rules. These parameters are used to tune the model to investigate different scenarios. For our model, most of these parameters are user-inputted, the user of the code can choose these values to their liking.

density: Describes the approximate percentage of sites occupied by bacteria at time = 0. Must be a value from 0 to 1 inclusive.
neighborhood: Defines the type of neighborhood to be used to see how many bacteria are around a site and how many phage are around each bacteria. Only 'vonNeumann' or 'Moore' are accepted values. The default is 'vonNeumann'.
radius: Defines how big the specified neighborhood is. Can either be 1 or 2. The default is 2.
nGen: The maximum generations the code will run for if possible. Must be an integer. If None, the animation stops once no more bacteria can become infected. The default is None.
pbc: If True, periodic boundary conditions are used. If False, dead zone boundaries are used. The default is False.
grid: If True, grid lines are shown on the plot. If False, grid lines are not displayed. The default is True.
gridSize: Defines the space with which to create the bacteria on. Essentially the grid dimensions. Must be a tuple of 2 ints or floats. The first value is xMax and the second yMax, so the grid will run from 0 to xMax and from 0 to yMax. The default is (50,50).
progeny: The number of phage released once a bacteria is lysed (killed). Must be an int or float. The default is 100.
growth: Probability of bacteria replicating. Must be a value from 0 to 1 inclusive. So a value of 0.5 would mean there is a 50% of a bacteria replicating. The default is 0.1.
search: Probability of phage finding and navigating to a bacteria. Also dictates the portion of surrounding phage that then move to surround the bacteriaMust be a value from 0 to 1 inclusive. So a value of 0.5 would mean that each phage have a 50% of finding a bacteria if it is in their neighborhood and that 50% of the phage in the bacteria's neighborhood surround it if the bacteria becomes surrounded. The default is 0.5.
lethality: Probability of phage killing a bacteria. Must be a value from 0 to 1 inclusive. So a value of 0.5 would mean that each phage, once surrounding a bacteria, have a 50% of then infecting that bacteria. The default is 0.5.
nbact: The number of bacteria types to be present in the model. Must be either 1 or 2. The default is 2.
nphage: The number of phage types to be present in the model. Must be either 1 or 2. The default is 2.
amountPhage1: The initial number of phage type 1 put into the model. The default is 100.
amountPhage2: The initial number of phage type 2 put into the model. The default is 100.

Having made assumptions, defined states, come up with rules, and tuned those rules with user defined parameters, the model was ready to be run. Its first step is to dose a set amount of phage at the center of the 2D space and then each rule is put into place and new states determined. Each time a rule is implemented and new states determined, a new plot is generated. This time step is termed generations and when run, the plot will be continually updated like a time lapse. The output plot has three different graphs. The left and center graphs are phage heatmaps, they show where phage are concentrated. A red color signifies the highest number of phage, descending to warmer colors orange and yellow, down to cooler colors until dark blue, the lowest number of phage. The right graph is the bacteria graph, showing the different states of the bacteria as signified by their color.

Figure 3. Example of a Full Plot Output Including Both Phage Heatmaps and Actual 2D Plane

This model can now be used to show the difference between a phage cocktail and just a normal phage. If we set the types of phage to 1 and types of bacteria to 2, we get the output below.

Figure 4. Resulting Output When Only 1 Phage Type is Introduced into the System with 2 Bacteria Types

Now if we pretend to administer a phage cocktail and have two types of phage and two types of bacteria, we get a much cleaner graph, showing a larger absence of bacteria.

Figure 5. Resulting Output When a Phage Cocktail (2 Phage Types) is Introduced into the System with 2 Bacteria Types

The goal of this model was to show differences in bacterial concentrations when comparing a phage cocktail to just phage. Despite many assumptions, our model does indeed show that with the 2 phage type system resulting in a drastic decrease in bacteria filled cells. That's how important phage cocktails are. If you have a single bacteria type, they are bound to mutate, essentially resulting in two total types of bacteria. Using a single phage type could prove ineffective if a receptor mutated and the phage can no longer infect. However, if a phage cocktail is used, two phage types in this example, there is an increased chance that at least one of the phage can still infect the mutated bacteria and insert its genetic material.

This model was constructed in Python, so you can try it out yourself. You can experiment with parameters or even add on to the code itself. A commentted pdf file containing the Python code can be downloaded here.

Overview

References

References

References