MoDeL Demo

Part1: toggle switch.

Toggle switch is characterized by mutual inhibitory network. That is, the protein product of one gene represses the promoter of the other, and vice versa. So with proper parameters, the system will stay at one stable steady state at equilibrium. In other words, one gene's product will be domanant. But with proper inducers (usually IPTG or ATC), the system manage to flop from one state to another, still keeping steady.

We tried an instance of a toggle switch composed of a promoter which can be repressed by laci protein, and is denoted as pLac(r0010). And the other promoter is tet-repressible and be denoted by pTet(r0040). The plac promoter promotes the expression of TetR(c0040), And the LacR(c0012) is initiated by pTet(r0040).

Without inducers, which protein will dominate depend on the binding affinity of the repressor to the promoter. In this case, it's the LacR protein take control at null input of inducers. The only input to the algorithm is a input file like this.

Here are a bit explanations to the input file

As to the seedspecies section:
# Medium iptg : IPTG is contained in the compartment named medium
# Ecoli dna1 d:r0040(tetr1,tetr2)-b0034(rib)-c0012(dna,iptg,dim)-b0014()
#dna1_init
# d:means it's a strand of dna composed of 4 parts.
# The promoters have two binding site(tetr1 and tetr2) for associative repressors.
# c0012(dna,iptg,dim) means the c0012 is a DNA sequence. Iptg means the product #of c0012, which is LacI protein, in terms of dimer, can be binded to iptg which #deactivate it's repression to placi promoter.
#b0014()is the terminator.
#dna1_init means the initial concentration of dna1 is dana1_init which is given in the parameter section of the input file.

Once the input file is provided, both the net and SBML file will be generated within a second. The species "discovered" by the algorithm is illustrated below, totally 18 substances in such a simple network.

A bit explanations are given to the table above.

1.
#nb: non-biobrick part, indicating a constraint of only one part constituting a sequence, making it significantly distinguished from DNA, RNA, Protein sequences.
#d: means the following thing is a sequence of dna with structure information.
# b0014() is a terminator with no binding site.
# the hyphen sign connects the parts on a dna.
#C0012*(dim,dna,iptg) means the c0012 is a dna sequence, and the product of c0012 which is #laci protein, can form dimers, iptg can bind to monomers of laci .

2.
As you can see from species 12 to 15, the advantage of this approach is that it contains information about the structure of the species. For instance, s12 is nb:i0001(laci!1).p:c0012(dim!2,dna,iptg!1).p:c0012(dim!2,dna,iptg)

The ! is used to denote binding , here iptg bind to the laci denoted by !1, and the . is used to separate different molecules of one species. The molecule that comprise the s12 species are nb:i0001(laci!1)，p:c0012(dim!2,dna,iptg!1)，p:c0012(dim!2,dna,iptg).

The two proteins bind to form dimers through bingding !2.
So now you can imagine the structure of s12 species.
The network generated by our software is as below. User can drag the nodes to the place they want to have a better clarity of the network.

