Introduction

To get an idea of how our gene circuit would work in an ideal situation, we explored the structure and dynamics of our system by creating a mathematical model of the reaction kinetics and a simulation that can be controlled in real time. We started working with the mathematical model without any detailed information about the system. We derived the differential equations to demonstrate how blue light and the changes in phosphorylated $YF1$ and $FixJ$ concentrations would affect 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 three genes encoding the proteins of choice.

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 via trial-and-error iteration we arrived to a visualized simulation that can be used to demonstrate the intended function of our system. This is an idealization. Based on present and future measurement data, the parameters can be adjusted to better the dynamics of our system.

Mathematical Model

Assumptions

We assumed that the species identified in our gene circuit are the only ones that affect the overall concentrations inside the bacteria. We further assumed that binding of $CI$ to $O_R$ operator site does not impact the overall $CI$ concentration and that the amount of $CI$ bound to the $O_R$ sites was proportional to $CI$ concentration. The model is also strictly deterministic and doesn’t take any noise into account. The phosphorylation, decay, binding and production of proteins are assumed to be linear functions of concentration. We further assumed that the phosphorylation of $FixJ$ by phosphorylated $YF1$ would not involve phosphate transfer between the reacting molecules.

The first model was constructed before our lab work had even begun and it contains 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_{YF1} + 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_{FixJ} + 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_{CI} - Deg_{CI}[CI] \\ \\ & & \frac{d[TetR]}{dt} = P_ARbs_{TetR_A} + P_BRbs_{TetR_B} - Deg_{TetR}[TetR] \end{eqnarray*}

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 and its depends on the concentration of protein in question.

Rate Coefficients

$P_1$, $P_2$, $P_A$ and $P_B$ denote the relative strengths of the promoters. $Rbs$s are the relative strengths of ribosome binding sites, which both affect the mRNA translation rate 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 the respective dephosphorylation coefficient $DP$(1&2 for $YF1$ and $FixJ$). Later on, we found out that non-phosphorylated YF1 acts as a phosphatase on FixJ. However, these effects are not taken into account in our model.

Equations for Promoter Activities

\begin{eqnarray*} & & P_2 = C_{P_2}N_1[CI] \\ \\ & & P_A = \begin{cases} C_{P_A}N_1[CI] \quad \text{if} \quad N_1[CI] \leq 1 \\ 0 \quad \text{if} \quad N_1[CI] > 1 \end{cases} \\ \\ & & P_B = \begin{cases} C_{P_B}N[CI] \quad \text{if} \quad N_1[CI] < 1 \\ C_{P_B}(1-(N_1[CI] - 1)) \quad \text{if} \quad 1 \leq N_1[CI] < 2 \\ 0 \quad \text{if} \quad N_1[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*}

Promoter Coefficients

Here the $C_{P_n}$s denote the respective promoter's maximum activity. The $N_1$ and $N_2$ in front of $CI$ and $TetR$ concentrations are normalization coefficients, which are needed to map the values of $[CI]$ to the interval $(0,3)$ and values of $[TetR]$ 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 functions definitions must also change so that they never take negative values, which would make no sense when it refers to promoter activity. We have simplified the model so that the promoters’ activity only depend on $[CI]$ and $[FixJ]$.

Simulation

Overview

Based on our mathematical model, we created an interactive simulation and a graphical user interface for it. This visualization, although idealized, is suitable for demonstrating the intended functioning of our gene circuit and the gene switch system. We included two sliders, one for red and one for blue light. With these, the user can see the effect of the light intensity to the simulated bacterial culture in real time. Proteins $A$, $B$ and $C$ are represented by GFP, RFP and BFP (green, red and blue fluorescent protein) and therefore the bacteria change color when lights’ intensities are changed.

Lights

In our system, the communication between user and the bacteria happens via illuminating the culture with blue and red light. Blue light phosphorylates the $YF1$-protein, which is the key to controlling the production of $A$, $B$ and $C$ proteins inside the bacteria. 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 the system behaves as described previously.

Our original design also had a transcription intensity switch, controlled by red light. Due to time constrains, 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 zero production state and the production at maximum promoter activity. With this, the user has control of all desired 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 and computed using 4th order Runge-Kutta method (RK4) 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 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 visualisation script to illustrate the dynamics in a simple and clear 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.

Here is a screenshot of the simulation. You can adjust the amount of the red and blue light and see how it affects the bacteria. You can also see how active each gene (A, B, C) is.

To demonstrate our work for the general public in an event, Summer of Startups Demo Day, we used the simulation to show our system in action. It can be accessed with a web browser and shows an animated bacterial plate with adjustable light intensity sliders to remotely control the bacteria. The proteins the bacteria produce in this simulation 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 graph that shows the protein levels in real time so you can see more clearly what's going on in the cell.

All the code (including the Python simulation with more detailed graphs) is available at the project's GitHub page.

Discussion

Our model doesn’t take any noise into consideration. Therefore all interactions produce smooth, good-looking curves. On the other hand, the clear graphics generated by the simulation are easily interpretable, so even someone not familiar with science can clearly see what's going on in our system.

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

Some unexpected observations were made after running the simulation several times. We noticed that when activating all the promoters while the $CI$ concentration was zero, both proteins $A$ and $C$ were produced simultaneously instead of just the anticipated $A$. When $TetR$ is further produced by the activation of $P_R$ promoter, $C$ production is repressed and the concentration drops back to zero. Secondly, when blue light intensity is set to a level that corresponds to the maximum concentration of either $A$ or $B$, the promoter activity adjustable using the red light and it should only affect the said concentration. Again, when lowering the activity enough, we noticed that a spike in production of protein $C$ appeared again. This seemed to be caused by the lowered concentration of $TetR$ that allowed a leak in $P_C$ promoter. We had no idea that a $C$ spike would appear based on the theoretical model of our gene circuit, so this phenomenon was discovered early thanks to our simulation.

We also noticed that going from directly producing the protein $A$ to protein $C$, or reversely, from $C$ to $A$ is virtually impossible without producing some protein $B$ along the way. We thus concluded that our Gene Switch is not entirely orthogonal between the three channels. The possible interactions with other products with protein $B$ are needed to be taken under consideration when designing applications that only use genes $A$ and $C$.

Upon later research, we found out the actual mechanism with which $FixJ$ was phosphorylated. In contrast to our model, the phosphate is actually transmitted from $YF1$ to $FixJ$. In their paper, Möglich et al. (2009, reference in research section) showed that in a two-step reaction, $FixL$ first undergoes autophosphorylation and then transfers the phosphate to its cognate, noncovalently bound, response regulator $FixJ$. The $YF1$ protein is a derivative of $FixL$ with different sensory domain, so it behaves the same way in this reaction. This wasn't however implemented in our model.

In it's current state the simulation gives a good idea on how the system should work. Making it realistic and accurate requires measuring the appropriate parameters, research on appropriate scientific publications, and tuning. Still, this version is ideal for demonstration of our idea, bringing the visual UI a significant marketing value.