Introduction

To get an idea of how our gene circuit would work on an ideal situation, we explored the structure and dynamics of our system by creating a mathematical model of the reaction kinetics and a real-time simulation. Wit the mathematical model, we started with no information ready whatsoever and derived differential equations to demonstrate how the use of blue light and the changes in phosphorylated $YF1$ and $FixJ$ concentrations would control the production of our three target proteins. We labeled them simply $A$, $B$, and $C$, because the system is intended to be used with any user defined three genes coding different proteins.

Using our equations we constructed a simulation showing the effects of red and blue light on our system in real time. The user can control the input of both lights to see how they affect the production of proteins $A$, $B$ and $C$. We experimented with different values for all constants and with trial-and-error iteration we arrived to a visualized simulation that can be used to demonstrate the intended funcion of our system. All in all, this is an idealization. Based on present and future measurement data, the parameters can be adjusted to better match the real world.

Mathematical model

Overview

Here, we will discuss the dynamics and interactions of different proteins and promoters introduced in the research section. We have developed a simplified mathematical model describing our system.

Simplifications

We assumed that the species identified from our gene circuit are the only ones that affect the overall concentrations inside the bacteria. The bonding of $CI$ to $O_R$ sites is assumed to be insignificant compared to the overall concentration. The model is also strictly deterministic and doesn’t take into account any noise. The phosphorylation, decay, bonding and production of proteins are assumed to be linear functions of concentration. We further assumed that the phsophorylation of $FixJ$ by phosphorylated $YF1$ would not involve phosphorous transition between the reacting molecules

The first model that was constructed before our lab work had even begun involves many harsh simplifications. Our aim was to get a general picture of how the system could work in ideal conditions and how stable it was.

Equations for dynamics

Based on the assumptions made before, we arrived at following differential equations to describe the idealized dynamics of our system:

\begin{eqnarray*} & & \frac{d[YF1]}{dt} = P_1Rbs_1 + DP_1[YF1]_{phos} - (Deg_{YF1}+I_B)[YF1] \\ \quad \\ & & \frac{d[YF1]_{phos}}{dt} = I_B[YF1] - (Deg_{\text{YF1}} + DP_1)[YF1]_{phos} \\ \ & & \frac{d[FixJ]}{dt} = P_1Rbs_1 + DP_2[FixJ]_{phos} - (C_{phos}[YF1]_{phos} + Deg_{FixJ})[FixJ] \\ \\ & & \frac{d[FixJ]_{phos}}{dt} = C_{phos}[YF1]_{phos}[FixJ] - (DP_2[FixJ]_{phos} + Deg_{FixJ}[FixJ]) \\ \\ & & \frac{d[CI]}{dt} = P_2Rbs_1 - Deg_{CI}[CI] \\ \\ & & \frac{d[TetR]}{dt} = (P_A + P_B)Rbs_2 - Deg_{TetR}[TetR] \end{eqnarray*}

Legend

• $[YF1]$ = concentration of $YF1$ protein
• $[YF1]_{phos}$ = concentration of phosphorylated $YF1$ protein
• $[FixJ]$ = concentration of $FixJ$ protein
• $[FixJ]_{phos}$ = concentration of phosphorylated $FixJ$ protein
• $[CI]$ = concentration of $CI$ protein
• $[TetR]$ = concentration of $TetR$ protein
• $P_1$ = relative strength of the first promoter in gene circuit
• $P_2$ = relative strength of the $FixK_2$ promoter
• $P_A$ = relative strength of the $P_R$ promoter, codes gene A
• $P_B$ = relative strength of the $P_{RM}$ promoter, codes gene B
• $Rbs$ = relative strengths of ribosome binding sites
• $Deg$ = degradation coefficient
• $C_{phos}$ = phosphorylation coefficient
• $DP$ = de-phosphorylation coefficient

These equations describe the essential proteins our system ($YF1$, $FixJ$, Phosphorylated $YF1$, Phosphorylated $FixJ$, $CI$, $TetR$). Proteins are produced with rates that depend on the strength of respective promoter and ribosome binding site, and also when phosphorylated protein (denoted with phos) is dephosphorylated back to its original form. The concentration of all proteins reduces by degradation, the rate of which depends on the concentration of protein in question.

Dynamics' coefficients

$P_1$, $P_2$, $P_A$ and $P_B$ denote the relative strengths of promoters. These change and $Rbs_1$&$Rbs_2$ the relative strengths of ribosome binding sites, which both affect the protein synthesis linearly. Each protein has its own degradation coefficient (denoted $Deg$). $I_B$ is the combined effect of blue light that affects the phosphorylation of $YF1$. The phosphorylation of $FixJ$ is assumed to depend on phosphorylation coefficient $C_{phos}$ and the concentration of phosphorylated $YF1$. The dephosphorylation here depends on a respective dephosphorylation coefficient $DP$(1&2 for $YF1$ and $FixJ$).

Equations for promoter activities

\begin{eqnarray*} & & P_2 = C_{P_2}N[CI] \\ \\ & & P_A = \begin{cases} C_{P_A}N[CI] \quad \text{if} \quad N[CI] \leq 1 \\ 0 \quad \text{if} \quad N[CI] > 1 \end{cases} \\ \\ & & P_B = \begin{cases} C_{P_B}N[CI] \quad \text{if} \quad N[CI] < 1 \\ C_{P_B}(1-(N[CI] - 1)) \quad \text{if} \quad 1 \leq N[CI] < 2 \\ 0 \quad \text{if} \quad N[CI] \geq 2 \end{cases} \\ \\ & & P_C = \begin{cases} C_{P_A}(1-N_2[TetR]) \quad \text{if} \quad N_2[TetR] \leq 1 \\ 0 \quad \text{if} \quad N_2[TetR] > 1 \end{cases} \end{eqnarray*}