One thing to note is the assumptions made by Chen LIAO:
1.Both LacR protein and TetR protein can form dimers. But only the dimers, rather than the monomers, can bind to the promoter regions to repress the expressions of the downstream genes.
(#note that this a an advantage of our algorithm, it is more close to biology reality)
2.IPTG molecules can bind to LacR proteins or each monomer of LacR Dimers but cannot bind to LacR dimers whose DNA-binding domains Have been occupied.
(#When there are less LacR proteins, there are less LacR dimers. When there are less LacR dimers, the complex of LacR dimer and DNA is more likely to disassociate.)
3.There are leaky expressions for repressed pLac and pTet promoters.

In order to test the ability of the system to switch from one state to the other, we conducted a time course simulation. The IPTG are added to the system at time 10000s.
The GFP expression level changes triggered by the pulse are illustrated by curves with different colors.

As you can see from the figure above, at initial time, the LacI protein soon takes control and reaches equilibrium. When IPTG is added to the system, there is a sharp decrease. Since the half life of LacR is short and IPTG is so favorable to bind with LacR, making it inactive (lose the ability to bind tightly to the promoter that helps to express TetR)

Part2: RTC two counter (RiboRegulated transcription cascade)

Recent years saw a emergence of the small RNA application in the gene regulation network. This is a timely way to regulate gene expression compared to other transcription level modulation.

As a mimic to the electronic digital circuit that can count the pulses or other events (a basic function of most MCU module), synthetic biologist designed ways to count biological events such as the adding of inducers. Here, we adopted the RTC-two counter as an example.

The constitutive promoter pl tet 0-1 , which corresponds to j23100 part, drives the T7 RNA polymerase T7RNAP (i2032 part), the RNAP binds to THE T7 promoter PT7 and the GFP will be expressed then.

Between them is the part j01010, which consist of the rbs sequence and the complementary sequence cr ( so the cr section of the mRna transcribed will form a loop with the Rbs section ,thus inhibiting translation)

B0014 is the T7 promoter. The RNAP binds to the T7 promoter to drive the expression of the downstream gene (GFP). But still, the translation is inhibited by the stem loop formed by the cr and RBS section of j01010.

To release the inhibition, the arabinose induced promoter PBAD can drive the expression of J01008, which corresponds to taRNA. But this promoter can be repressed by c0080, which is the AraC protein dimer. When arabinoses are added to the system, the arabinoses bind to the AraC monomers and prevent them from forming dimers. Thus release the repression to the i01008 promoter.

The taRNA is specific to the cr sequence and it bind to the cr sequence, thus open the stem loop, allowing the expression of RNAP and GFP.

The input file to this case is as below:

The network view generated by our software is as below (this visualization is realized by Junyuan Xie, Fangming Liu, and Chuocheng He. You can drag the nodes to your preferred place.

The list of species generated by the algorithm is listed the table below.

Below is some assumptions made by Liao Chen who worked out this example with his algorithm.
1.AraC proteins cannot degrade so as to keep forever repression of pBAD before any pulse of arabinose.
2.AraC protein can form dimers. But only the dimers, rather than the monomers,can bind to pBAD and repress the expressions of the downstream genes.
3.Arabinose molecules can bind to AraC proteins or each monomer of AraC dimers but cannot bind to AraC dimers whose DNA-binding Domains have been occupied.
4.No leaky expressions is possible for both pBAD and T7 promoter, Which can reduce the leaky GFP expression (noise) before the second.

Pulse comes and make the leap of GFP amount more impressive and distinguishable.

We also attest the validity of our approach using the time course simulation in this case.

As you can see from above, two Arabinose pulses are added to the system at time 2500.01s and 2561.01s. The black curve shows the normalized GFP expression of the cells that received both the two pulses. The red line refers to the cells that only receive the first arabinose pulse and the blue line refers to the cells that only receive the second one.

It is clear from the figure, that at the point of the pulse, the GFP level begin to leap to a higher state. Since the inducer is removed from the environment after it generate the pulse, the expression of RNAP will not last long so there is only a small leap if the cell received only one pulse, if with the second, the GFP level will leap markedly as indicated by the blackline.

Extended work and future plan

We tried on the following prospect and finished the test code with Matlab (ready and eager to implement the test codes into Lachesis)

• Automatic parameterfitting, adjustment, and estimation with a given mathematical structure. User selects or describes the standard or criteria according to wet LAB data or just describe. Minimize the cost function using simulated annealing.

• Assess a design with user input and pre-defined criteria from different aspect: Analyze the sensitivity and robustness of a given network.

1 Parameter section:

1.1 Background

Almost ideal tools that already exist in electronics: User provide the part’s name, the software fetch the standard model and its associated parameter, and return a complete mathematical model with all the parameters known (for instance, the PROTEL-software can easily get the released device from the NI company with all the response parameters standardized from a device library).

But temporally, the knowledge and cognition level of synthetic biology do not support this: Most parts do not have a standard model associated with it. Some parameters of the parts haven’t been measured yet. For some parts with quantitativeparameter value, which is highly context dependent, is hard to transfer with defined parameter across different hosts.

1.2 More considerationM

Why do we try automatic parameter fitting(or estimation, adjustment)？

• If it’s a big network with too many parameters undefined, it would get the user exhausted adjusting all the empty parameters by hand according to his/her experience and web resource. Networks generated by rule based modeling is ordinarily large

• Compared to manually parameter adjustment, in silicon auto parameter adjustment is more convenient.It givesan estimation of the interest parameter value according to wet lab data, thus free the user from spending too much time on estimating a good parameter.

1.3

We adopted a demo from a Tinkercell tutorial website on implementing a simplerepressilator. After the network is generated by hand, the parameters are left behind to the user to adjust themselves. But as you can see from the following four figures, the process is tough and time consuming. (more snapshots of the process are eliminated )

zwx

………………（more trial and error）

The four figures above show a difficult parameter adjustment process by hand.

By contrast, parametersestimated by our approach in silicon can be done in a timely fashion. It’s fast and convenient(see figure below)

1.4 Algorithm

For PSO and SA algorithm detail, please go to our algorithm and technology wiki part.

The user input behavior is compared with the result generated by solving odes using cross correlation or absolute distance or fourier transform (user can select the criteria)

Part 2 Assessment section

2.1

Why do we try to evaluate a design with sensitivity and robustness analysis?

• Computer descriptions of biological mechanisms are by necessity,very simplified.

• The biology system is so complex that current mathematical description is far from accurate

• Parameters of biological parts are context dependent and some may vary drastically when host cell or environment change

• Noise come from various sources and can drastically influence the performance of a circuit

• Things that can be implemented through ode parameter adjustment might not necessarily be observed in real experiment.

• If the parameter range is slim at the working region, then most probably the desired phenomena will never be observed in wet experiment

• If we get the most sensitive parameter, we probably know what to do to refine a circuit to quick tune it to a desired behavior or just try to keep that most sensitive parameter as constant as possible.

Part3 Model reduction and network modulation (havn’t done yet)

Below is an aflow chart that demonstrates the whole process: