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 ($\mathrm{TF}_\mathrm{on}$) and off ($\mathrm{TF}_\mathrm{off}$), and that EL222 in the state of $\mathrm{TF}_\mathrm{off}$ will transition to $\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 $I$, it can be written as $I\cdot k_\mathrm{on}$, where $k_\mathrm{on}$ is a constant. Also, $\mathrm{TF}_\mathrm{on}$ changes to $\mathrm{TF}_\mathrm{off}$ at a constant rate $k_\mathrm{off}$ regardless of the light intensity. Putting these together, we can write the following differential equations.

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

the differential equations can be rewritten as follows.

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 $\mathrm{TF}_\mathrm{on}$ as a variable and translated at a constant rate of $k_\mathrm{trans}$. Also, assuming that mRNA is transcribed at a constant rate of $k_\mathrm{basal}$ even if TFs in the on state is not present, and mRNA and protein are degraded at a rate of $k_\mathrm{degR}$ and $k_\mathrm{degP}$, respectively, the following differential equation can be written.

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

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").

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=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 $73 \mathrm{nM}$, which is about ($\frac{1}{300}$) of that in the presence of light.

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.

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.

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].

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 $S_0$ and $S_1$ in Figure 2 are as follows:

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

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

It is needed to adjust the $\Delta t$ so that $p_i$ satisfies $0 \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]$ is constant as long as there are no external influences.

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 $I$ transitions to state $P$ where two DNA sequences ($P_a$ and $P_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.

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: $p_1$, $p_{-1}$, $p_2$, $p_{-2}$, $p_3$, $p_{-3}$, $p_4$, and $p_{-4}$, described in Figure 3 below.

The state after recombination is now denoted by $P$.

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

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].

$k_3$ represents the value at 3044 bp, and $k_{-3}$ and $k_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 $P$). As time passes, the recombination probability increases. In other words, the "probability of being the post-recombination DNA sequence" increases.

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 $S_3$ and $S_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.

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

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.

(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 $p_1$ and $p_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.

$p_1$ and $p_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 $p'_1$ and $p'_2$ can be defined by

$p_{-1},p_{-2}$ are assumed to be unchanged here. Since $p_1 \gg p_{-1}$ and $p_2 \gg p_{-2}$ in the situation we are considering here, we assume that the changes in $p_{-1}$ and $p_{-2}$ are negligible compared to those of $p_1$ and $p_2$.

We considered that $p_3$ and $p_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.

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].

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 [\mathrm{\mu W/cm^2}]$), indicating that recombination can occur sufficiently with PA-Cre as well. (Figure 9)

Then, we simulated the case when about 1% intensity of the blue light ($I=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.

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 $p_3$, $p_{-3}$ and $p_4$ was large. This is confirmed by the fact that the stochastic model in Figure 5 often stays in state $S_4$ and $I$, 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 $k_3$, $k_{-3}$, and $k_4$ (note: $k_{-4}$ in the paper's orientation), corresponding to $p_3$, $p_{-3}$, and $p_4$, depend on the sequence length between recognition sequences.

$k_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 $k_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 $k_{-3}$ and $k_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 $k_{-3}$ and $k_4$ at each sequence length in the paper.

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

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 $p_3$ transition (from $S_4$ to $I$), 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 $p_3$ of the transition from the state $S_4$ to the state $I$ can be written as $p_3=\beta\cdot q(N)\Delta t$ using the probability $q(N)$ of the recombinases attached to the two recognition sequences approaching within a certain distance and the constant $\beta$. $N$ is the sequence length between the recognition sequences.

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

Therefore,

Fitting to the $k_3$ literature values for the two $N$, we obtained $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).

We first simulated in Figure 6 using the parameter $k_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.

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.

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."

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]:

where $N$ stands for population, $t$ for time, $r$ for proliferation rate, and $K$ 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 $N_1$ as the population of microorganism of the dummy strain and $N_2$ as the population of that of the target strain, the Competitive Lotka-Volterra’s equations are obtained as the following simultaneous differential equations:

where $r_1$, and $r_2$ stand for proliferation rates, $K_1$ and $K_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.

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.

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.

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.

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:

This can be rewritten as follows by arranging the coefficients:

Figure 3 shows that microorganisms are in an equilibrium state after $t=400$, so $t=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 $r'_1=0$ since it cannot multiply any more. For simplicity, the following assumption was made.

The setting of $r'_1=0$ at $t=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 6 shows a linear increase after $t=400$. Two days after the irradiation, the amount of target is about $e^{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.

Parameters used before, such as $r_1$, $r_2$, $K_1$, $K_1$, $\alpha_{12}$, and $\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 $N_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, $r_\mathrm{kill}=0 (t < 400), 0.64 \mathrm{[/h]} (t\geq 400)$.

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.

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.

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.

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)).

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):

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

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

$k_4[\mathrm{Sa}][\mathrm{Sb}]$ and $-k_5[\mathrm{Sab}]$ are generated from the interaction between Sensor a and b. $\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.

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

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.

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

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