In GRN, there are two kinds of files: TF to TF and TF to Gene. Since the database about the regulation between TFs and Genes contains only one-way interaction, the matrix of GRN is a rectangle.

First of all, read the regulation relationship of TFs. Our software filters the documentation of RegulonDB on the head of all files and then reads the name of regulate and regulated TF, which is also the name of its genes, one by one. In the same time, our software numerates the genes and stores their names into an objects’ array of genetic data.

  The format of regulation database:
    TF_name    TF_name    +/-/+-

The regulation of TFs has been put into a square matrix whose row is the regulator and column is the one regulated by. To make our GRN as complete as possible, the regulation between TF and genes has joined into the matrix. The one-way interaction results that we must read the TF in order to fulfill the regulator before completing the TF to gene’s regulation in the same way of TF to TF.

  The format of regulation database:
    TF_name    Gene_name    +/-/+-

At last, a regulatory matrix whose row represents regulate gene (TF) and whose column represents gene regulated by (TF+Gene) has been output into a file called “old_GRN” in root directory. The values in GRN matrix are regulations in which “1” means positive activation, “-1” means repression and “0” means no relationship. There have been some regulations both positive and negative identified regulations are determined by the experimental environment. As a result, our software picks out those uncertain genes and stores them into a file named “uncertain_database”.

  The format of uncertain database:
    ?    Gene_name->Gene_name

The question mark represents the unknown regulation between regulator and regulated-by whose names presented afterward. Users could replace the question mark with the data known in past experiment. (“+” rep positive, “-” rep negative). Our software will replace the values in matrix automatically. But if not rewrote, our software will regard those regulation as unknown.

All gene information has been deposited into a file named gene_info which could be downloaded here. In order of picking out the genes in GRN as fast as possible, all genetic information are stored in a “map”. “Map” is just like a dictionary yet its words are names of genes and its descriptions of words are replaced by genetic information. By using binary tree method, it is very fast to search the “word” wanted in the “dictionary”. As tested, the speed of binary tree method built-in “map” function is 720 times faster than traversal method.

  The format of Gene Info database:
    ID_assigned_by_RegulonDB    Gene_name    Left_end_position    Right_end_position    DNA_strand    Product_type    Product_name    Start_codon_sequence    Stop_codon_sequence    Gene_sequence

The label of the map vector is gene name which will be picked out based on the names read in regulation matrix before. It is really fast using the binary tree method to find the specific genetic information and store them into a specific object. Those information includes gene ID, left position, right position, gene description and gene sequence. The gene ID is used to link to RegulonDB’s gene details; The left position is used to find its specific transcription unit; The right position is used to figure out the base amount; The description of genes is used to distinguish the RNA and protein; The sequence is used to predict the regulation by alignment.

All promoter information has been deposited into a file named promoter_info which could be downloaded here. But we also need transcription unit information because the information files about promoter do not contain all genes’ names backward. “TU Info” file, which can be downloaded here, contains the starting position of each TU and its promoter name. Our software picks out the starting position into a integer array. Using the left position picked out in gene info, our software would find out which unit the gene belongs to through dichotomy method and then stores the name of promoter into corresponding object.

  The format of TU info database:
    Operon_name    Unit_name    promoter_name    Transcription_start_site ......

The principle of fetching information of promoters is same as fetching genes’s. Our software stores the promoter information from the file named “promoter_info” in a “map” which could be used to pick out the promoter sequence by searching promoter name through binary tree method.

  The format of Promoter Info database:
    Promoter_ID_assigned_by_RegulonDB    Promoter_name

The sequence of promoter will be used in the alignment method in next module which could make a prediction of exogenous genes’ regulation pattern.

Our software integrates all information we picked out about genes and generates a file named “all_info” —— all information about genes —— for the output graphical interface’s reading. In the meanwhile, the array of objects containing all information has been stored in computer memory which greatly improve the computing speed of our software.

  The format of all_info database:
    No.    promoter_sequence    gene_sequence    gene_name    ID    left_position    right_position    promoter_name    description
The fetching module generates three files: old_GRN, all_info and uncertain_database.

We would like to start with a simple example. A cell operates like a basketball team. Imagine you are a manager of a team who wants to bring in some talent players making up a “big three” and build a champion-potential team this season. Before you pay the sky-high bills for the “big three”, you can evaluate the effect of the talent introduction. New members‘ records are good reference, but not the whole thing.

