Volume 10 Supplement 3

## Eighth International Conference on Bioinformatics (InCoB2009): Computational Biology

# Detecting robust time-delayed regulation in *Mycobacterium tuberculosis*

- Iti Chaturvedi
^{1}Email author and - Jagath C Rajapakse
^{1, 2, 3}Email author

**10(Suppl 3)**:S28

https://doi.org/10.1186/1471-2164-10-S3-S28

© Chaturvedi and Rajapakse; licensee BioMed Central Ltd. 2009

**Published: **3 December 2009

## Abstract

### Background

Time delays are often found in gene regulation though most techniques of building gene regulatory networks are not capable of capturing such phenomena. Here we look at the delays in the DNA repair system of *Mycobacterium tuberculosis* which is unusually slow in the bacteria. We propose a method based on a skip-chain model to study this phenomena in gene networks. The Viterbi paths of the underlying Markov chains find the most likely regulatory interactions among genes, taking care of very long delays. Using the derived networks, we discuss the delayed regulations and robustness of the DNA damage seen in the bacterium.

### Results

We evaluated our method on time-course gene expressions after DNA damage with Mitocyin C. Several time-delayed interactions were observed with our analysis. The presence of hubs in the networks indicates that a small number of transcriptional factors regulate the rest of the system. We demonstrate the use of priors to overcome over-fitting problem in the generation of networks. We compare our results with the gene networks derived with dynamic Bayesian networks (DBN).

### Conclusion

Different transcription networks are active at different stages, and constant feedback and regulation is maintained throughout the activities of a biological pathway. Skip-chain models are capable of capturing, long distant and the time-delayed regulations. Use of a Dirichlet prior over parameters and Gibbs prior over structure can greatly reduce the over-fitting in the new model.

## Keywords

## Background

Cellular activities of genes and gene products represented in gene regulatory networks (GRN) provide a basis for signal transduction pathways. Since the signal transduction is transient, the study of dynamics of the transduction is essential. Further, the *distributed* nature of cell fate regulation events manifest's itself as intense crosstalk between the nominal pathways. States of gene networks are often presumed to be stable, meaning that slight changes in the state of a few parents do not change the expression state of the child gene. This phenomena relates to the redundancy of biological systems which are to ensure that the system retains functioning inspite of the perturbations.

In this work, we use Bayesian networks (BN) in the stochastic framework to represent GRN. Pathways have a natural representation of BN, where genes are nodes in the network and edges are causal interactions among them. The causal dependencies are given as conditional probabilities which infer 'cause and effect' relationships among genes in the network. A BN being acyclic is not able to model feedbacks and self-regulation events. The dynamic Bayesian network (DBN) is defined by a pair of structures (*S*_{
t
}, *S*_{t+1}) each corresponding to time instances *t* and *t* + 1 and a transition network of interactions between the two networks [1]. DBN assumes that the genetic regulation process is first-order Markovian where parents are from the previous time point and can allow cyclic events.

However, several time-delayed interactions are known to exist in biological systems. DBN was extended to a higher-order where mutual information (MI) has been used to determine the best time-delay of an interaction [2]. However, these generative models become intractable at very high orders, so we resort to a conditional skip-chain model. In a skip-chain model, the linear features model the lower-order delays and the skip features model long-distant delays [3].

The linear feature attempts to model interactions which occur instantly or with little delay. The skip feature model interactions occurring much later in the pathway, for example, a gene *g*_{
i
}inhibits a gene *g*_{
j
}to start a process, and later *g*_{
i
}regulates another gene *g*_{
k
}towards the end of the process. The skip-feature probability is decomposed into a sum of terms for consecutive pairs of genes in the time-course and the most likely interactions are found using the Viterbi algorithm. The Viterbi skip-feature can automatically determine the best time delay in a higher-order Markov chain representing the instantaneous network of DBN.

Our approach consists of three stages: first, our method involves identifying time-delayed interaction features and predicting the optimal GRN by using a GA. The fitness function of the GA is modified to include Viterbi scores of time-delayed interactions by using the skip-chain model. Next, an application to DNA repair system of *Mycobacterium tuberculosis* has been performed. This bacteria causes tuberculosis in man and is known to have a very slow growth rate *in vitro*. In particular, we consider the DNA repair pathway which is activated when a damage to the DNA occurs. The system consists of proteins *LexA* and *RecA* as well as up to 40 genes that are regulated by these two proteins together. Lastly, we discuss our findings and directions for future work.

## Methods

where *x* = (*x*_{1}, *x*_{2}, ...., *x*_{
n
}), the conditional probability of gene expression *x*_{
i
}given its parents *a*_{
i
}is *p*(*x*_{
i
}|*a*_{
i
}, *θ*_{
i
}), and *θ*_{
i
}denotes the parameters of the conditional probabilities.

