Circuit Modeling

Overview

Biological modeling is an indispensable element in scientific research, which deals with massive size of data, analysis and prediction that helps to optimize parts, especially in synthetic biology, optimizing circuit parts is always one of the goals of DRY LAB. Mathematical modeling of biological elements gives us insight on how a system might work before or during the experimental works, which protocols can be refined according to the model.

Engineering achievements in short

  • Predict the time for result
  • Shorten the waiting time for showing results
  • Optimize and create more controllable circuits

    • Potassium Ionophore Modeling

    Introduction

    The very first step of our circuit is the triggering of efflux of potassium efflux by the cereulide, which is our interested toxin and potassium ionophore. Because of this, we are interested in modelling the potassium efflux caused by cereulide:

    1. The rate of efflux versus time with different initial cereulide concentration
    2. The concentration of intracellular & extracellular potassium ion versus time with different initial cereulide concentration

    Unfortunately, we couldn't find literature finding the relationship between cereulide and potassium ion efflux. However, based on the assumption that cereulide and valinomycin(which is also a potassium ionophore) are similar, we found a paper[1] that describes the kinetics of valinomycin-mediated K+ ion transport. By assumption, we modeled the kinetics according to the figure below:

    Figure 1 Kinetics of valinomycin-mediated K+ ion transport

    We wrote the ode and input the constants given from the paper:

    ODE

    • Results

    Figure 2 Rate of potassium efflux/influx that is due to the K+ ionophore([Val]=0.4e-3M)
    Figure 3 [K+] versus time

    The higher the initial concentration of valinomycin, the higher the rate of potassium efflux is, the faster it achieves equilibrium.(Fig 5)

    Figure 4 Overall rate change of [K+]
    Figure 5 Changes of the rate of potassium efflux under different initial concentration of valinomycin.

    Phosphorelay Modeling

    Introduction

    Assumptions of our model

  • The activation of KinC is directly proportional to the rate of potassium efflux.
  • Parameters are found based on 25 degree celsius.
  • B subtilis is not under starvation stage. Therefore, KinA won't be activated.

    • When the K+ efflux occurs, it will trigger the phosphorylation of KinC, a member of histidine kinases family. KinC~P will continue triggering a series of downstream phosphotransfer reactions and produce Spo0A~P, the transcription factor which we'll be using to trigger our circuit, in the end. We modify the model The model is referenced from the paper[3] for the phosphorelay modeling in this section.

    Figure 6 Para_Post-translational reactions and parameters of the phosphorelay model

    Results

    Figure 7 K+ efflux and KinC phosphorelay

    The concentration of Spo0A and Spo0A~P in unit (uM) versus time when input of cereulide is 0.4mM, which is our toxicity threshold.

    Circuit Modeling

    Introduction

    Assumptions of our model

  • Assuming first-order kinetics as the governing principle.
  • The model has been reduced using the Quasi-Steady-State Approximation (QSSA) and Rate Equations Approximation (REA).
  • The system operates at 25 degree celsius since parameters are found mostly at this temperature.

    • Cycles

    Figure 8 Cycle0 circuit design

    ODE

    • Modeling Results
    Figure 9 Cycle0 - RFP concentration versus time
    Model Fitting with experimental data

    We consulted Professor Henry Lam before we started the model fitting as the actual data is too far from our expected (check drylab external help). From the meeting, we are inspired to use genetic algorithms (GA) to search for the possible parameter space for our model first and do model fitting with those parameters. Then in the future we would fit the model with local optimization algorithms. We proposed the idea of using 2-step fit for model fitting in the future.

    Here we show the circuit model at cereulide=0.3uM fitted with experimental data by using GA.

    Figure 10 KinC Phosphorylation System versus time after model fitting
    Figure 11 [Spo0A] & [Spo0A~P] versus time after model fitting
    Figure 12 [RFP] versus time after model fitting

    Our fitted model output (red) did not fit well to the experimental data (blue). Our predicted time for signal outputting is 2 hours earlier than the actual one. Not only did it show the signaling cascade of KinC to Spo0A~P takes time, it also pointed out further tuning of the assumptions and parameters of our model is needed, so that we can get to our modeling goals, which are shortening our result output time and further improving our circuit with the fixed and tuned model.

    The result of model fitting is not good, with Mean squared error (MSE)=7.33. It might due to:

    1. The assumptions of Potassium efflux and KinC relationships.
    2. The assumptions of first-order kinetics.
    3. The assumptions of neglecting the other mechanisms that is related to Spo0A.
    4. Not enough experimental data.
    5. Fitting too many parameters at the same time.
    6. Our GA reaches a local minimum.

    Future steps for model fitting to fine tune the assumptions:

    1. Analyse circuit experimental data with fold-change, try to find out the relationship of KinC and Potassium efflux by analysing the input toxin concentration to the [RFP] at equilibrium.
    2. Conducting experiments with a larger range of cereulide concentration and performing model fitting again.
    3. Do the fitting with GA with a larger range of boundaries.

    Cycle0 only has RFP as result, which doesn't have a negative control. In this state, our test kit is only able to warn the user when the food contains excess cereulide, which is not safe to eat. However, the user might not be able to know if the test kit is functioning well and whether the rice is really safe to eat. Due to that reason, we have Cycle1 with both RFP and GFP as a result. GFP will act as negative control signal showing that the our test kit is working normally and the food doesn't contain cereulide.

    Figure 13 Cycle1 circuit design

    We have added a conditional degradation system for cycle 1. Promoter hyper spank is a constitutively expressing promoter, which would keep expressing GFP-ssrAtag at the beginning until repressed by LacI which is expressed from Promoter sinI. At the same time when the promoter hyperspank is activated, it would express sspB and RFP also. The sspB is an adapter that targets the ssrA-tagged GFP and brings it to the protease clpXP for degradation. Eventually, the GFP would be degraded in a higher rate and RFP would have a higher fluorescence than GFP, minimizing the color mixing issue and shortening the result time.

    ODE & Model reduction

    Repressible mechanism

    Modeling Results
    Figure 14 Cycle1 modeling graph, [Cereulide] = 0.3uM ( < toxicity threshold)
    Figure 15 Cycle1 modeling graph, [Cereulide] = 0.5uM ( > toxicity threshold)
    Conditional degradation tag performance

    GFP fluorescence attained equilibrium at around 1 hour, which is much faster than RFP. As GFP is degraded by a strong conditional degradation system, which could degrade most of the GFP in a system around 20 to 30 mins according to the paper[4].

    Interpretation on our circuit

    With higher input (higher concentration of toxin), the rate of RFP production and GFP degradation is increased.

    While the goal of our test kit is to show a clear cut of red and green signal when the rice is unsafe or safe, our current model shows our circuit failed to achieve that by cycle 1 circuit. Both of them are showing the GFP would eventually be lower than RFP (with a time difference of having RFP>GFP = 3 mins), and show low responsiveness to the input.

    One thing to notice as well is the degradation of GFP by sspB, which is the adaptor that performs conditional degradation. SspB is not present naturally in B. subtilis, but expressed by our circuit through the promoter sinI. While RFP is being produced, sspB and LacI (repressor of Promoter hyperspank) are also being produced, which both sspB and repressor would cause the reduction of GFP concentration, one by degrading and another one by repressing the expressing of GFP respectively. This means that GFP concentration is strongly correlated to RFP expression (including the basal expression), once RFP is produced, GFP would be degraded with higher rate, but RFP is not affected by mechanisms related to GFP. One of the reasons our circuit failed might be due to the leaky promoter or the strong expression of the promoter.

    We would want the RFP signal to be smaller when the circuit is under the toxicity threshold. Possible way to reduce RFP signal when the input is below the toxicity threshold is by reducing the RFP signal when GFP is produced. As with higher concentration of cereulide, the faster our circuit produces RFP. This increased RFP production is associated with the faster production of sspB and LacI, which in turn affects the degradation and reduction of GFP. By introducing more parts to our circuit, It is possible that we can manipulate the circuit to control the threshold of expressing GFP and RFP signals and increase the responsiveness of our circuit (check circuit modeling cycle 2). Further analysis on the model is needed in order to design additional parts for circuit manipulation. (Check Cycle 1 SA)

    Comparing our model to circuit design data

    From our circuit design data, the data shows low responsiveness of the circuit to the input as well. (check circut experiments cycle 1)

    Figure 16 Cycle 1 circuit experiment result 1
    Figure 17 Cycle 1 circuit experiment result 1
    Engineering Cycle Proposal

    In order to achieve our goal according to the input concentration of cereulide, we would first:

    1. Reduce the RFP signal
      1. Why? How does it affect us?
        1. As one of our goals is to have a clear cut for our test kit according to the threshold of toxicity, the current RFP signal from the leaky promoter is not our desired output. We are looking for the signal that depends on the cereulide concentration but not due to the basal expression.
      2. Our plan:
        1. Perform Sensitivity Analysis (SA) on the circuit and check which parameters (mechanisms) contribute most to the RFP output.
        2. Change the parts in the model according to the SA result and analyse the impacts of modifying most sensitive parts (Check Cycle 2).
    2. Reduce the basal expression of PsinI
      1. Why? How does it affect us?
        1. At the beginning of the circuit, GFP is constitutively produced. High basal expression of PsinI would express sspB and LacI at the very beginning, which the GFP would be degraded due to the basal expression of PsinI but not the activation effect of our input (Cereulide concentration).
      2. Our plan: Changing the promoter which is less leaky
        1. Dry Lab exploring different promoter with different basal expression.
        2. Replace the PsinI with possible candidates to our model and choose the most suitable on for our circuit.
    Sensitivity Analysis
    Figure 18 Sensitivity analysis of Cycle1 RFP
    Figure 19 Sensitivity analysis of Cycle1 GFP

    From our sobol analysis on the model, we discovered that several parameters are significantly more sensitive to the output of RFP (alpha0, alpha, dr1, kt1, d1). In this case, we could utilise the result from SA of Cycle1 RFP and reduce the RFP signal. As our circuit is using the ssrA tag for conditional degradation, dry lab proposed a new idea to circuit design which is to add one more conditional degradation system to RFP, utilising the high sensitivity index of the degradation rate of RFP, we should be able to reduce the impact of the strong basal expression of RFP. Before the actual decision and experiment for our proposed circuit, we started with modelling our idea to prevent wasting manpower and time on constructing a new circuit, and also find out possible errors of the circuit.

    According to another paper studying conditional degradation[5] they introduced CCsspB, which only targets the substrate linked with CCssrA. We first tried this in our model for RFP degradation. In cycle2 modeling, we assumed that CCsspB / CCssrA kinetics and sspB / ssrA kinetics and parameters are similar.

    Figure 20 Cycle2 circuit design
    Additional ODE

    Repressible mechanism

    Modeling Results
    Figure 21 Cycle2 modeling graph, [Cereulide] = 0.3uM ( < toxicity threshold)
    Figure 22 Cycle2 modeling graph, [Cereulide] = 0.5uM ( > toxicity threshold)

    After we integrate CCsspB into our circuit model, we can see that for the concentration (0.3uM) which is below the threshold, the GFP is always larger than RFP signal. While for (0.5mM) which is above the toxicity threshold, the RFP signal would exceed the GFP signal at around certain time.

    Cycle 2 successfully shows that our circuit might be able to produced a clear cut signal when the rice is safe or not with the dual conditional degradation system.

    Future plan:
    1. How can we further tune the circuit so that it works like a sigmoid function (very low RFP & high GFP when cereulide is below the threshold; very high RFP & low GFP when cereulide is above the threshold)?
      1. How?
        1. Manipulate the degradation rate of RFP or GFP
        2. By changing the ssrA tags sequence, there is different degradation performance[5]. We proposed an idea of collaborating with circuit design and wet lab to build a degradation tag library for more controllable degradation. With this collaboration, we can explore potential ssrA tags and utilise their performance to optimise our circuit. (Check Degradation tag library)
    2. Another problem that has come to our eyes is: Time. From the graph above, we can see that when cereulide is larger than the threshold, the RFP signal takes around 1.3 hours to be higher than the GFP signal. The time to produce an unsafe signal (result for users) is too long compared to the survey result from IHP (Acceptable waiting time of the test kits for users is around 45 mins).
      1. How?
        1. By also manipulating the ssrA tag sequence in both GFP and RFP, we could find the best pair that could shorten the result time, while having the 'sigmoid' output according to our input.
        2. Using optimization algorithms, we could find the optimised degradation rate. With the degradation tag library, we can pick the tag with similar degradation performance with our optimised rate. Another way is by machine learning (check degradation tag library-dry lab future computation work), it is possible that we could input the desired degradation rate and output the potential ssrA sequences for us.

    Future Direction and Improvement

    1. Experimental data validation for cycle 2
    2. Research and model possible circuit parts for speeding up time for signalling cascade in phosphorylation
    3. Changing promoter and RBS to replace PsinI and Phyperspank
    4. Develop Cycle 3 while developing degradation tag library
    5. Explore other possible parts/ mechanisms we can use to manipulate our circuit to give the desired output

    Discussions & Reflections on circuit modelling

    “All models are wrong, some are useful” - George Box

    Model is never accurate as there is always uncertainty or assumptions made due to complexity. At the beginning of constructing a model for our circuit, we encountered many issues such as unknown relationship of the kinetics and unknown parameters. Data is required to understand a relationship. Currently in this stage, we are still under research and gathering the data for better model constructing.

    How can our model be useful to the project in the future? In cycle 1, we found out that the results (RFP and GFP signals) are not sensitive to the input concentration. We then analyse our model by using Sensitivity Analysis to search for possible sensitive parameters for future development of our circuit. By using the result from SA, we proposed a dual degradation system for our circuit in order to achieve our goal: significant difference of the signal when the rice is safe to use or not. Before any physical prototyping, we verified our idea by doing modelling. By using the interpretation from the graph, we could see the possibility of controlling our threshold of showing red or green signal by using a dual degradation tag system. Though there must be differences when comparing to real life data, we see the potential of using modelling to further study the kinetics and interplay of different circuit components. In the future, we are going to design and build cycle 3 by playing with more parameters and circuit mechanisms to reach our goals.

    References

      [1] R. Naumann, D. Walz, S. M. Schiller, and W. Knoll, “Kinetics of valinomycin-mediated K+ ion transport through tethered bilayer lipid membranes,” Journal of Electroanalytical Chemistry, https://www.sciencedirect.com/science/article/pii/S0022072803000135 (accessed Oct. 12, 2023).
      [2] D. López, E. Gontang, and R. Kolter, “Potassium sensing histidine kinase in Bacillus subtilis,” Methods in Enzymology, https://www.sciencedirect.com/science/article/pii/S0076687910710132 (accessed Oct. 12, 2023).
      [3] Z. Chen et al., “Bacillus subtilis histidine kinase KinC activates biofilm formation by controlling heterogeneity of single-cell responses,” mBio, https://pubmed.ncbi.nlm.nih.gov/35012345/ (accessed Oct. 12, 2023).
      [4] K. L. Griffith and A. D. Grossman, “Inducible protein degradation in bacillus subtilis using heterologous peptide tags and adaptor proteins to target substrates to the protease ClpXP,” Molecular microbiology, https://pubmed.ncbi.nlm.nih.gov/18811726/ (accessed Oct. 12, 2023).
      [5] J. M. Flynn et al., “Overlapping recognition determinants within the ssrA degradation tag allow modulation of proteolysis,” The Proceedings of the National Academy of Sciences (PNAS), https://www.pnas.org/doi/full/10.1073/pnas.191375298 (accessed Oct. 12, 2023).