There are various factors influencing the effect of introduction. Let’s carefully choose one of the most profound aspects and focus on the relationship of the members. In the original team, you are familiar with all players’ characteristics, their roles in the team and the coach’s style, i.e. you have the information of the original player interaction network.

You know Alex is a good shooter and Bob is a strong centre forward. Carl is your target player. Carl is famous for his shooting skills and appears dominant in the court. In other words, Carl shows more obvious similarity with Alex but a low level of similarity with Bob. Then, in the new player interaction network, Carl might play a role 80% like Alex and 20% like Bob. He is similar with Alex and Bob, however different. That’s an analysis and prediction at the global point of view.

[ Pic.1 Alex, Bob and new member Carl]

Just like the basketball team example, researchers often need to insert exogenous genes(new players) into a cell(original team) to achieve a specific goal(win the champion). In the past, the behaviors of exogenous genes are mainly speculated by wet lab experiments. Now we are trying to give an answer before wearing laboratory gloves.

Regulatory Model
Regulation of gene expression includes 4 levels:
•Level of DNA rearrangement.
•Level of transcriptional regulation.
•Level of translation.
•Level of post-translation
This year we focus on the level of transcriptional regulation both for the importance of the level and model simplification. By carefully examining the lac operon system, which is widely considered as the first discovery of the gene regulation system, we constructed our regulation model with functional units called “Regulation Unit” [FXIME: regulation or regulatory?]. A regulation unit consists of two segments. The first one is a promoter sequence and the second one is a protein coding sequence.
[Pic. 3 Promoter Sequence, Protein Coding Sequence and Regulation Unit]
In our model, the promoter segment is regarded as the main regulated region. A transcription factor is a protein that binds to specific DNA sequences, thereby controlling the flow of genetic information from DNA to mRNA. The binding sites are promoter regions of DNA adjacent to the genes that they regulate. At first, according to the lac operon system, regulation units as a regulatory target with the same promoter are supposed to have same behavior. But we found it is insufficient because there are units with the same promoter showing different properties. Then we took the genes regulated by transcription factors into consideration. The different properties of two units are first owing to their promoters. If they have the same promoter, their protein coding sequences are supposed to make the difference. By taking this method, it turns out that this model works better.
[Pic. 4 Regulated Region]
A unit in the network regulates another through the transcription factor. That is, the product of the protein coding sequence of the unit is a transcription factor and the transcription factor regulates the promoter of the another unit.

The basic idea behind the prediction model is deceptively simple: the more similar two sequences are, the more likely they have similar behaviors. In fact, it is extremely difficult to predict an exogenous gene’s behavior because of the complexity of the problem, random noise of the system and the coupling of biosystems.
Advanced alignment algorithm is selected to reduce the complexity. Sequences which contain all the information of the species are the entity of the gene regulatory network. Sequence similarity is an essential concept to the prediction model. The selected alignment algorithm can significantly reduce the complexity of the problem and makes it possible to give a reliable prediction from a global point of view.
We designed a random method to filter the noise in sequence alignment. There are no totally different sequences. Even the similarity of any two random sequences is not zero. Filtered results are more significant and reliable to the following steps.
Coupling of biosystem is also simulated at some level. When predicting exogenous gene’s behavior, all the units in the original gene regulatory network are taken into consideration.
Given that the exogenous gene may have never been inserted into E. coli before, all possible reactions in gene regulatory network are reserved to be filtered.
Using the innovated methods above, we are trying to challenge the difficulties and obtain a global perspective of the relationship between the exogenous gene and the original gene regulatory network.

The GRN matrix is the mathematical description of gene regulatory network in which “1” represents “enhance”, “-1” represents “repress” and “0” represents “no regulatory relationship”. The units(RU) in x-axis regulate the units in y-axis. A row can be seen as a vector containing all the information of the target(corresponding unit in the y-axis). Similarly, a column can be seen as a vector containing all the information of the regulator(corresponding unit in the x-axis).