*S*

_{ t },

*S*

_{t+1}) corresponding to time instances

*t*and

*t*+ 1. The DBN structure is obtained by unrolling the transition network over time. In time instance

*t*+ 1, the parents of genes are those specified in the time instant

*t*. The likelihood of transition network

*S*of interactions between time instances

*t*and

*t*+ 1 is given by

where
correspond to the number of instances of *θ*_{
ijk
}= *p*(*x*_{i, t+1}= *k|a*_{i, t}= *j*), *k* is the discretized gene expression level, and *j* is the discrete state combination of parent genes. The first-order DBN has two layers of genes, and therefore 2*n* nodes.

The classical DBN is unable to capture complex time-dependencies and is extended to an *o*-order Markov chain (*o* ≥ 2). It predicts the expression levels of a set of genes based on expression upto previous *o* time points. However, such an approach cannot handle long range dependencies because as *o* increases the search space becomes intractable. Instead, we employ skip-chain models which augments linear chain features that represent local features, with skip-features representing long range dependencies [4, 5]. It then simply factorizes the prediction probabilities into linear and skip features.

*f*(

*x*

_{ i },

*a*

_{i(t-o:t)},

*t*) represent local dependencies that are consistent with an

*o*- order Markov assumption of gene expressions. But for long distant interactions, we relax this assumption by using skip-chain feature functions

*h*(

*x*

_{ i },

*a*

_{ i },

*s*

_{ t },

*t*) which exploit dependencies between genes that are arbitrarily distant at time instances

*s*

_{ t }and

*t*, respectively (Fig. 1). Such a skip-feature models variable length Markov chain upto

*m*- 1 order where

*m*is the number of time points.

*x*

_{ i }is a weighted sum of linear and skip-edge scores:

where *λ* ≤ 1 is a weight determined heuristically.

For interactions, we look for causal effects of regulated genes as features. We can use the Viterbi algorithm to find a maximum likelihood (ML) path between two genes at distant time points in a hidden Markov model (HMM) [7]. The ML can then be used to make a choice between different time-delayed interactions of the same pair of genes. For any two genes *g*_{
i
}and *g*_{
j
}, we choose the highest Viterbi score among all the possible interaction features.

A genetic algorithm is used to find the optimal network structure. Here an individual is defined by matrix {*c*_{i, j}}_{n × n}with dimension *n* × *n*. Each cell *c*_{i, j}is randomly initialized with interactions which have MI at a time lag *o* above a threshold. Here *g*_{
j
}is the parent of *g*_{
i
}. The GA then finds the structure with the highest posterior probability (Eq. 3). The GA provides an optimal structure maximizing the likelihood asymptotically. We also explored the use of two priors over the network.

### Dirichlet prior over parameters

Most higher-order Markov models are far from optimal. They are extremely sensitive to change in pathways and associated data. This happens as most of the data is general rather than feature specific for an interaction. The goal of adaption has been to make good use of available feature data and reduce the over-fitting in the model. Our adaption model combined the reliable general DBN with a volatile feature specific HMM for long delays. We further extend the MLE to a Bayesian learning where a Dirichlet conjugate prior is imposed on each of the parameters.

*θ*= {

*θ*

_{ i }:

*i*= 1, 2, ...

*n*}, the likelihood can be written as

where *h'*(*x*_{
i
}) is total probability of the skip-path, *α* is a weighting factor between the linear and skip features.

*g*

_{ i }based on linear and skip-edges.

here, instead of using a constant, *λ* is estimated using prior linear feature and the total probability of the skip path.

### Gibbs prior over graph

*P*(

*S*) of the gene network. A Gibbs distribution takes the form of

*P*(

*S*) ∝

*e*-

^{E(S)}where energy of the graph

*E*(

*S*) can be factorized into a sum of interaction potentials

*U*

_{ ij }between genes

*g*

_{ i }and

*g*

_{ j }. If an interaction exists in the target network, we set

*U*

_{ ij }=

*σ*

_{1}otherwise

*U*

_{ ij }=

*σ*

_{2}. The total energy of the graph over existing edges is

*E*(

*S*) = ∑

_{{i, j} ∈ S}

*U*

_{ ij }. The posterior probability of the graph is then given by

A small *σ*_{1} and a large *σ*_{2} will reflect the prior target network more in the GRN and vice-versa.

## Experiments and results

