Modeling

Overview

Optopass system was intended to improve safety and security in synthetic biology. This makes highly accurate gene expression control essential. Therefore, it was important to evaluate and optimize our system quantitatively. We modeled the blue light-inducible promoter and the recombination system, which are fundamental parts of Optopass. We also modeled the Dummy System, an application of the basics of Optopass. In addition, we found a way to build differential equations from a network of "actors," and developed Genochemy, an educational tool that would allow for a realistic simulation of synthetic biology.

Modeling of Blue light inducible promoter

Purpose

The light-dependent promoter is the input part of Optopass and the very beginning of the grand system. Therefore, modeling this part is very important in evaluating the overall design and feasibility of Optopass. However, the experimental data from the wet lab does not provide sufficient information on the parameters, so simulations were performed and applied to the model for recombination.

Method and Modeling

EL222, which is used in blue light-activated promoters, undergoes a complex process of conformational change when exposed to blue light, forming homodimers and binding to DNA (see Blue light section in Description page for details), but the model employed in this study ([1]) simplifies this part of the process. In this model, we assumed that EL222 (simply TF), the transcription factor that activates the blue light promoter, has two states: on (TFon\mathrm{TF}_\mathrm{on}) and off (TFoff\mathrm{TF}_\mathrm{off}), and that EL222 in the state of TFoff\mathrm{TF}_\mathrm{off} will transition to TFon\mathrm{TF}_\mathrm{on} at a certain rate when blue light hits it. Since the rate constant is proportional to the intensity of blue light II, it can be written as IkonI\cdot k_\mathrm{on}, where konk_\mathrm{on} is a constant. Also, TFon\mathrm{TF}_\mathrm{on} changes to TFoff\mathrm{TF}_\mathrm{off} at a constant rate koffk_\mathrm{off} regardless of the light intensity. Putting these together, we can write the following differential equations.

dTFondt=IkonTFoffkoffTFondTFoffdt=IkonTFoff+koffTFon\frac{d\mathrm{TF}_\mathrm{on}}{dt}=I\cdot k_\mathrm{on}\mathrm{TF}_\mathrm{off}-k_\mathrm{off}\mathrm{TF}_\mathrm{on}\\ [0.5em] \frac{d\mathrm{TF}_\mathrm{off}}{dt}=-I\cdot k_\mathrm{on}\mathrm{TF}_\mathrm{off}+k_\mathrm{off}\mathrm{TF}_\mathrm{on}\\

Assuming that the sum of TF concentrations is constant, represented by the equation below,

TFtot=TFon+TFoff\mathrm{TF}_\mathrm{tot} = \mathrm{TF}_\mathrm{on}+\mathrm{TF}_\mathrm{off}

the differential equations can be rewritten as follows.

dTFondt=Ikon(TFtotTFon)koffTFon\frac{d\mathrm{TF}_\mathrm{on}}{dt}=I\cdot k_\mathrm{on}(\mathrm{TF}_\mathrm{tot}-\mathrm{TF}_\mathrm{on})-k_\mathrm{off}\mathrm{TF}_\mathrm{on}\\

Here, suppose that only the TFs in the on state can promote downstream gene expression. The mRNA of the downstream gene is transcribed in proportion to the Hill function with TFon\mathrm{TF}_\mathrm{on} as a variable and translated at a constant rate of ktransk_\mathrm{trans}. Also, assuming that mRNA is transcribed at a constant rate of kbasalk_\mathrm{basal} even if TFs in the on state is not present, and mRNA and protein are degraded at a rate of kdegRk_\mathrm{degR} and kdegPk_\mathrm{degP}, respectively, the following differential equation can be written.

dmRNAdt=kbasal+kmaxTFonnKdn+TFonnkdegRmRNAdProteindt=ktransmRNAkdegPProtein\frac{d\mathrm{mRNA}}{dt}=k_\mathrm{basal}+k_\mathrm{max}\frac{\mathrm{TF}_\mathrm{on}^n}{K_d^n+\mathrm{TF}_\mathrm{on}^n}-k_\mathrm{degR}\mathrm{mRNA}\\ [0.5em] \frac{d\mathrm{Protein}}{dt}=k_\mathrm{trans}\cdot\mathrm{mRNA}-k_\mathrm{degP}\cdot\mathrm{Protein}

Putting the previous equations together, the complete differential equation can be written as below.

dTFondt=Ikon(TFtotTFon)koffTFondmRNAdt=kbasal+kmaxTFonnKdn+TFonnkdegRmRNAdProteindt=ktransmRNAkdegPProtein\frac{d\mathrm{TF}_\mathrm{on}}{dt}=I\cdot k_\mathrm{on}(\mathrm{TF}_\mathrm{tot}-\mathrm{TF}_\mathrm{on})-k_\mathrm{off}\mathrm{TF}_\mathrm{on}\\ [0.5em] \frac{d\mathrm{mRNA}}{dt}=k_\mathrm{basal}+k_\mathrm{max}\frac{\mathrm{TF}_\mathrm{on}^n}{K_d^n+\mathrm{TF}_\mathrm{on}^n}-k_\mathrm{degR}\cdot\mathrm{mRNA}\\ [0.5em] \frac{d\mathrm{Protein}}{dt}=k_\mathrm{trans}\cdot\mathrm{mRNA}-k_\mathrm{degP}\cdot\mathrm{Protein}

The values of each parameter are as follows. Basically, the values in the paper [1] were used as they are, some of them appropriately multiplied by the unit (labeled "converted").

Table 1. Set of parameters.

DescriptionValueCited from
TFtot\mathrm{TF}_\mathrm{tot}total cellular TF1.142×103nM1.142\times 10^3\mathrm{nM}[1], converted
konk_\mathrm{on}light dependant TF activation rate1.6399×103s1(μWcm2)11.6399\times 10^{-3}\mathrm{s^{-1}(\mu Wcm^{-2})^{-1}}[1]
koffk_\mathrm{off}TF dark-state reversion rate3.4393×101s13.4393\times 10^{-1}\mathrm{s^{-1}}[1]
kbasalk_\mathrm{basal}basal transcription rate1.491×102nMs11.491\times 10^{-2}\mathrm{nMs^{-1}}[1], converted
kmaxk_\mathrm{max}maximal induced transcription rate7.7565nMs17.7565\mathrm{nMs^{-1}}[1], converted
ktransk_\mathrm{trans}translation rate1.4514s11.4514\mathrm{s^{-1}}[1]
kdegRk_\mathrm{degR}mRNA degradation rate4.2116×102s14.2116\times 10^{-2}\mathrm{s^{-1}}[1]
kdegPk_\mathrm{degP}protein degradation rate1.166×104s11.166\times 10^{-4}\mathrm{s^{-1}}[1]
KdK_dTFon\mathrm{TF}_\mathrm{on} level required for achieving kmax/2k_\mathrm{max} / 25.4614×104nM5.4614\times10^{-4}\mathrm{nM}[1], converted
nnhill coefficient4.2034.203[1]

We employed this model to simulate a blue light-dependent promoter.

Results

Figure 1 shows a comparison of continuous exposure to blue light at a constant intensity(I=200I=200) and no exposure to blue light. The continuous exposure to blue light initially increases the protein concentration, but it saturates after about 600 minutes, after which it remains constant. The amount of protein expressed in the absence of light is about 73nM73 \mathrm{nM}, which is about (1300\frac{1}{300}) of that in the presence of light.

Figure 1

Figure 1. Change in protein concentration over time with continuous exposure to blue light (left) and without exposure to blue light (right).