The sequence similarity is obtained by sequence alignment based on Needleman-Wunsch algorithm. The Needleman-Wunsch algorithm performs a global alignment on two protein sequences or nucleotide sequences. It was the first application of dynamic programming to biological sequence comparison.
When dynamic programming is applicable, the method takes far less time than naive methods. Using a naive method, many of the subproblems are generated and sovled many times. The dynamic programming approach seeks to solve each subproblem only once. Once the solution to a given subproblem has been computed, it is stored to be looked up next time.
[Pic. 5 Dynamic programming and naive method]
Like the Needleman-Wunsch algorithm, of which it is a variation, Smith-Waterman is also a dynamic programming algorithm. But it is a local sequence alignment algorithm. The famous BLAST(Basic Local Alignment Search Tool) is improved from Smith-Waterman algorithm. Although local algorithm has the desirable property that it is guaranteed to find the optimal local alignment, we decided to choose the global one because we regarded the segment sequence as a unit.
Sequences are aligned with different detailed methods in different situations. In the regulated side, what we care about is the DNA sequence. In the regulating side, it is the amino acid sequence. When it comes to predict the regulated behavior, we use a DNA substitution matrix to align promoter and protein coding sequences. In the prediction of regulating behavior, the substitution matrix BLOSUM_50 is used to align the amino acid sequences translated from protein coding sequences.
The promoter similarities of the query unit and subject units are stored in a vector. The protein coding similarities are stored in another vector. These vectors are prepared to be used in the new network construction.

Once the similarity vectors are calculated, the next step is to filter them. As explained in the previous part, there is random noise in sequence alignment. In order to filter these meaningless values, a certain amount of random sequences are generated for each query-subject alignment. Normally, 100 is sufficient. Because the sequence length will influence alignment result, random sequences are fixed at the same length as the query one. Then align these random sequences with the subject sequence. The statistic result of these random similarities will be used as a threshold. If the original similarity is lower than the threshold, it is abandoned. In this case, the original value is usually short of statistical significance.

If there are units in the original network having the same promoter as the exogenous one, the first step is to pick them out. Positive and negative regulations of these units are counted separately and distinguished into “positive group” and “negative group”. Then compare the exogenous one with these units. The similarities have already been calculated and stored in the corresponding positions in the similarity vector. The similarity mentioned here is the similarity of protein coding sequences as explained in the model part. The next step is to calculate the average similarity of each group. The exogenous unit is supposed to have the larger one’s direction(positive or negative). The weighted average regulation value of the chosen group whose weight is the sequence similarity is the new element’s value. It means regulatory intensity.
If there is no unit having the same promoter as the exogenous one, given that the promoter is the main regulatory region, the promoter similarity is used as the weight. And the weighted average of the regulation of the whole column is the new element’s value

The construction of the new regulating vector is achieved in a way similar to the one described above. By calculated the weighted average regulation of a row, the program gives the regulatory intensity that the exogenous unit regulate the corresponding unit in the network.

The behavior similarity of two units can be described by the dot product of two regulated vectors or two regulating vectors. A more intuitive way is using the vectorial angle to measured the similarity of two behaviors. But there are some zero vectors in the gene regulatory network which usually means the units either play the role of target or the regulator.
[Pic. 4 GRN matrix, target vector, regulator vector and their dot product]
We have tested the hypothesis by analyzing all 1748 regulation units of Escherichia coli, K-12, recorded in RegulonDB. By pairwise comparison of all these units, about 1.6 million sets of data was obtained. Each set of data consists of promoter sequence similarity, protein coding sequence similarity and behavior similarity. We hope to find some structure in the data that supports our hypothesis. And it is lucky enough to find there is a tendency showing the relationship between sequence similarity and behavior similarity(Pic. 2).
[Pic. 2 Sequence similarity and behavior similarity]
Sequence similarity is set as x axis and behavior similarity is set as y axis. Obviously sequence similarity is continuous-valued (from 0 to 1) and behavior similarity is discrete-valued. Values of behavior similarity determined by the dimension(N) of the vector are between -N and N. According to the result, promoter sequence similarity mainly distributes from 0.4 to 0.6, protein coding sequence similarity mainly distributes from 0 to 0.7 and behavior similarity mainly distributes from -3 to 5. As it is shown in Picture 4, high behavior similarity is partial to high sequence similarity. Peak value of behavior similarity, 17, appears where sequence similarity is 0.537. When behavior similarity value is fixed, for example, set behavior similarity as 8, it is obvious that the higher the sequence similarity is, the more intensive the dots are.

Network analysis includes finding stable condition of network, adding new gene, finding new stable condition and changes from original condition to new condition. We use densities of materials to describe network condition. If all material densities are time-invariant, we can say the network condition is stable.

