In order to find the best way to implement the toxin-antitoxin system, we resort to modelling. We use the gillespie algorithm to model the interactions of the toxin antitoxin system.
We find that when we implement enhanced relE production as a tool to make the bacteria dormant, an additional implementation of relB to ensure don’t stay dormant when in light again.
The model found that the system is sensitive to the relE:relB ratio as well as the total production, and that an implementation with production rates in the vicinity of 50 and 35 molecules pr. min for relB and relE respectively yields close to the wished for effect: THe bacteria goes dormant in an hour and wakes up quickly.
Gillespie Algorithm
Approach
The Gillespie algorithm is a way to get calculate the evolution of stochastic functions; in this case cell concentrations. To use the algorithm, two things are needed:
reaction rates of the system (a1, a2, a3) at a given configuration.
a random number generator.
For each time step two things are calculated, using the random number generator: The time before next reaction and which reaction occurs.
The time before next step is given by ( Figure (well... equation): reaction time insert here (tex doc))
Where S is the sum of the reaction rates and r1 is a random number between 0 and 1. This gives the time as if the system was one reaction with reaction rate S, using the random number to give a statistic distribution.
The reaction is chosen proportionally to each individual reaction rate, using another random number. This way, reaction with high rates compared to other reactions will happen the most.
The reaction is carried out, the time is moved the calculated amount, and new reaction rates can now be calculated, to repeat the whole process. This continues until the time reaches the wanted limit or a specific number of reactions have occurred. It is necessary to have a limit on the number of reactions as it is possible for the time steps to grow smaller and smaller, in which case the calculation time quickly become either immense or impossible.
Notice that all values are integers as the gillespie algorithm works with discrete numbers of molecules. The values were chosen based on a stable equilibrium found by letting the model run a simulation of the inherent system over 450 minutes, with different starting values.
The implemented total production rates shown in the model might seem too high as they range from 1-350 molecules pr. min, while the rates in the inherit system is effectively around 80-100 for relB and 2-5 for relE. The possibility of placing the system on high-copy plasmids, however, makes the high total production values reasonable, as the individual relE promoter only need a production of 0.01-2 (assuming a few hundred copies).
Considering the inherent toxin-antoxin system activating under starvation, we see that the the magnitude of relE copies around 50-80 molecules pr. cell. This makes it reasonable to believe that the cell enters dormancy when a few tens of free relE copies.
Running the model
The model is using the gillespie algorithm to give a stochastic view of the system and is run on matlab. The code uses an implementation by matlab user Nezar (https://se.mathworks.com/matlabcentral/fileexchange/34707-gillespie-stochastic-simulation-algorithm).
For the dormancy runs we input deterministic initial values, and then let the system run 30 min without activating the inserted toxin promoter. This results in a stochastic distribution of initial values mimicking variations between cells. Analysis show that 30 minutes is enough for the model to find a stable distribution, which is realistic considering the growth cycle of an e.coli.
For wakeup runs, we use the data generated at the end of a sleep run as initial value and remove the production of relE. We consider a cell to be woken up when the concentration of free relE decreases below 15 copies. This is perhaps a low value, but tests show marginal difference between 15 and 45 copies, where the low bound is chosen to decrease uncertainty of the cells state.
All runs simulate 1000 cells, which should be sufficient to get stable averages. The model assumes well-mixed conditions in each cell, but considers each cell independently.
The model has no cut off for maximum values of relE, as we don’t know the exact relation between relE concentration and hibernation state, yet as a functional cutoff is found through wake-up times, it is not necessary.