We also simulated what happens to the protein expression when the cells are repeatedly exposed to blue light and not exposed to blue light. Protein expression was simulated in 5 different types of cycles. The cycle time was set to 30 minutes, and the duty ratio, defined as the percentage of time the light was on, was

  • 0% = always off
  • 25% = 7.5 minutes on, 22.5 minutes off cycle
  • 50% = 15 minutes on, 15 minutes off cycle
  • 75% = 22.5 minutes on, 7.5 minutes off cycle
  • 100% = always on

The results showed that the maximum protein expression increased in proportion to the duty ratio.

Figure 2

Figure 2. Time variation of protein concentration at different duty ratios.

Modeling of Recombination

Purpose

In our project, we use a site-specific recombinase system to control yeasts by the order of the light. To achieve this, we needed to quantitatively comprehend the circumstances in which the recombination takes place. For further details on our recombinase system, please see the Why Order section on Description page.

In Optopass, recombinase is expressed at the downstream gene of the light-inducible promoter, but a little amount of it can still be expressed due to leakage even when the promoter is not activated. Therefore, the recombination may take place also in the dark if it occurs even with little amounts of recombinase, which could lead to the unexpected disarming of the security mechanism.

This modeling tries to improve the project by determining at what rate recombination happens under various conditions. Specifically, assuming leakage, we assessed how fast the recombination occurs at low recombinase concentrations and proposed ways to make the system more leakage tolerant. In addition, we modeled how rapidly the recombination takes place at a high concentration of recombinase, to estimate the time it takes for the recombination to occur when the light activates the promoter. In order to reduce the time for each light stimulus, we found a way from simulation to increase the likelihood of the recombination, since the time it takes for recombinase to process the recombination for most individuals is the basis for the time required to shine the light. We employed the worm-like chain model developed by the iGEM Chalmers-Gothenburg team, our partner team for this year.

As a result of our model, it was found that the system will be more resilient to the promoter leakage by altering the recombinase binding sequence to the one which is less likely to bind to the recombinase or by using Split-Cre. Additionally, the worm-like chain model demonstrated that, in situations where recombinase is expressed at high concentrations, the probability of the recombination might be improved by decreasing the length of sequences between the recognition sequences.

Preparation

Consider a linear DNA strand with two Recombinase Recognition Sequences in the same orientation. Two recombinases will eventually bind to each recognition sequence.

Figure 1

Figure 1. The assumed model of the binding of the recombinase to the DNA strand.

The steps in the process of recombinase binding to the two recognition sequences on the DNA strand, recombinase association, and recombination are listed with reference to [2].

Figure 2

Figure 2. Reaction steps before and after recombination, where SiS_i for each ii represents the states before recombination, where the recombinase binds to the linear DNA strand; II represents the state where the two recombinase binding sites meet to form a junction; P0P_0 represents the state after recombination, where one loop of DNA and one linear DNA are obtained. The states of recombinase binding after P0P_0, as in the SiS_i steps, are assumed but omitted in this figure.

In vitro studies on the degree of recombination reaction progression according to the amount of recombinase are common, and many of those studies have simulated the DNA strand using its concentration ([2][3]). While many of iGEM's previous projects have modeled recombinases in this manner using the concentration of DNA sequences ([4][5]), we assume that the designed sequences will be incorporated into the genome, which means that there will be one DNA strand, or at most two. In this situation, it would be more appropriate to reconsider each step of the reaction process as a "state" that a single DNA sequence can take, rather than modeling the DNA sequence as a concentration.

We constructed a model which consists of the stochastic state of a single DNA strand, referring to the model which uses DNA concentration to determine the parameters. As an example, the differential equations for the concentration of DNA in the state of S0S_0 and S1S_1 in Figure 2 are as follows:

d[S0]dt=k1[S0][E]+k1[S1]d[S1]dt=k1[S0][E]k1[S1]+\frac{d[S_0]}{dt} = -k_1[S_0][E]+k_{-1}[S_1] \\ [0.5em] \frac{d[S_1]}{dt} = k_1[S_0][E]-k_{-1}[S_1]+\cdots

[E][E] represents the concentration of recombinase. The \cdots in d[S1]dt\frac{d[S_1]}{dt} are the terms which represent the reaction where S1S_1 becomes S21S_{21} or S22S_{22}. We omitted these terms here because we wanted to focus on the reaction between S0S_0 and S1S_1.

From this, we can obtain the probability for "a DNA in the S0S_0 state to transition to S1S_1 state in Δt\Delta t", which we named p1p_1 and the probability for "a DNA in the S1S_1 state to transition to S0S_0 in Δt\Delta t", which we named p1p_{-1}. The probabilities were calculated from the proportions, that is, we calculated p1p_1 as "the proportion of DNAs in the S0S_0 state that transitions to the S1S_1 state in Δt\Delta t". p1p_{-1} was calculated in the same way.

p1=k1[S0][E]Δt[S0]=k1[E]Δtp1=k1[S1]Δt[S1]=k1Δtp_1 = \frac{k_1[S_0][E]\Delta t}{[S_0]} = k_1[E]\Delta t \\ [0.5em] p_{-1} = \frac{k_{-1}[S_1]\Delta t}{[S_1]} = k_{-1}\Delta t

It is needed to adjust the Δt\Delta t so that pip_i satisfies 0pi10 \leq p_i \leq 1.

Thereafter, since the number of molecules of recombinase is very large relative to the number of molecules of DNA sequence, we assume that the concentration of recombinase does not change before and after the reaction. In other words, we assume that [E][E] is constant as long as there are no external influences.

Assumption 1.

Recombinase concentration [E][E] is assumed to be constant before and after the reaction.

Detailed explanation of Assumption 1.

Click to open details

Many reactions can be expressed as "a recombinase attaches to or detaches from a DNA sequence," but in the recombination part, the reaction "a linear DNA sequence becomes two DNA sequences (a loop and a linear sequence with the loop part skipped)" and the reverse reaction, "a linear DNA sequence is obtained from two DNA sequences" is exceptional. The former (positive reaction) can be regarded as "the DNA in state II transitions to state PP where two DNA sequences (PaP_a and PbP_b) coexist" as in the above, but the latter (reverse reaction) is a bimolecular association reaction of DNA, and information on probability cannot be easily obtained from concentration parameters as in the above. However, the reverse reaction is hardly necessary to consider in situations where the DNA concentration is very low. This point was ensured by interviewing a specialist in recombinase in yeast (See Integrated Human Practices page for details). Based on these considerations, the reverse reaction of recombination is not considered here.

Assumption 2.

Assume that the reverse reaction of recombination does not occur.

Detail: Forward and backward reactions of recombination.

Click to open details

At this point, we are ready to model by probability. From now on, we will describe the model based on probability instead of concentration.

The transition probabilities between states are expressed using eight probabilities: p1p_1, p1p_{-1}, p2p_2, p2p_{-2}, p3p_3, p3p_{-3}, p4p_4, and p4p_{-4}, described in Figure 3 below.

Figure 3

Figure 3. Transition probabilities between states. Post-recombination states are now consolidated to PP from Assumption 2.

The state after recombination is now denoted by PP.

The following table shows the meaning of each probability and its calculation formula based on the concentration parameters.

Table 1. Details of each probability.