Regulation relationship in genetic network includes positive regulation, negative regulation, positive-or-negative regulation and no regulation. We store regulation relationship in matrix R. Rji means the unit in line j and row i. For the material of original network, Rji=1 means material i enhance material j, Rji=-1 means material i repress material j, Rji=0 means material i has no influence on material j, Rji=2 means material i enhance or repress material j. For the new material, Rji ranges from -1 to 1. Rji<0 means the possibility of positive regulation is Rji; Rji>0 means the possibility of negative regulation is –Rji; Rji=0 means there is no regulation from i to j. We use Hill equations to describe intensity of regulation. Equations are like following:



The left side of the equation is the derivative x(density) on t(time).”qi”,”pi”,”ri”,”mi”,”ni” are parameters, which determine the intensity of regulation."ri" is degradation rate. Mji is exponent. M is a matrix whose dimensions are equivalent to R's. Mji is 0 or ranges from 0.5 to 1.2 or ranges from -1.2 to 0.5. For the material of original network, if Rji=1,Mji ranges from 0.5 to 1.2;if Rji=-1, Mji ranges from -1.2 to -0.5; if Rji=2;Mji ranges from -1.2 to 0.5 or 0.5 to 1. These Mjis’ absolute values are given randomly by program. If Rji=0, Mji=0.
For the new material,

Stable condition is the condition in which densities are time-invariant. We store material densities in a vector and solve the differential equations with Euler’s formula, which is like below



We know the network will be stable at last, so every material density has a limitation.

Record the original stable condition, set new material density to 0 and this is the new initial density vector. Solve new equations and record density vectors before the new condition is stable and store these data in a text file.

To evaluate the new network, we introduce the grading system.



"xi" and "Xi" are densities of material i, which is not the new material."ny" is the number of materials. The more new densities are close to the original, the less the influence the cell endues. In general, cells close to the original cell are more likely to survive than those who are far different from the original cell. That is the thought of the grading system.

We did a lot of running and found that the “AbsValue” ranges from 0 to 370, so "ScoreA" ranges from 0 to 4.9.We get the integer part and store it in an array, which has five sections. Generate 100 or 200 matrix M from matrix R and run the original and new network for each M, so we can get 100 or 200 of "ScoreA"s. The section which has maximum "ScoreA"s is the eventual score.

In some cases, importing exogenous gene is for enhancing or suppressing the expression of some specific genes in engineered bacteria itself. But it is hard to choose an appropriate regulatory gene. Our software analyzes the GRN forward as well as simulates by optimization algorithm backward for giving a reference of choosing to the users. Our software not only focused on the direct regulation but also focused on the global GRN. In the same time, controlling the expression of multiple genes in network has been realized by global prediction. What’s more, Particle Swarm Optimization (PSO) Algorithm makes it possible.

Before prediction, the expression of specific genes which the experimenter needs should be input into our software as well as the improvement or depression. The number of target gene is SEVEN at most.

It is a must that figuring out the strongest and weakest expression strength before inputting the extreme cases into the target expression. The way to find out the strongest and weakest expression is modeling the GRN’s steady state by a large amount of random regulation from -1 and 1. On the other hand, the expression of genes unpicked by the users should be stable as much as possible. The initial strength of expression is calculated by modeling the original GRN with Hill’s equation.

For getting the best regulation, our software uses PSO algorithm based on 30 particles to simulate the GRN’s changing. First of all, the interactions of regulator and regulated-by have been put into those particles in random so that each particle will have the whole set of regulation. Secondly, the variance between target expressions and stable expression of new GRN have been regarded as the optimize requirements in PSO algorithm. As a result, the minimal variance of 30 particles is the global optimum and the minimal variance of the procession in one particle is the local optimum. Then, taking a step towards global and local optimum as well as considering the inertia and perturbation avoids falling into the sub-optimal condition.

At last, when the variance of expression reaches an acceptable range, our software picks out and saves the best global optimum particle following by the movement of those particles stop.

We constantly revises the factors in PSO algorithm by machine learning method for accurate simulation with a fast PSO particle-motion equation. At the same time, our software also filter the result of regulatory value which is more intuitive.

To improve the efficiency of choosing a suitable gene after getting a series of regulatory value, our software picks out some obvious regulation. The value of regulation is between -1 to 1 in which -1 means negative effect and 1 means positive effect. As a result, what our software has done is filtering out the absolute value which is lower than 0.9. Because the difference of regulatory intensity lower than 0.1 has very little effect to the stable expression, the final result of regulation is indicated by Boolean value.

The format of regulatory prediction in “Result”:
Gene_name->Gene_name regulation(+/-)