Promoters' coefficients

Here, the $C_{P_n}$'s denote the respective promoter's maximum activity. The $N$ in front of $CI$ and $TetR$ concentrations is a normalization coefficient, which is needed to map the fitting part of concentrations' values to the interval $(0,1)$. This way, when multiplied by the promoters' maximum activity values, we get values in the interval $(0, C_{P_n})$ The fuctions definitions must also change so that they never take negative values, which would make no sense when talking about promoter activities. We have simplified the model so that the promoter activities only depend on the concentration of $[CI]$ protein. We also assumed that the binding of $CI$ had insignificant effect on it's concentration, and that the amount of $CI$ binded to the $O_R$ sites was propotional to $CI$ concentration.

Simulation

Overview

Based on our mathematical model, we created an interactive simulation and a graphical user interface for it. This visualization, although idealized, is ideal for demonstrating the intended functioning of our gene circuit and gene switch system. We included two swithces, one for red and one for blue light. With these, the user can see the effect of light intensity to a bacterial coulture in real time. Here, proteins A, B and C are represented by GFP, RFP and BFP, so the bacteria change color in different circumstances.

Lights

In our system, the communication between user and the bacteria happens via illumination of blue light to the coulture. Blue light phosphorylates the $YF1$-protein, which is the key to controlling the overall protein production inside the bacterium. In the simulation, this is represented by change in the $I_B$ parameter from the mathematical model. This takes values between 0 and 1, and the rest of system behaves as described earlier.

Our original design had also an intensity switch, operated by red light. Due to time constraits, this wasn't yet implemented in our gene circuit. In the simulation, we added a second user controlled parameter in front of every promoter. This takes values between 0 and 1, representing the no production at all -state and the production at maximum promoter activity. With this, the user has control of all protein concentrations. The assumed mechanism is idealized and has a linear effect on the activity.

Runge-Kutta method

The dynamics of our system were approximated using 4th order Runge-Kutta (RK4) method for the differential equations in our mathematical model. The point of this method is to approximate the function in question by it's derivatives without having to solve the function itself. The starting values of each concentration are assumed to be zero, so $y(0) = 0$. The simulation computes the next datapoint adding the derivative times a timestep $h$ to previous concentrations. The method uses a mean value of different derivatives (the different k's below) during this timestep h to get a more accurate approximation.

\begin{eqnarray*} & & y' = f(t,y(t)), \quad y(t_0) = y_0 \\ \\ \\ & & y_{n+1} = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) \\ & & t_{n+1} = t_n + h \\ \\ & & k_1 = f(t_n,y_n) \\ & & k_2 = f(t_n +\frac{h}{2}, y_n + \frac{h}{2}k_1) \\ & & k_3 = f(t_n + \frac{h}{2}, y_n + \frac{h}{2}k_2) \\ & & k_4 = f(t_n + \frac{h}{2}, y_n +hk_3) \end{eqnarray*}

Software implementation

A computational model was created based on our mathematical model and the RK4 approximation. We made a real-time plotting function to illustrate the dynamics in a simple and clead graphic UI. We added two light switches so that the user can have an impact on our simulation in real time. This all was then further developed into a presentable, user-friendly form that is accessible from our website. The simulation itself was created using Python and translated into Javascript for web implementation.

The UI of our simulation

To demonstrate our work to the general public in an event called Summer of Startups Demo Day, we made a simulation that shows our system in action. It shows an animated bacteria plate with adjustable light intensity sliders to remotely control the bacteria. The proteins the bacteria produce are colors, so you can see how the changes in light intensity correlate to the color of the colonies on the plate. The simulation also has a nice grap that shows the protein levels in one second intervals so you can see more clearly what's going on in the cell.

All of the code (including a more in-detail Python graph simulation) is available at the project's GitHub page.

Discussion

In biological systems, a huge number of variables affect the outcome of any process. Most of them are unknown, as in our case, so models generally include some noise to approximate these unknown factors. However, our model is simplified: it takes no noice in consideration. Therefore all interactions produce smooth, good-looking curves in the graphic output. From another point of view, the clear graphipcs produced by the simulation are easily interpretable, so even someone not familiar with science can clearly see what's going on "behind the scenes" of our system.

So far we have also used arbitrary parameters, simplified reaction pathways and reaction equations. The parameters were acquired by estimation and empirically testing. Full experimental data wasn't available when the simulation was created, so derivation of differential equations by using the law of mass-action wasn't possible. All reaction mechanisms are our own estimations of what's going on inside the bacteria.

In it's current state the simulation gives a good idea on how the system should work. Making it realistic and really accurate requires a lot of measurement of appropriate parameters, research on apprpriate scientific publications, and loads of tuning. Still, this version is ideal for demonstration of our idea. bringing the visual UI a lot of business value. Assuming we had a working three-channel switch, marketing and advertising our technologies and solutions would be considerably more fruitful thanks to the credibility brought by the client trying the technology himself in a real-time simulation.