We evaluated our method on a DNA repair system of *Mycobacterium tuberculosis* by building regulatory networks with DBN, HDBN, and skip-chain model. Here we looked at the response of bacteria to drug-induced stress. Treatment with Mitomycin C caused DNA damage and hence led to the upregulation of associated repair genes. Eight time points are available at NCBI Gene Expression Omnibus (GSE1642-GPL1396 series) 0.33 hr, 0.75 hr, 1.5 hr, 2 hr, 4 hr, 6 hr, 8 hr and 12 hr after DNA damage. The data was discretized into 1 for upregulation and 0 for downregulation by using an approach described previously [9].

The corresponding skip probabilities were calculated as described in methods. Upto seven time points of delays were allowed. Firstly, we used 9 genes previously specified [10]. In order to get an expanded dataset, the original dataset was subjected to ICA and the components closest to 9 genes were identified [11, 12].

Time-delayed interactions in predicted network

Higher-order edges | |||||||
---|---|---|---|---|---|---|---|

# Genes | Model: o | ML | 1 | 2 | 3 | 4 | 5 |

9 | DBN:1 | -14.7 | 9 | ||||

HDBN:3 | -8.69 | 8 | 2 | 7 | |||

SKIP-CHAIN:1 | -6.05 | 13 | (3) | ||||

32 | DBN:1 | -48.9 | 36 | ||||

HDBN:4 | -39.4 | 20 | 6 | 14 | 20 | ||

SKIP-CHAIN:2 | -37.2 | 54 | 18 | (41) | (4) |

Time-delayed interactions in predicted network using prior

Higher-order edges | |||||||
---|---|---|---|---|---|---|---|

# Genes | Model: o | ML | 1 | 2 | 3 | 4 | 5 |

9 | SKIP-CHAIN:1 | -6.05 | 13 | (3) | |||

SKIP-CHAIN(Gibbs):1 | -5.8 | 11 | (2) | ||||

SKIP-CHAIN(Dirichlet):2 | -5.2 | 7 | 13 | (11) | (1) | ||

SKIP-CHAIN(Gibbs and Dirichlet):3 | -3.27 | 2 | 7 | (4) | (5) | ||

32 | SKIP-CHAIN:2 | -37.2 | 54 | 18 | (41) | (4) | |

SKIP-CHAIN(Gibbs):3 | -35.7 | 37 | 16 | 24 | (40) | (3) | |

SKIP-CHAIN(Dirichlet):2 | -35.05 | 54 | 16 | (37) | (4) | ||

SKIP-CHAIN(Gibbs and Dirichlet):2 | -34.54 | 50 | 15 | (41) | (4) |

The presence of hubs or single genes regulating several other genes are also seen in the network. These networks can buffer environmental variations. It can be seen that a small number of transcription factors (TF) regulate the rest of the repair system. At the same time the in-degree is low, as each gene is regulated by just one TF. RecA causes inactivation of lexA which suppresses DNA repair genes. We observe binding of recA(DNA repair) to dnaB(DNA replication) helicase. RecA also activates linB which causes dehalogenation needed for transformation events in dna repair. The Fadd genes initiate apoptosis and are also required for cell-wall formation.

The second dataset of 32 genes indicated that our method is good for identifying core genes (Fig. 4). RecA and lexA are shown to be critical hub by both DBN and HDBN. The HDBN showed several time-delayed interactions at 2 and 4 hrs. The skip-chain gave a fewer interactions though it also showed interactions at 6 hrs. Use of prior gives better networks with few hubs in Fig. 4(f). They could detect new hubs like ruvC, fadD21 and fadD23.

## Discussion and conclusion

An organism responds to changes in its environment by altering the level of expression of critical genes. The virulence of *Mycobacterium tuberculosis* depends on the ability of the bacilli to switch between replicative (growth) and non-replicative (dormancy) states in response to host immunity. Different transcription networks are active at different stages of the response. The coordinated repression of genes are likely to contribute to survival by conserving energy and precursors under nutrient-limiting conditions and/or minimizing expression of potential antigens.

*M. tuberculosis* is known to have an unusally long period of 10 hrs for the DNA replication fork to traverse the chromosome. Our results showed several interactions at 4 hrs in the DNA repair pathway. An order-4 HDBN with skip-chain dependencies was shown to outperform ordinary HDBN's. For genes to interact they both have to be upregulated. We use this property to select events where a pair of genes are both upregulated at similar or delayed time points. It is well established that interacting genes have correlated expression patterns. To this end, we add the interactions at non-consecutive time points. This is because a DBN assumes a first-order network and is not able to model complex time-delayed interactions. We assumed that all interactions had equal priors. However our method is able to distinguish between short- and long-term interactions and hence allow us to make a better judgement on DNA repair.

