Team:TAU Israel/Model
type="text/css">
"What distinguishes a mathematical model from, say, a poem, a song, a portrait or any other kind of "model", is that the mathematical model is an image or picture of reality painted with logical symbols instead of with words, sounds or watercolors."
John L. Casti
Motivation
The role of the models in biological experiments (and in other fields in science) is to help us understand the experiments' outcomes, explain the system and predict different results based on changes we make in the system.
In our project, we created a mathematical model that predicts the pyocins-bacteria interactions and dynamics in 2D. As part of the project's proof of concept, the model will help us track and optimize the interaction and diffusion of pyocins and bacteria, and bacterial death as a result of these interactions in various 'real-life' scenarios. The model, together with a relevant algorithm, is also part of our suggested 'package', which includes not only designed pyocins, but also tools for optimization the efficiency and cost of their usage.
In the Diffusion-interaction model section more information about our mathematical model will be provided. In addition, in the Bioinformatics models section you can read more about the project's part improvement and the discovery of a new tail process.
Diffusion-Interaction Model
Background
We decided to model the amount of bacterial and pyocin populations over time and space in a medium (e.g agar), while taking into consideration the diffusion and the effect of each population on the other.
In the diffusion factor, we used Fick's second law.
Fick's second law is a partial differential equation, which describes density fluctuations in a compound undergoing diffusion[1]:
$$ \frac{\partial c}{\partial t} = D \frac{\partial^2 c}{\partial x^2} $$
While examining the different factors that effect both populations (bacteria and pyocins), besides diffusion, we decided to use a set of similar problems from the dynamics of ecological systems world - "Lotka-Volttera" equations, also known as the "prey-predator" equations[2][3]:
Where x is the number of prey, y is the number of predator, $ \alpha $ is the growth factor of the prey, $ \beta $ is the predation rate coefficient, $ \delta $ is the reproduction rate of predators per 1 prey eaten, and $ \gamma $ is the predator mortality rate. In our similar problem, the pyocins are the predators and the bacteria are the prey.
The predator mortality rate is equivalent to the pyocins (a protein) degradation rate.
The pyocin-bacteria interaction can be described as a semi prey-predator interaction, due to the fact that the pyocin's population doesn't increase after the binding moment - the pyocin can't be re-used to bind another bacterium, and therefore the population decreases in a manner that is relevant to the interaction we are focusing on.
Therefore, the "Lotka-Volterra" equations were altered to a more suitable set of equations:
Furthermore, because of the resemblance of pyocins' tails to bacteriophage tail structures, especially having the same nature of interaction, we focused on similar studies that investigated the interaction of phages with their target and the effect on both populations due to this interaction[4][5].
BackModel's parameters
It is important to mention that the equations described are generic and can describe any set of parameters and initial conditions (e.g. diffusion coefficients in different parts of the environment, initial number of bacteria and pyocins, growth rates, predation rate coefficient, etc.) defined by the user. Specifically, the values of the model parameters can be altered with measurements to suit specific environments[5] or can be found in the literature[6]. The user can perform sensitivity analysis of the simulation results by using various sets of parameters.
1D Model
After evaluating the different factors that affect the pyocins and bacterial growth, we used Fick's second law and the growth-related parameters to create a 1D diffusion model.
Each equation describes a different population, taking into consideration the second population's amount.
Initial conditions
For bacteria population on interval$ \quad x\in [0,L]$ with initial condition: $$ B(x,0) = B_0, \quad x\in [0,L] $$ And Dirichlet boundary condition: $$ B(0,t)=B(L,t)=B_0, t>0 $$ For pyocin population on interval$ \quad x\in [0,L]$ with initial condition: $$ P(x,0) = f(x), \quad x\in [0,L] $$ And Dirichlet boundary condition: $$ P(0,t)=P(L,t)=0, t>0 $$
Pyocins equation
$$ \frac{\partial P}{\partial t} = D_p \frac{\partial^2 P(x,t)}{\partial x^2} -m P(x,t) - \alpha u(t-d) B(x,t-d) P(x,t-d) $$Bacteria equation
$$ \frac{\partial B}{\partial t} = D_B \frac{\partial^2 B(x,t)}{\partial x^2} +r B(x,t)\left(1-\frac{B(x,t)}{K}\right) - \alpha u(t-d) B(x,t-d) P(x,t-d) $$Where P - pyocins amount, B - bacteria amount, $D_B$ – bacterial diffusion coefficient $[m^2/𝑠𝑒𝑐]$, r - bacterial growth [1/𝑠𝑒𝑐], K - bacterial carrying capacity [amount], $\alpha$ – a parameter that represents the nature of the pyocin-bacteria connection and the killing efficiency, $D_P$ – pyocin diffusion coefficient $[m^2/𝑠𝑒𝑐]$, m – pyocin degradation rate [1/𝑠𝑒𝑐], d – the delay time between the moment of pyocin binding to a bacteria until bacterial death [𝑠𝑒𝑐].
The pyocins are not usable after binding. Therefore, the part in the equations that represents the decline in population because of the binding is equal in both equations.
These equations are not solvable in analytic ways and therefore were converted into discrete equations. Consider the pyocins and bacteria equations with the initial conditions. The discrete equation is according to Forward in Time, Central in Space derivative numerical method (FTCS)[7]:
Bacteria discrete equation
$$ \frac{B[x,n+n_0]-B[x,n]}{n_0} = D_B \frac{B[x+h_0,n]-2B[x,n]+B[x-h_0,n]}{h_0^2} $$ $$ +rB[x,n] (1-\frac{B[x,n]}{K})- \alpha u[n-d] B[x,n-d] P[x,n-d] $$Pyocins discrete equation
$$ \frac{P[x,n+n_0]-P[x,n]}{n_0} = D_P \frac{P[x+h_0,n]-2P[x,n]+P[x-h_0,n]}{h_0^2} $$ $$ -m P[x,n] - \alpha u[n-d] B[x,n-d] P[x,n-d] $$Where u - step function, $h_0$ - step in space, $n_0$ - step in time.
Back1D simulation and the model's sensitivity to parameters
To ensure that the constructed model fits the described pyocins-bacteria interaction, the model's outcome was observed over time using a 1D model simulation in MATLAB, based on different initial conditions. In the figure below, we can see a heat map that portrays the amount of pyocins and bacteria for a period of one second, in time steps of 10 msec. Each row is a time step, while the row itself is our 1D space. The top row is the starting point, time t=0 [sec]. With the progress of every row downwards - we advance in a time step until we reach our final time step of time t=1 [sec] in the bottom row. The right heat map shows the amount of bacteria and the left heat map shows the amount of pyocins. As can be seen, at time t=0 all the pyocins start from a single point, in the middle of the 1D space (left heat map), while there is an even spread of bacteria in the right heat map. The bacteria continue to grow until we reach the delay time (20 minutes)[8], in which the interaction starts to affect both populations, resulting in the populations' decrease. The more we advance in time, the more the pyocins diffuse in space and kill more bacteria.
As an additional tool to validate the model, we changed the model's parameters and checked the effect of the change on the time it takes for the pyocins to kill all the bacteria. The increase or decrease of each parameter can shorten or lengthen the kill time.
The left graph shows the relationship between the total bacteria elimination time and the pyocins degradation constant - the larger the degradation rate (pyocins degrade faster without meeting a bacterium), the longer it takes to kill all the bacteria. The right graph shows the relationship between the total bacteria elimination time and the pyocins killing efficiency. We can see that the connection is reverse in comparison to the one shown on the left - the higher the efficiency of killing, the faster the entire bacterial population dies.
We emphasize that the killing time always turned out greater than 20 minutes as expected, since this is the delay time (d), the time it takes for the pyocins to kill a bacterium from the moment of binding.
2D Model and simulation
In order to simulate the pyocin-bacteria interaction on a plate, the FTCS method was used to model the diffusion and the interaction on a 2D surface:
2D bacteria equation
$$ \frac{B[x,y,n+n_0]-B[x,y,n]}{n_0} = D_B (\frac{B[x+h_0,y,n]-2B[x,y,n]+B[x-h_0,y,n]}{h_0^2}+ $$ $$\frac{B[x,y+m_0,n]-2B[x,y,n]+B[x,y-m_0,n]}{m_0^2}) +rB[x,n] (1-\frac{B[x,n]}{K})$$ $$- \alpha u[n-d] B[x,y,n-d] P[x,y,n-d]$$2D pyocins equation
$$ \frac{P[x,y,n+n_0]-P[x,y,n]}{n_0} = D_P (\frac{P[x+h_0,y,n]-2P[x,y,n]+P[x-h_0,y,n]}{h_0^2}+ $$ $$\frac{P[x,y+m_0,n]-2P[x,y,n]+P[x,y-m_0,n]}{m_0^2})-m P[x,y,n] $$ $$- \alpha u[n-d] B[x,y,n-d] P[x,y,n-d]$$Where $h_0$ - step in x axis, $m_0$ - step in y axis, $n_0$ - step in time.
The populations' interaction and growth were simulated over time. The 2D model was validated with an article's[5] results, which describes a similar nature of interaction with phages.
The pyocin's amount over time is plotted on the left side of the video and the bacterial amount on the right side. The Z axis and the heat map colors represent the populations' amount in the current time step, and X, Y axes represent the space. In the beginning of the simulation, the time is still smaller than the delay time (20 minutes). Therefore, the populations are not affected by the interaction - the bacteria continue to grow and the pyocins diffuse in space. After crossing the delay time, the interaction is starting to affect both populations and is translated to a decline in the populations' amount.
2D simulation video
Back
Pyocins distribution optimization
In addition to the model, we developed an algorithm that enables optimizing the distribution of pyocins in the environment for efficient usage (i.e. minimizing the time and the number of pyocins needed for bacterial elimination). The algorithm relies on our 2D simulation, which is used to calculate the time it takes to eliminate all bacteria in the environment.
The number of options for the position of the pyocin dots can sum up to hundreds of millions of options (even for a small grid and a small number of drops). Therefore, to reduce the running time of the algorithm we used a greedy heuristic algorithm, which practically converges to a local optimum (in terms of bacteria elimination time) efficiently. At each step, the algorithm performs a change on the position of one of the pyocin dots, such that the improvement of the objective function will be maximal. In addition, for a smooth convergence to a global or local optimum, the search steps are decreased gradually during the iterations.
We developed an easy to use software for easy usage of the algorithm mentioned above. When using this software, the user can define the resolution of the problem. In both X-axis and Y-axis, the step size is 10 um. The user can also define the diffusion coefficient in the different regions in the environment for both bacteria and pyocins; The user can choose the number of pyocins and number of pyocin dots, and define the problem's grid.
In addition, the user can limit the algorithm's running time, in order to get the best solution found until the specified time. If after a few iterations the initial pyocins amount isn't enough to eliminate all bacteria in the user's surface, the user is requested to enter a new number of initial pyocins amount. Below are two examples of the algorithm's output.
1. Equal diffusion in space for both pyocins and bacteria; 3 pyocin drops; 20*20 grid:
Drop #1: x location-6, y location-11
Drop #2: x location-11, y location-3
Drop #3: x location-18, y location-17
2. Equal diffusion in space for bacteria, and higher pyocins' diffusion coefficient in the lower half of the grid; 4 pyocin drops; 50*50 grid:
Drop #1: x location-4, y location-21
Drop #2: x location-18, y location-25
Drop #3: x location-18, y location-11
Drop #4: x location-31, y location-31
As expected, the pyocin drops' density is higher in the upper half of the grid, because of the difference in the pyocins' diffusion coefficient.
The software guide is attached under the appendix section.
Back
Bioinformatics Models
Part Improvement
The improvement of an existing part is an important aspect of our iGEM project. In order to do so, we decided to take a part that has been used by many teams and increase its expression. For this purpose, we used bioinformatics tools that allowed us to obtain synonymous mutations that improve the expression without altering the amino acid sequence. To obtain this, we used a genetic algorithm, which models the folding energy across several generations. In each one of its generations, the genetic algorithm inserted mutations and looked for the best variants. You can find out more about this algorithm and how it works in our bioinformatics work page. Later on, we created (in real life) the best variants and tested them in the lab in comparison to the original part. You can read more about the wet lab experiment of the part improvement in part improvement page.
BackDiscovery of a new tail
In order to target new bacteria, we needed to find new tails which will bind to the pyocin’s baseplate but will target new type of bacteria. To do so we used our bioinformatics software ‘Tail-or Swift’. The software uses our own alignment algorithm and models possible new tails (you can read all about it in our software page).
We ran the software and looked for possible candidates. We found out that the protein ‘P2 gpH-like protein’ of the virus Salmonella virus Fels2 has a high similarity in its N-terminal part to the original Pyocin’s tail (about 43%). Therefore, according to our model, it has a high probability of binding to the pyocin’s baseplate. In addition, the Salmonella virus Fels2 is a known bacteriophage, which targets different strains of Salmonella. Due to the high similarity and the fact the Salmonella is a known pathogen, harming many people annually, we decided to use this tail.
To test this new tail, we had to acquire a Salmonella strain that is sensitive to this bacteriophage. We believed that if the Salmonella is sensitive to this bacteriophage, and we will use its tail in our pyocin system, our engineered pyocin will be able to target these bacteria. We contacted Prof. Ohad Gal-Mor from Shiba hospital who provided us Salmonella typhimurium LT2. Salmonella typhimurium LT2 is a Salmonella strain which has Salmonella virus Fels2 as a prophage in its genome, and therefore is supposed to be sensitive to this bacteriophage.[9][10]
At the same time, we created our new tail and inserted it to our pyocin system.
This tail and target duo was one the tails and targets that we tested our system on.
Appendix
Refrences
[1]: J. Crank. The mathematics of diffusion (1975), 3-4
[2]: N. Keyfitz. International Encyclopedia of the Social & Behavioral Sciences (2001).
[3]: K.S. Al Noufaey, T. R. Marchant and M. P.Edwards. The diffusive Lotka-Volterra predator-prey system with delay (2015), 1-2
[4]: Stephen a. Gourley and YangKuang. A delay reaction-diffusion model of the spread of bacteriophage infection (2005).
[5]: D. Ping, T. Wang et al. Hitchhiking, collapse and contingency in phage infections of migrating bacterial populations (2018).
[6]: https://bionumbers.hms.harvard.edu/search.aspx
[7]: M. Thongmoon and R.McKibbin. A Comparison of Some Numerical Methods for the Advection-Diffusion Equation (2006),51
[8]: Michel-Briand Y, Baysse C. The pyocins of Pseudomonas aeruginosa (2002),3
[9]: Mohammed, M. & Cormican, M. Whole genome sequencing provides possible explanations for the difference in phage susceptibility among two Salmonella Typhimurium phage types (DT8 and DT30) associated with a single foodborne outbreak. BMC Research Notes 8, 728 (2015)
[10]: Frye, J. G., Porwollik, S., Blackmer, F., Cheng, P. & McClelland, M. Host Gene Expression Changes and DNA Amplification during Temperate Phage Induction. J. Bacteriol. 187, 1485 (2005)
Diffusion-Interaction Model MATLAB codes
Pyocins distribution optimization software guide
Back