DescriptionFormula (using parameters in Table 2)
p1p_1Probability that a transition in which a recombinase newly binds to recognition sequence with no recombinase occurs in Δt\Delta tp1=k1[E]Δtp_1=k_1[E]\Delta t
p1p_{-1}Probability that a reverse transition for transition corresponding p1p_1 occurs in Δt\Delta tp1=k1Δtp_{-1}=k_{-1}\Delta t
p2p_2Probability that a transition that a recombinase binds to recognition sequence with one recombinase occurs in Δt\Delta tp2=k2[E]Δtp_2=k_2[E]\Delta t
p2p_{-2}Probability that a reverse transition for transition corresponding p2p_2 occurs in Δt\Delta tp2=k2Δtp_{-2}=k_{-2}\Delta t
p3p_3Probability that S4S_4 transitions to II in Δt\Delta tp3=k3Δtp_3=k_3\Delta t
p3p_{-3}Probability that II transitions to S4S_4 in Δt\Delta tp3=k3Δtp_{-3}=k_{-3}\Delta t
p4p_4Probability that II transitions to PP in Δt\Delta tp4=k4Δtp_4=k_4\Delta t
p4p_{-4}Probability that PP transitions to II in Δt\Delta tp4=0p_{-4}=0 (from Assumption 2)

The parameters appearing in the formulas in Table 1 were not obtained from our experimental results, so we referred to literature values obtained in in vitro systems [2][3].

Table 2. Parameters from literature.

ValueCited from
k1k_12.2×101(nMs)12.2\times10^{-1}\mathrm{(nM\cdot s)^{-1}}[3]
k1k_{-1}6.6×102s16.6\times10^{-2}\mathrm{s^{-1}}[3]
k2k_22.3×101(nMs)12.3\times10^{-1}\mathrm{ (nM\cdot s)^{-1}}[3]
k2k_{-2}4.8×103s14.8\times10^{-3}\mathrm{s^{-1}}[3]
k3k_31.23×103s11.23\times10^{-3}\mathrm{s^{-1}}Value of 3044 bp in [2]. min1\mathrm{min^{-1}} was rewritten as s1\mathrm{s^{-1}}.
k3k_{-3}1.15×101s11.15\times10^{-1}\mathrm{s^{-1}}Average (see Assumption 3) of the 3044 bp and 870 bp values in [2], where min1\mathrm{min^{-1}} is rewritten as s1\mathrm{s^{-1}}.
k4k_42.92×104s12.92\times10^{-4}\mathrm{s^{-1}}Average (see Assumption 3) of the 3044 bp and 870 bp values of k4k_{-4} in [2] (for clarity, the direction of inversing was reversed from [2]), where min1\mathrm{min^{-1}} is rewritten as s1\mathrm{s^{-1}}.

k3k_3 represents the value at 3044 bp, and k3k_{-3} and k4k_4 is the average of the value at 3044 bp and 870 bp. This point will be discussed later (at (ii) Increasing the maximum probability).

Simulation

Using these parameters, simulations were performed up to 12 hours in the situation where recombinase concentration is 100 nM (roughly assumed a situation where recombinase is expressed at high levels).

The legend of "Product" represents the DNA sequence after recombination (state PP). As time passes, the recombination probability increases. In other words, the "probability of being the post-recombination DNA sequence" increases.

Figure 4

Figure 4. the states of each state up to 12 hours later. Product represents the state PP. Intermediate represents state II. Substrate includes all the states SiS_i for each ii.

We actually simulated the state transition of each DNA strand at each probability. The figure below shows the results. It was found that DNA was in states S3S_3 and S4S_4 for much of the time. As a result, the recombination actually occurred (reached state P) in 17% of all DNA strands by 12 h, which is certainly close to the result in Figure 4.

Figure 5

Figure 5. Simulation of the transition between states of DNA strand (n = 100). Left: The typical transition of a DNA strand. Recombination did not occur until 12 hours so state PP is not shown. Right: All transition (n = 100).

We were interested in the probability of recombination after a certain time at a certain concentration of recombinase. We obtained the following graph.

Figure 6

Figure 6. Recombination probability after 12 hours for each constant recombinase concentration. Horizontal axis is recombinase concentration (constant, see Assumption 1) and the vertical axis is the recombination probability after 12 hours.

The graph shows that the recombination probability is low below a certain threshold and converges to a constant value above it.

There are two points where we would like to modify the parameters in this model.

Figure 7

Figure 7. Explanation of Figure 6.

(1) Increase the threshold value. This increases the tolerance to leakage. When using the parameters in the paper without optimizing, the threshold value was around 0.1 nM, which is quite severe considering the expression level of the light-dependent promoter in the dark as described earlier. This may be due to the fact that the experimental environment in the paper was in vitro, which is very different from the actual in vivo situation, but our Results suggest that the recombinase leakage happened in vivo. It is thought that increasing the threshold value will be important in the actual in vivo situation as well.

In this part, the measurements of our experiments led us to the decision of continuing modeling with these parameters.

Alternatively, it may be that leakage is mitigated by reducing the number of active recombinases in the first place. Since this is also essentially similar to the threshold, we will consider this together with (1).

(2) Increase the maximum probability. When the recombinase concentration is greater than the threshold, the probability of recombination converges to a certain probability no matter how high the concentration is. 12 hours later, when the recombinase is legitimately expressed under the correct light, we would like the probability of recombination to be higher, in order to shorten the time of light exposure and make the system easier to use. As in (1), the experimental environment in the paper is very different from that in vivo, but we believe that the following measures to increase the maximum probability will play a role to some extent in vivo as well.

(1) Increase threshold and reduce the amount of active recombinase

The results of the simulation using the blue light promoter described earlier indicate that about 73 nM of recombinase would be expressed at leakage. We would like to be able to set the threshold greater than 73 nM by combining some ways.

As an approach to increase the threshold value, we found the following measure to be effective as a result of simulation from several simulations.

  • Make p1p_1 and p2p_2 smaller. In other words, make it harder for the recombinase to attach to the DNA sequence. For this purpose, the recombinase recognition sequence may be changed.