To include time-delays, we used a skip-chain model. The Viterbi shortest path allowed us to choose between time delayed interactions of two genes of same and different time delays. This lets us identify the best interaction information from the dataset. By using a single parent Viterbi path to model the upregulated events, we were able to focus on special cases in the DBN. This significantly reduces the search space for the GA. Our search is however constrained by various parameters like MI and number of parents.

Skip-chain models address the difficulties of a DBN by easily incorporating overlapping input features. We also see that using approximate inference leads to lower total training time without loss in accuracy. The skip-chain BN is not an HDBN because usually different long-distance dependencies are used by skipping the intermediate time points. We proposed a method that can extract long distant regulations and demonstrated it on DNA repair of tuberculosis. Our approach may be useful for understanding complex gene regulation mechanisms.

Lastly, using priors gave us higher likelihood and improved the over-fitting in building the regulatory networks. The Dirichlet prior gave fewer hubs as compared to the Gibbs prior and gave a higher likelihood. The combination of the two priors gave us the best regulatory networks. We can see that the prediction with prior allows higher-orders of linear model aswell.

## Note

Other papers from the meeting have been published as part of *BMC Bioinformatics* Volume 10 Supplement 15, 2009: Eighth International Conference on Bioinformatics (InCoB2009): Bioinformatics, available online at http://www.biomedcentral.com/1471-2105/10?issue=S15.

## Declarations

### Acknowledgements

The authors wish to acknowledge the partial support from Interdisciplinary Research Group on Infectious Diseases at the Singapore-MIT Alliance Research and Technology (SMART) Center.

This article has been published as part of *BMC Genomics* Volume 10 Supplement 3, 2009: Eighth International Conference on Bioinformatics (InCoB2009): Computational Biology. The full contents of the supplement are available online at http://www.biomedcentral.com/1471-2164/10?issue=S3.

## Authors’ Affiliations

## References

- Friedman N, Murphy K, Russell S: Learning the Structure of Dynamic Probabilistic Networks. Proceedings of the 14th Annual Conference on Uncertainty in Artificial Intelligence (UAI-98). 1998, 139-14.Google Scholar
- Zhengzheng X, Dan W: Modeling Multiple Time Units Delayed Gene Regulatory Network Using Dynamic Bayesian Network. Data Mining Workshops, 2006. ICDM Workshops 2006. Sixth IEEE International Conference on. 2006, 190-195.Google Scholar
- Chaturvedi I, Rajapakse J: Fusion of Gene Regulatory and Protein Interaction Networks Using Skip-Chain Models. Pattern Recognition in Bioinformatics. 2008, Lecture Notes in Computer Science, Springer Berlin/Heidelberg, 5265: 214-224. full_text.View ArticleGoogle Scholar
- Galley M: A Skip-Chain Conditional Random Field for Ranking Meeting Utterances by Importance. Proceedings of the 2006 Conference on Empirical Methods in Natural Language Processing (EMNLP 2006). 2006, Sydney: Association for Computational Linguistics, 364-372.View ArticleGoogle Scholar
- Sutton C, McCallum A: Collective Segmentation and Labeling of Distant Entities in Information Extraction. Presented at ICML 2004 Workshop on Statistical Relational Learning. 2004Google Scholar
- Fink GA, service SO: Markov Models for Pattern Recognition From Theory to Applications. 2008Google Scholar
- Hao T, Huang TS: Improved Graphical Model for Audiovisual Object Tracking. Multimedia and Expo, 2006 IEEE International Conference on. 2006, 997-1000.Google Scholar
- Shuanhu B, Haizhou L: Bayesian Learning of N-Gram Statistical Language Modeling. Acoustics, Speech and Signal Processing, 2006. ICASSP 2006 Proceedings. 2006 IEEE International Conference on. 2006, 1: I-IGoogle Scholar
- Shmulevich I, Zhang W: Binary analysis and optimization-based normalization of gene expression data. Bioinformatics. 18 (4): 555-65. 10.1093/bioinformatics/18.4.555.Google Scholar
- Gebert J, Motameny S, Faigle U, Forst CV, Schrader R: Identifying Genes of Gene Regulatory Networks Using Formal Concept Analysis. Journal of Computational Biology. 2008, 15 (2): 185-194. 10.1089/cmb.2007.0107.View ArticlePubMedGoogle Scholar
- Francis RB, Michael IJ: Kernel independent component analysis. J Mach Learn Res. 2003, 3: 1-48. 10.1162/153244303768966085.Google Scholar
- Suri RE: Application of independent component analysis to microarray data. Integration of Knowledge Intensive Multi-Agent Systems, 2003. International Conference on. 2003, 375-378.Google Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.