Another possible method to reduce the amount of active (have recombination ability) recombinase is as follows.

  • Attach a degradation tag to the recombinase to accelerate degradation.
  • Modify the mRNA of the recombinase (e.g. attach something like a riboswitch on the 5' side) to reduce the amount of translation.
  • Light-dependent Split Cre can be used to reduce the amount of recombinase that are active (has recombination ability) in the absence of light.

The model for using degradation tag was omitted because it is common. The method using Split Cre is quite practical and unique, so this will be discussed later.

p1p_1 and p2p_2 are the probabilities of the transition of the recombinase binding to DNA. It is certainly conceivable that making these probabilities smaller will make it more difficult for the recombinase to bind to DNA, thereby increasing the threshold at which the recombination probability improves.

Consider the ease of binding of a recombinase to a recognition sequence to be varied by a constant ratio. That is, using the constant α\alpha, we assume that new p1p'_1 and p2p'_2 can be defined by

p1=αp1p2=αp2.p'_1=\alpha p_1 \\ [0.5em] p'_2=\alpha p_2.

p1,p2p_{-1},p_{-2} are assumed to be unchanged here. Since p1p1p_1 \gg p_{-1} and p2p2p_2 \gg p_{-2} in the situation we are considering here, we assume that the changes in p1p_{-1} and p2p_{-2} are negligible compared to those of p1p_1 and p2p_2.

We considered that p3p_3 and p4p_4 would also remain unchanged. Though the probability of recombination may change due to the change of recognition sequences, we assumed that the effect of mutation of the recognition sequences on these parameters would be small. The simulation was performed.

Figure 8

Figure 8. Recombination probability for each recombinase concentration when α\alpha is changed.

While the maximum probability has not changed, we could certainly shift the threshold to the right, indicating an increased tolerance for leaks.

The value of α\alpha could be adjusted by using a recognition sequence that is less likely to bind to the recombinase than loxP with respect to Cre recombinase. For example, changing loxP to loxAA would reduce the final recombination probability by about 1/2 ([6]). Although the information in the paper [6] did not provide information on the specific value of α\alpha in the loxP variant, it can be assumed that only α\alpha, the ease of binding between the recombinase and the recognition sequence, has changed because the other conditions are assumed to be unchanged. Then, in this case, as seen in the figure above, the threshold could be shifted.

The α\alpha in Figure 8 is still not sufficient to deal with actual leaks. Of course, the experimental environment may be very different from that of actual organisms, but we can deal with the problem by combining some ways, including the use of PA-Cre.

PA-Cre Modeling

PA-Cre is a Cre with recombination activity only when exposed to blue light. It is composed of an N-terminal domain of Cre, CreN, to which a photoswitch called nMag is attached, and a C-terminal domain CreC to which a photoswitch called pMag is attached. These two domains are separated in the dark, but when exposed to blue light, nMag and pMag bind to each other, resulting in the binding of the two complexes and the acquisition of recombinase activity [7].

Leakage can be reduced by expressing PA-Cre instead of Cre in Optopass. Due to the basal expression of promoters and the wavelength characteristics of LED lights, it is difficult to reduce promoter expression to strictly zero when not exposed to blue light. Considering the fact that recombination can occur even at low levels of recombinase, it is essential to reduce the expression of active Cre when not exposed to blue light. Replacing Cre with PA-Cre reduces the recombinase activity when blue light is not shone, since PA-Cre does not have recombinase activity until it is activated by further exposure to blue light after being expressed.

Modeling on PA-Cre was performed. We used a combination of papers on modeling EL222, a blue light promoter [1], and PA-Cre [8].

Details of modeling

Click to open details

The results of the simulation using this model are as follows.

First, there is not much difference between Cre and PA-Cre when exposed to light of sufficient intensity (I=200[μW/cm2]I=200 [\mathrm{\mu W/cm^2}]), indicating that recombination can occur sufficiently with PA-Cre as well. (Figure 9)

Figure 9

Figure 9. Change in concentration of Cre and PA-Cre when exposed to light of normal intensity (I=200I=200).

Then, we simulated the case when about 1% intensity of the blue light (I=2I=2) hits the yeast. This may occur when it is illuminated with the light from a green LED, due to the characteristics of the LEDs used [9].The result is shown in Figure 10. The concentration of Cre is converged to about 73 nM, as described above, but that of PA-Cre is converged to about 4.9 nM, which is about 6.7% of that of Cre.

Figure 10

Figure 10. Changes in Cre and PA-Cre concentrations under weak light (I=2I=2).

In this way, the use of PA-Cre can reduce the concentration of active recombinase expressed in the absence of blue light.

(2) Increase the maximum probability

For the maximum probability, we tried several simulations and found that the influence of the values of p3p_3, p3p_{-3} and p4p_4 was large. This is confirmed by the fact that the stochastic model in Figure 5 often stays in state S4S_4 and II, and these transitions are the rate-limiting factors of the overall recombination. Simulations show that changes in the value of these probabilities contribute little to the threshold in (1). This made it possible to perform the optimization of (2) independently of (1).

The paper [2] shows that k3k_3, k3k_{-3}, and k4k_4 (note: k4k_{-4} in the paper's orientation), corresponding to p3p_3, p3p_{-3}, and p4p_4, depend on the sequence length between recognition sequences.

k3k_3 is "the rate constant for the reaction in which the recombinases attached to the two recognition sequences of a single DNA strand meet to form a complex". Since the paper only measured the sequence length between recognition sequences in two cases, we continuously simulated the dependence of k3k_3 on the sequence length between recognition sequences using the "worm-like chain model" used by Chalmers-Gothenburg, a partnership team of us.

Since we could not clarify how k3k_{-3} and k4k_4 are affected by the sequence length between two recognition sequences, we assumed that they are constant regardless of the sequence length. Therefore, we took the average of k3k_{-3} and k4k_4 at each sequence length in the paper.

The paper's values are denoted with '. Since k4k_4 corresponds to k4k'_{-4} in the paper, we also took that into account.

k3=12(k3(870bp)+k3(3044bp))k4=12(k4(870bp)+k4(3044bp))k_{-3}=\frac{1}{2}\left(k'_{-3(\rm{870 bp})}+k'_{-3(\rm{3044 bp})}\right) \\ [0.5em] k_4=\frac{1}{2}\left(k'_{-4(\rm{870 bp})}+k'_{-4(\rm{3044 bp})}\right)
Assumption 3.

Assumed k3k_{-3} and k4k_4 are constant regardless of sequence length and averaged from the data in the paper.

Worm-like chain model

Worm-like chain model treats DNA as a wire with a flexibility depending on its bending constant [10]. This model makes it possible to calculate the possibility that a DNA takes a certain form. By calculating how the DNA chain is deformed, we can calculate the probability that the DNA takes a shape such that "recombinases attached to two recognition sequences on a single DNA strand" come within a certain distance.

For the p3p_3 transition (from S4S_4 to II), the recombinases attached to the two recognition sequences must come close together. We assumed that once they got close, the reaction to form the complex would proceed. Then, the probability p3p_3 of the transition from the state S4S_4 to the state II can be written as p3=βq(N)Δtp_3=\beta\cdot q(N)\Delta t using the probability q(N)q(N) of the recombinases attached to the two recognition sequences approaching within a certain distance and the constant β\beta. NN is the sequence length between the recognition sequences.

Figure 11

Figure 11. Relationship between the length of sequences between recognition sequences and the probability that two recombinases come close to each other (both logarithmic).

The results of the simulation using the worm-like chain model showed that q(N)q(N) is proportional to the length of the sequence between the recognition sequences, NN, to the power of approximately -1.3814. Thus we assumed that q(N)q(N) can be expressed as follows:

q(N)=aN1.3814q(N)=aN^{-1.3814}

Therefore,

p3=aβN1.3814Δtp_3=a\beta N^{-1.3814}\Delta t

Fitting to the k3k_3 literature values for the two NN, we obtained aβ=448.3761032[s1]a\beta=448.3761032\mathrm{[s^{-1}]}. Simulating recombination while changing the sequence length, the recombination probability was found to be as follows. The recombinase concentration is adjusted to be above a threshold value (i.e., a concentration that produces a value close to the maximum probability of convergence destination in the previous figure).

Figure 12

Figure 12. Relationship between sequence length between recognition sequences and probability of recombination after 12 hours.

We first simulated in Figure 6 using the parameter k3k_3 for the 3044 bp sequence length. The data in the paper is only 3044 bp and 870 bp, but we need to put the recombinase and terminator between the recognition sequences, which requires about 1250 bp at minimum.

The graph below is the comparison between 3044 bp and 1250 bp.

Figure 13

Figure 13. The comparison of recombination after 12 hours between 3044 bp and 1250 bp. p3p_3 is now calculated from p3=aβN1.3814Δtp_3=a\beta N^{-1.3814}\Delta t, so the probability of 3044 bp is changed from Figure 6.

The worm-like chain model allows us to calculate the probability of recombination with various sequence lengths. This allowed us to learn that we can increase the probability of recombination by a factor of 1.75 in the case of 1250 bp, for example, compared to 3044 bp.

Simulation using optimized parameters

Simulation was performed using the parameters after the threshold was raised by the optimization in (1) and the maximum probability was raised by setting the sequence length between recognition sites to 1250 bp in (2). The parameter α\alpha is set to be 0.01, although it is too small compared to Figure 8.

Figure 14

Figure 14. Simulation with optimized parameters. The point shows the probability at Recombinase is 4.9 nM.

When PA-Cre is used, the active recombinase concentration is assumed to be 4.9 nM at leakage. In this case, the concentration is below the threshold and recombination is less likely to occur (about 6%, from Figure 14). Moreover, if the light hits, as shown in Figure 9, PA-Cre averages more than 3000 nM, which is far above the threshold, and we can say that there is more than an 80% probability of recombination.

Conclusion

In this modeling, we confirmed that the recombination threshold can be raised to be more tolerant to leakage by reducing the α\alpha, which represents the ease of recombinase attachment to the recognition sequence, and by reducing the amount of active recombinase using PA-Cre. This ensures that security can be improved by reducing the possibility of being deciphered in the dark. We also confirmed that the maximum recombination probability can be increased by shortening the length of sequences between recognition sites to 1250 bp. This is expected to shorten the time it takes for recombination to occur and increase the convenience of the system.

However, there are still some unresolved issues. Although α\alpha is set to 0.01 for convenience, there are not always such recognition sequences in reality. In addition, we have not been able to investigate in detail how much the recombination probability of 6% at leakage, which was obtained in Figure 14, would have an impact. For example, when using the Dummy System described in the Description page (also described below), it will be necessary to verify in detail whether a 6% probability of accidental release of the killer factor will kill the target completely. Furthermore, here we focused only on leaks in the dark, but there is also the problem of recombinase being expressed by exposure to light of another color. But it can be solved in the same way as we did here, as long as we know how much recombinase will be expressed at that time.

The parameters used in this model were obtained from in vitro experiments. Although more accurate modeling requires simulations with in vivo parameters, our method, which was able to improve relative thresholds and probabilities, is likely to be effective even with different parameters.

Modeling of the Dummy System

Abstract

The Dummy System is a new containment technique based on the fundamental concept of Optopass. It exploits the structure of a dynamic system in which two microorganisms coexist in a device in a competitive relationship and their equilibrium is shifted by external inputs. Here, a microorganism carrying the gene of interest (target) coexists with a corresponding GMO (dummy) designed to respond to the order of light exposure. If the light is cast in the correct order, the dummy will kill itself, but if the light is cast in the wrong order or it is exposed to sunlight, it will kill the target organism, using killer factors. As a result, the target microorganism can be obtained from the device only if the light code is entered in the correct order, while only the dummy can be obtained if the device is destroyed or the wrong code is entered. Thereby, the target microorganism can be protected.

Purpose

We came up with a new idea to protect strains using the basic concept of Optopass through the discussion with experts we held with Tohoku Medical Megabank Organization (see Integrated Human Practices page). We call this "the Dummy System."

Figure 1

Figure 1. Schematic diagram of the Dummy System.

When transporting genetically modified strains with special or dangerous characteristics, it is necessary to avoid leakage of the strains and their sequences, particularly in terms of maintaining confidentiality.

To fulfill the purpose, the functions required in the Dummy System are as follows:

  1. In the initial state, a mixture of the target and dummy microorganisms is cultured, making it nearly impossible to read the sequence of the target even with DNA sequencing.
  2. When the cryptic input is made, the balance of the mixtures in the system is disrupted, and if the correct input is made, the target will occupy the entire system, while if the wrong input is made, the dummy will occupy the entire system and the target will not be obtainable.

To achieve these characteristics, we have conducted the following examination.

First, we verified whether it is possible for two species to coexist in the same environment through research and simulation, and searched for parameters in the differential equations representing the competition between the two species as a means to achieve this goal.

Second, we evaluated both of the situations where the correct light is applied and where the wrong light or sunlight is applied, using the parameters obtained in the first examination. We verified that the cut off of the Sod1 gene and the transcription of the killer gene would cause changes that would make the target and the dummy dominant, respectively.

Method and Modeling

Coexistence conditions

When a strain is cultivated alone from a very small population, the change in the population follows a sigmoid curve and can be represented by the following logistic equation [14]:

dNdt=rN(1NK)\frac{dN}{dt}=rN\left(1-\frac{N}{K}\right)

where NN stands for population, tt for time, rr for proliferation rate, and KK for environmental carrying capacity.

Next, we present a model to simulate the coexistence of the dummy and target strains. When simulating the population growth of two strains in a competitive relationship, we need to consider the influence both inside and between them. That means, the growth of each strain is hindered by the combination of both the density effect within the strain and the negative effect from the other strain. These complex dynamics can be simply modeled with the help of the Competitive Lotka-Volterra’s equations.

Taking N1N_1 as the population of microorganism of the dummy strain and N2N_2 as the population of that of the target strain, the Competitive Lotka-Volterra’s equations are obtained as the following simultaneous differential equations:

dN1dt=r1N1(1N1+α12N2K1)dN2dt=r2N2(1N2+α21N1K2)\frac{dN_1}{dt}=r_1N_1\left(1-\frac{N_1+\alpha_{12}N_2}{K_1}\right) \\ [0.5em] \frac{dN_2}{dt}=r_2N_2\left(1-\frac{N_2+\alpha_{21}N_1}{K_2}\right)

where r1r_1, and r2r_2 stand for proliferation rates, K1K_1 and K2K_2 for environmental carrying capacities, and α\alpha for competition coefficient.

After a certain period of time, the population of each species converges to an equilibrium point.

The equilibrium point is obtained by solving the simultaneous differential equations with the isocline method.

Details of the isocline method

Click to open details

In Gause's paper [11], it is said that Saccharomyces cerevisiae and Schizosaccharomyces cerevisiae coexist in the same medium. Considering the coexistence of these two species as a model of the coexisting state of the Dummy System, the former can be regarded as a dummy strain and the latter as a target strain. The necessary parameters for the Lotka-Volterra's competition equation were calculated from the experimental data shown in the paper. They were determined by fitting the theoretical simulation results with the least squares method. To learn about the details of this fitting, see the explanation below.

Fitting by the least-squares method

Click to open details

Table 1. The parameter

ParameterValue
r1[/h]r_1\mathrm{[/h]}0.110649
r2[/h]r_2\mathrm{[/h]}0.020761
K1K_111.318389
K2K_25.455101
α12\alpha_{12}0.929203
α21\alpha_{21}0.117767

In the simulation, we assumed that the total population of the initial Target and Dummy were 1 and that they would mix in the proportion of 1:1, so initial values were set to 0.5 and 0.5. We used the Euler method for the numerical computation of differential equations. The calculations were performed in a Python program using the SciPy package.

Here are the simulation results based on parameters we got through optimization.

Figure 3

Figure 3. The variation in the population until the Dummy and Target reach an equilibrium point when the initial size is (0.5,0.5) (relative values). The plots indicate the experimental data obtained from Gause’s paper.

Here, we confirmed that equilibrium between two species of yeast is achievable. Gause (1934) [11] did an experiment about the coexistence of them, but it was insufficient in that he only collected data at the beginning of translation, and did not collect the data until they reached equilibrium or one of them eliminated the other. Looking at the concrete transition, the dummy with the higher proliferation rate increases first but decreases later on and gradually moves toward equilibrium. Such changes are often seen as variations in the population of two species, which suggests that proper fitting has been done. It is desirable that there are very few targets compared to the dummy at the equilibrium, because in this way, the target sequence could be difficult to read off even if the cocultured system is inappropriately taken out maintaining the first coexistence state. We could not realize this in this simulation, which would be a next challenge. At least we have succeeded in realizing equilibrium itself.

Does this result show the suitability to the social implementation?

Click to open details

When light is applied in the correct order

Next, we examined the shift in this coexistent state when exposed to the correct passcode. In the Dummy System, we realize the target occupying the entire system with the correct light input by cutting off Sod1, one of the tetracycline resistance genes [12][13], in the dummy to eliminate the dummies.

Since tetracycline resistance genes are present in yeasts [15], both strains can grow in tetracycline medium although tetracycline inhibits protein synthesis. However, when the correct input is applied, Sod1 in the dummy is cut off by the Cre-loxP system, and the dummies would not be able to grow anymore and are competitively eliminated.

Here, we continued the simulation after the light exposure by modifying the competitive Lotka-Volterra’s equations with some assumptions and settings.

The competition equation of Lotka-Volterra, shown before, is as follows:

dN1dt=r1N1(1N1+α12N2K1)dN2dt=r2N2(1N2+α21N1K2)\frac{dN_1}{dt}=r_1N_1\left(1-\frac{N_1+\alpha_{12}N_2}{K_1}\right) \\[0.5em] \frac{dN_2}{dt}=r_2N_2\left(1-\frac{N_2+\alpha_{21}N_1}{K_2}\right)

This can be rewritten as follows by arranging the coefficients:

dN1dt=N1(r1β1N1γ12N2)dN2dt=N2(r2β2N2γ21N1)\frac{dN_1}{dt}=N_1(r'_1-\beta_1N_1-\gamma_{12}N_2) \\ [0.5em] \frac{dN_2}{dt}=N_2(r'_2-\beta_2N_2-\gamma_{21}N_1)

Figure 3 shows that microorganisms are in an equilibrium state after t=400t=400, so t=400t=400 was used as the time of correct input.

When the correct input is applied, Sod1 would be cut off. Here we can say that r1=0r'_1=0 since it cannot multiply any more. For simplicity, the following assumption was made.

Assumption 1.

Sod1 was cut off immediately after the correct light was cast.

The setting of r1=0r'_1=0 at t=400t=400 does not consider the possibility of some dummy cells to still proliferate after tetracycline starts to affect them. However, we judged that this assumption was acceptable because here we do not have to seriously discuss the speed at which the dummy disappears. It was sufficient to confirm that the dummy eventually dies when obtaining the target.

Figure 4

Figure 4. Population changes when exposed to the correct light at t=400t=400 (relative values).

Figure 5

Figure 5. Enlarged view around t=400t=400 in Figure 4 (relative values).

Figure 6

Figure 6. Logarithmic changes in the ratio of the population of the target to that of the dummy when the correct light was applied.

Figure 6 shows a linear increase after t=400t=400. Two days after the irradiation, the amount of target is about e2.512e^{2.5} \approx 12 times larger than the dummy, which would be sufficient to use the target. This shows our system is practical because the time required isn’t so long compared to the time needed when yeast is cultured by itself.

When light is applied in the wrong order

Finally, we examined the shift in this coexistent state when exposed to the wrong passcode. In the Dummy System, when the light is applied in the wrong order, the dummy strain occupies the whole system, by releasing K28, one of the killer factors which derives from the so-called killer yeast. The killer factor is a protein that kills sensitive yeasts that do not have resistance to it, and killer yeasts have resistance to the factor of their own. In our system, K28 kills the target strain by inducing apoptosis, while the dummy strain is not affected because it is equipped with the resistance gene the killer yeast has. Since K28 destroys DNA into fragments, sequence information cannot be obtained from the dead target yeast.

Modeling in previous studies [16] has confirmed that the killer factor K28 reduces the target population exponentially, so the following differential equation was formulated.

dN1dt=r1N1(1N1+α12N2K1)dN2dt=r2N2(1N2+α21N1K2)rkillN2\frac{dN_1}{dt}=r_1N_1\left(1-\frac{N_1+\alpha_{12}N_2}{K_1}\right) \\ [0.5em] \frac{dN_2}{dt}=r_2N_2\left(1-\frac{N_2+\alpha_{21}N_1}{K_2}\right) - r_\mathrm{kill}N_2

Parameters used before, such as r1r_1, r2r_2, K1K_1, K1K_1, α12\alpha_{12}, and α21\alpha_{21} are the same as above. This is because the competition between the two strains is expected to remain unchanged even after K28 is released. The only difference is that a term representing apoptosis by the killer factor is added at the end of the differential equation of the N2N_2 variable.

The effect of the killer factor was assumed to be 1.7%/min reduction rate i.e. 64%/hour reduction rate, using the reference. [16] Thus, rkill=0(t<400),0.64[/h](t400)r_\mathrm{kill}=0 (t < 400), 0.64 \mathrm{[/h]} (t\geq 400).

Figure 7

Figure 7. Population changes when exposed to the wrong light at t=400t=400 (relative values).

Figure 8

Figure 8. Enlarged view around t=400t=400, when light is cast in Figure 7 (relative values).

Figure 9

Figure 9. Changes in the ratio of the population of the target to the dummy when the wrong light was applied.

Since the target is nearly eliminated (to 1/40 of the dummy) 5 hours after exposure to the wrong light, it seems to be difficult to steal the sequence of the target. Moreover, it is remarkable that the elimination of the target strain is quite fast, compared to that of the dummy strain shown in Figure 4, which ensures the protection of the target strain.

To sum up, due to the swift shift in equilibrium and the effect of K28 to cause apoptosis, it is difficult to isolate a living target cell from the mixture, and also to obtain the DNA sequence from a dead cell because the DNA is destroyed into fragments by K28. Therefore, it is impossible to steal the sequence of the target cells when a wrong passcode is entered.

Conclusion

In this modeling page, Saccharomyces cerevisiae was used as a model organism for the dummy strain to confirm the behavior of the Dummy System in silico. Initially, the dummy and target strains coexist in a tetracycline medium, and when exposed to light in the correct order, the dummy is deprived of its tetracycline resistance, leaving only the target. When the dummy strain is exposed to light in the wrong order, the target strain undergoes apoptosis due to a killer factor called K28 released by the dummy strain, and only the dummy strain remains.

We first checked the coexistence state of the model organisms. We obtained the parameters of the differential equations by fitting the experimental data of the beginning of the coexistence to a curve that follows Lotka-Volterra's competition equation. Then, based on these parameters, we simulated the whole process of microorganisms reaching the equilibrium state. By this simulation, we confirmed that equilibrium is achieved from the parameters we obtained from existing experimental data. One point that needs to be improved would be that the population ratio in the equilibrium state is not so overwhelming, since the dummy is on a scale of 1.53 times the target.

The simulation with the correct input showed that most of the dummies will die in about two days after the input. This result can be judged to be on a practical scale. One of the reasons why it is on a practical scale is that the same scale of time is required when the yeast is cultured by itself. Therefore, the time required doesn’t change whether you use strains in the Dummy System or non-consolidated strains. Maybe you would think that we could do more to shorten the shift span, but we assume that this is enough. In the case of a wrong light order input, the period until the shift should be shortened for security reasons, but when the light input is correct, there is no need to do so since the urgency of the shift in this case is not directly related to security. In the event of a wrong input, the mixed group must be quickly transformed into a dummy strain-only group to prevent the target strain from leaking out by means of device destruction, for example. In fact, in this simulation, most of the target strains are eliminated within 10 hours, ensuring the protection of the target strain.

From the above, numerical simulations confirmed that this Dummy System has a mechanism that can withstand practical use.

Future Prospects

In this modeling, we tested the conditions under which the two species coexisted before decoding. However, in general, it is not obvious or easy to assume two similar organisms to coexist in the same environment for a long period of time from the viewpoint of the competition exclusion rule, an ecological theory. Therefore, if we can create a general framework in which the Dummy System can be applied to two species for which it is difficult to find or establish such coexistence conditions under experimental conditions, it will be a more versatile technology.

One way to realize this framework is to keep two species coexist for as long as possible, if it is difficult for them to reach an equilibrium state in the end. If both of them manage to survive until they reach the receiver, who is supposed to unlock the system, the mechanism of coexistence and the shift of populations will still work. Therefore, the hurdle here is how to lengthen the coexisting time and what to do if the dummy dies out before reaching the receiver, which means the loss of the security function.

About the first hurdle, since the time until one of the two species is eliminated from the medium can be extended by adjusting the initial population ratio between the two species, it should be possible through experimental verification to extend this time limit somewhat to suit one's purposes.

About the second hurdle, the solution we have in mind at this point is to set a "time limit for decode" on the hardware side. For example, let us consider a situation in which dummies fall below the target population in 360 hours without any operation from the start of incubation, and then eventually die out. In this case, if the "time limit" is set to less than two weeks, shorter than 360 hours, and natural light is irradiated to the culture medium from the hardware when the expiration date arrives to kill the target, the two species can be treated as coexisting even if it is not in an equilibrium state, until the expiration date. Since the final equilibrium state resulting from the operations performed before the expiration date does not depend on the concentration of each strain and the timing of the cipher input, it does not matter whether the two species coexist in equilibrium or nonequilibrium before the cipher input, and in any case, the shift by the cipher input works correctly. Using these ideas, the Dummy System will be able to maintain its security function and lengthen the term of validity of itself at the same time.

In this way, a new and more versatile Dummy System can be obtained.

Modeling in Genochemy - Building Equations from Network of Actors

Purpose

Genochemy, a block programming tool for synthetic biology, has a full-scale differential equation simulation implemented in the background, as we wanted to pursue realism while making it an educational tool.

In the generation part of the differential equation (how to convert the sequence of connected blocks into a simulation of mRNA/protein concentration), a network-like structure, which can be called a "model," was used, and is described here. The model is intended to efficiently and ideally generate differential equations.

Figure 1

Figure 1. The diagram of the problem this model deals with.

Method and Modeling

It is easy to list the available mRNAs and proteins from the connections of the blocks on Genochemy. They can be enumerated by examining the position of promoters, terminators, and each block one by one. However, it is not as easy as expected to automatically calculate differential equations from mRNAs and proteins, ideally and efficiently.

Figure 2

Figure 2. The detailed diagram of the problem this model deals with.

There are several problems, but most arise because of the difficulty in modeling the situation where matter (protein or mRNA) interacts with each other. For example, Red-light-sensor a and Red-light-sensor b form a heterodimer when exposed to red light. Therefore, we have to consider the monomers a and b and the dimer ab as a system of differential equations. We can decide that the output of the differential equation for the dimer ab should be taken by the monomer a protein, but this is not beautiful and will be difficult to manage as Genochemy is updated and more and more blocks are added. We have a new idea for a "two-protein interaction". We thought that it would be most ideal if the "interaction of the two proteins" itself would be responsible for the output of the dimer ab differential equation.

Figure 3

Figure 3. The example problem, which could be solved by dealing with the interaction between two proteins.

Other issues that needed to be considered included what matter (protein or mRNA) should take the role of outputting the term that determines the amount of transcription and the degradation in the differential equation, which is another issue arising from the question of how to model the relationship between matters. Therefore, we modeled "actions from things to things" and "interactions between things" so that they could be expressed in the program. Specifically, in addition to mRNA and protein, promoters and degraders (assumed to be responsible for degradation), which are important when considering corresponding terms in differential equations, were all considered as actors, and a network of actors was first constructed from the connections of blocks. Then, edges (and some actors) in the network output terms of differential equations for each of them. (The name "actor" is named after Actor Network Theory (ANT)).

Figure 4

Figure 4. From blocks to the network of actors.

Figure 5

Figure 5. From the network to the system of differential equations.

Example

In the circuits shown in Figure 4 and Figure 5, the part of final system of differential equations looks like this (The equations regarding GFP is not included):

d[mRNASaSb]dt=k1k2[mRNASaSb]d[Sa]dt=k3[mRNASaSb]k4[Sa][Sb]+k5[Sab]k6[Sa]d[Sb]dt=k7[mRNASaSb]k4[Sa][Sb]+k5[Sab]k8[Sa]d[Sab]dt=k4[Sa][Sb]k5[Sab]k10[Sab] \frac{d[\mathrm{mRNA_{SaSb}}]}{dt}=k_1-k_2[\mathrm{mRNA_{SaSb}}]\\ [0.5em] \frac{d[\mathrm{Sa}]}{dt}=k_3[\mathrm{mRNA_{SaSb}}]-k_4[\mathrm{Sa}][\mathrm{Sb}]+k_5[\mathrm{Sab}]-k_6[\mathrm{Sa}]\\ [0.5em] \frac{d[\mathrm{Sb}]}{dt}=k_7[\mathrm{mRNA_{SaSb}}]-k_4[\mathrm{Sa}][\mathrm{Sb}]+k_5[\mathrm{Sab}]-k_8[\mathrm{Sa}]\\ [0.5em] \frac{d[\mathrm{Sab}]}{dt}=k_4[\mathrm{Sa}][\mathrm{Sb}]-k_5[\mathrm{Sab}]-k_{10}[\mathrm{Sab}]

The term k1k_1 of the equation d[mRNASaSb]dt=\frac{d[\mathrm{mRNA_{SaSb}}]}{dt}=\cdots is generated from the edge from Constitutive Promoter to mRNA-Sensor a, b. The term k2[mRNASaSb]-k_2[\mathrm{mRNA_{SaSb}}] of the same equation is generated from the edge from Degrader to mRNA-Sensor a, b.

The term k3[mRNASaSb]k_3[\mathrm{mRNA_{SaSb}}] of the equation d[Sa]dt=\frac{d[\mathrm{Sa}]}{dt}=\cdots is generated from the edge from mRNA-Sensor a, b to Sa. the term k4[Sa][Sb]-k_4[\mathrm{Sa}][\mathrm{Sb}] and k5[Sab]k_5[\mathrm{Sab}] are generated from the interaction between Sensor a and b. the degradation term is generated from the edge from Degrader. d[Sb]dt=\frac{d[\mathrm{Sb}]}{dt}=\cdots is almost the same.

k4[Sa][Sb]k_4[\mathrm{Sa}][\mathrm{Sb}] and k5[Sab]-k_5[\mathrm{Sab}] are generated from the interaction between Sensor a and b. d[Sab]dt=\frac{d[\mathrm{Sab}]}{dt}=\cdots The decomposition term is difficult to capture, but it was implemented in such a way that it can be thought of as "the interaction created a new actor named Sab, and an edge from Degrader to actor Sab generated the decomposition term."

The terms generated by each edge

We'll show the terms generated by each edge of the network shown in Figure 5. The parameters are assumed in the way described in the below section.

The edge from Const. Promoter to mRNA-Sensor a and b adds the term below to the equation of mRNA-Sensor a and b.

11

The edge from Sensor ab Activated Promoter to mRNA-GFP adds the terms below to the equation of mRNA-GFP.

0.05+[Sab]1+[Sab]0.05+\frac{[\mathrm{Sab}]}{1+[\mathrm{Sab}]}

The edges from mRNA-** to proteins like Sensor a, Sensor b, or GFP (denoted by "encodes" in Figure 5) adds the term below to the equation of the corresponding protein.

1[mRNAself]1\cdot[\mathrm{mRNA-self}]

The edges from Degrader to actors (denoted by "degrades" in Figure 5) adds the term below to the equation of the corresponding actor.

1[Actorself]-1\cdot[\mathrm{Actor-self}]

The edge between Sensor a and Sensor b adds the equation of d[Sab]dt=\frac{d[\mathrm{Sab}]}{dt}=\cdots and the term 0.5[Sab]0.5[\mathrm{Sab}] to the equations of Sensor a and Sensor b, 0.5[Sab]-0.5[\mathrm{Sab}] to that of Sensor ab, [Sa][Sb]-[\mathrm{Sa}][\mathrm{Sb}] to those of Sensor a and Sensor b, and [Sa][Sb][\mathrm{Sa}][\mathrm{Sb}] to that of Sensor ab.

Conclusion

It works well.

When the circuit shown in Figure 4 and Figure 5 is executed, the system is correctly implemented as follows: "When red light strikes, Sensor a and b form a heterodimer, thereby promoting transcription downstream of the red light sensor dimer induced promoter, resulting in increased expression of GFP."

The fact that it works well has been demonstrated in numerous opportunities to use Genochemy.

Assumptions made in the Genochemy simulations

Although Genochemy generates differential equations in this way, the output system is simplified with some assumptions. This is because it was difficult to extrapolate all the parameters from real organisms, and also because it was considered satisfactory in appearance even if the parameters were not those of real organisms.

  • Parameters are simple values (e.g., 1 and 0.05).
  • The coefficient of the degradation term is assumed to be constant regardless of the type of protein or mRNA. The actual coefficients of the degradation term vary greatly depending on the presence or absence of degradation tags on proteins, etc., and also vary greatly depending on the modification of mRNA.
  • The coefficients of the transcription strength and translation terms were also assumed to be constant without considering differences in mRNA or protein. This is a point that was pointed out by the general public who actually played Genochemy. They said in reality the speed of transcription can vary depending on the length of the mRNA. The rate of translation also varies greatly depending on the modification of the mRNA.
  • The concentration of drug A was assumed to change instantaneously when the slide bar was adjusted. (In reality, it is impossible to add or subtract drugs instantaneously.)
  • Recombination of recombinase is a stochastic event, but we did not do a full-scale modeling as shown in the Recombination Modeling section above, but based our calculations on a concentration-dependent probability, assuming that recombination occurs with constant probability (if the concentration is the same) at any time. In reality, it is possible that recombination does not occur for a long time, but as an education tool, we want recombination to occur as soon as practicable, so we adjust the probability of recombination to be higher when the concentration remains high for some time.
  • For the kill switch, the specific metabolic pathway leading to cell death was not considered, but rather a concentration-dependent probability was considered and adjusted so that cell death occurs at constant probability (if the concentration is the same) at each time. As with recombinase, the probability of cell death is adjusted so that the probability of cell death increases after a certain period of time, because it was judged that no kill after a long period of time is not appropriate as an educational tool.
References

[1] Benzinger, D., & Khammash, M. (2018). Pulsatile inputs achieve tunable attenuation of gene expression variability and graded multi-gene regulation. Nature communications, 9(1), 1-10.

[2] Shoura, M. J., Vetcher, A. A., Giovan, S. M., Bardai, F., Bharadwaj, A., Kesinger, M. R., & Levene, S. D. (2012). Measurements of DNA-loop formation via Cre-mediated recombination. Nucleic Acids Research, 40(15), 7452–7464. doi: 10.1093/nar/gks430

[3] Ringrose, L., Lounnas, V., Ehrlich, L., Buchholz, F., Wade, R., & Stewart, F. A. (1998). Comparative kinetic analysis of FLP and cre recombinases: mathematical models for DNA binding and recombination. Journal of Molecular Biology, 284(2), 363–384. doi: 10.1006/jmbi.1998.2149

[4] "Team:Hong_Kong_HKUST/Overview - 2017.igem.org." https://2017.igem.org/Team:Hong_Kong_HKUST/Overview (accessed Oct. 9, 2022).

[5] "Team:Edinburgh_UG/Model - 2017.igem.org." https://2017.igem.org/Team:Edinburgh_UG/Model (accessed Oct. 9, 2022).

[6] Ghosh, K., Lau, CK., Gupta, K., & Van Duyne, G.D. (2005). Preferential synapsis of loxP sites drives ordered strand exchange in Cre-loxP site-specific recombination. Nature Chemical Biology, 1, 275–282. doi: 10.1038/nchembio733

[7] Kawano, F., Okazaki, R., Yazawa, M., & Sato, M. (2016). A photoactivatable Cre-loxP recombination system for optogenetic genome engineering. Nature Chemical Biology, 12(12), 1059-1064. doi: 10.1038/nchembio.2205

[8] Kishi, K., Koyama, H., Oka, S., Kato, A., Sato, M., & Fujimori, T. (2021). Repetitive short-pulsed illumination efficiently activates photoactivatable-Cre as continuous illumination in embryonic stem cells and pre‐implantation embryos of transgenic mouse. genesis, 59(12), e23457. doi: 10.1002/dvg.23457

[9] Optosupply. 5mm Pure Green LED. [Data set]. https://akizukidenshi.com/download/ds/optosupply/OSPG5161A-RS.pdf (accessed Oct. 9, 2022)

[10] Lyklema, J. (1995). 5 - Adsorption of Polymers and Polyelectrolytes. Fundamentals of Interface and Colloid Science, 2, 5-1-5-100. doi: 10.1016/S1874-5679(06)80008-5.

[11] Gause, G. F. (1932). Experimental Studies on the Struggle for Existence : I. Mixed population of Two Species of Yeast. Journal of Experimental Biology, 9(4), 389–402. doi: 10.1242/jeb.9.4.389

[12] Avery, S. V., Malkapuram, S., Mateus, C., & Babb, K. S. (2000). Copper/zinc-Superoxide dismutase is required for oxytetracycline resistance of Saccharomyces cerevisiae. Journal of Bacteriology, 182(1), 76–80. doi: 10.1128/JB.182.1.76-80.2000

[13] Tarhan, C., Pekmez, M., Karaer, S., Arda, N., & Sarikaya, A. T. (2007). The effect of superoxide dismutase deficiency on zinc toxicity in Schizosaccharomyces pombe. Journal of Basic Microbiology, 47(6), 506–512. doi: 10.1002/jobm.200700220

[14] Koike, H. [小池文人]. (2004). Seitaigaku nyumon [生態学入門. Introduction to Ecology]. Tokyo, Japan: Tokyokagakudojin [東京化学同人].

[15] Chopra, I., & Roberts, M. (2001). Tetracycline antibiotics: mode of action, applications, molecular biology, and epidemiology of bacterial resistance. Microbiology and Molecular Biology Reviews : MMBR, 65(2), 232–260. doi: 10.1128/MMBR.65.2.232-260.2001

[16] Sheppard, S., & Dikicioglu, D. (2019). Dynamic modelling of the killing mechanism of action by virus-infected yeasts. Journal of the Royal Society Interface, 16, 20190064. doi: 10.1098/rsif.2019.0064