Modeling and prediction of gene-expression patterns reconsidered with Runge - Kutta discretization Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 9, № 6, 2004

MODELING AND PREDICTION OF GENE-EXPRESSION PATTERNS RECONSIDERED WITH RUNGE - KUTTA DISCRETIZATION

T. Ergenc Department of Mathematics METU Middle East Technical University, Ankara, Turkey e-mail: tanil@metu.edu.tr

G.-W. Weber Institute of Applied Mathematics METU Middle East Technical University, Ankara, Turkey e-mail: gweber@metu.edu.tr

Многие проблемы в вычислительной генетике и биоинформатике состоят в моделировании и прогнозировании характера экспрессии генов. С помощью конечной серии экспериментов и примеров и с использованием аппроксимации по методу наименьших квадратов мы получаем нелинейную систему обыкновенных дифференциальных уравнений, которая представляет непрерывную по времени динамику генетического процесса, до того как мы преобразуем ее в дискретную по времени систему. В данной работе мы частично усложняем новаторскую работу [11, 12, 30], используя дискретизацию Рунге — Кутты вместо дискретизации Эйлера, особенно правило трапеций. Таким образом мы модифицируем дискретизацию Эйлера, основываясь на дискретной по времени динамике и анализе устойчивости так, что она становится более подходящей для лежащего в основе генетического или медицинского процесса.

Introduction

In this contribution, we focus on the organization of cell-metabolism, i.e., on the concerted action of biochemical reactions [27], being an important expression of organization and adaptation capability of life [20, 28]. Let us for a better understanding recall some genetical foundations

(cf. [12]).

Background from Biology. Information in living cells is considered to flow from DNA (genome) via RNA (transcriptome) to protein (proteome). DNA can be viewed as a linear macro-molecule composed of nucleotides, carrying one of four different organic bases: A, C, G and T. The linear sequence of the four different nucleotides carries the information. Defined sections on the DNA string, the genes, encode for proteins. During transcription the gene sequences on the DNA are transcribed to linear mRNA macro-molecules. RNA widely resembles DNA chemically. Three nucleotides, comprising a codon, encode for one of the twenty amino

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2004.

acids, the building blocks of proteins. The transformation of the nucleic acid code into an amino acid sequence is called translation. Most of the resulting proteins (the group of enzymes) catalyze chemical reactions within the cell. The sum of all these reactions is defined as cell-metabolism. By metabolites we understand the educts and products of enzymatic catalysis, by metabolome the entity of all metabolites.

Cell Metabolism Traced. A general scientific opinion says that the proteome forms a tight network which lays the functional basis for a living organism. Unfortunately, the proteome is still not satisfactorily accessible by experiments. Due to recent advances, however, in whole genome DNA-sequencing and DNA-microarray technology, both the genome and transcriptome can be readily analyzed. Mining the genome sequence for all known enzymes allows to predict the metabolic potential of a cell. The current metabolic state of a cell is reflected in its metabolome (according to biological dogma). Both the entire metabolome and the entire proteome are not readily accessible for investigation with current methodology. Our approach to learn more about how cell-metabolism is based on the transcriptome. Since all proteins are encoded by mRNAs, the knowledge of all expressed genes (the transcriptome) well approximates the knowledge of all synthesized proteins and enzymes.

"Biochips". Molecular biological methods for analyzing the transcriptome traditionally focus on one gene per experiment. This implies that the throughput is very slow. However, DNA-microarray technology allows the presence of thousands of genes to become monitored in one experiment [25]. Herewith, we obtain a detailled and simultaneous picture of the co-expression of thousands of genes. An array is an ordered arrangement of samples. In the case of DNA-microarrays, the sample is DNA, usually attached to a glass support [25]. The sample spot diameter is less than 200 ^m. On a glass support of the size of a microscope-slide up to 25,000 DNA-samples can be spotted. Herewith, the expression of 25,000 genes could be traced in one experiment (unique probes assumed). There are two major applications for the DNA-microarray technology: (a) identification of sequences (genes / gene mutations), and (b) determination of the gene-expression level, i. e., abundance of transcripts. In any of both cases, the underlying principle is matching known to unknown DNA-samples. The physicochemical matching process hybridization is based on complementary DNA base-pairing rules (i.e., AT and G-C). The basis for a genome wide gene-expression analysis is a completely annotated genome sequence. Based on this information gene specific DNA-probes can be generated (dozens to several thousand nucleotides long). These probes are then spotted onto the glass support. To analyze the abundance of gene transcripts (mRNA or ncRNA), they are isolated and reverse transcribed into cDNA. During the latter step, in general, the cDNA is labelled with a fluorescent dye. The labelled cDNA-population is then hybridized to the DNA-probes on the microchip. Then, the DNA-microchip becomes laser scanned. Thus, depending on the fluorescence intensity of individual spots (usually, a 32 bit message), conclusions about the abundance of the corresponding transcript can be drawn. Finally, all expression data from a single DNA-microarray experiment are represented in a vector, containing as many elements as genes were analyzed.

All mRNA species of the test material (bacteria, tissue, etc.) are isolated, reverse transcribed into cDNA, and labelled with a fluorescent dye. This sample is hybridized with the DNA-microarray (also called DNA-chip). The probes on the DNA-microarray are derived from gene sequence data of the test material. Each individual spot on the DNA-microarray carries millions of copies of one specific probe. Only the cDNA with the same sequence as the corresponding probe can hybridize. Since the samples are labelled with a fluorescent dye, hybridization becomes visualized by laser scanning. Computer assisted image processing, including data

mining, finally reveals the quantity of each individual mRNA species in the test material.

1. Biological Gene-Expression Data Analysis and Genomic Stability

Stability in its mathematical meaning from dynamical systems' theory refers to stationary points (equilibria). Here, it means that being in some sufficiently small neighborhood of the equilibrium, the states never escape [2, 18]. Asymptotic stability is stronger by implying attraction by the equilibrium. In our research, the difference between stability and asymptotic stability is represented by the largest absolute value of a matrix eigenvalues being < 1 or even < 1. In these matrix terms, the step-wise multiplicative form of our time-discrete system will be evaluated by the dynamics' boundedness.

In population biology, stability is well defined [17]. In molecular biology, stability usually refers to the resistance of chemical compounds towards conformational changes and is associated with thermodynamics. Genomic stability is defined as the capability of an organism to repair physical and chemical damages and changes of the genome. Cancer is usually initiated by genome instability, ranging from elevated mutation rates to gross chromosomal rearrangements and alterations in chromosome number. Generally, healthy cells have access to a sophisticated molecular toolbox to fight against instability [16]. For our model we consider gene-expression to be stable, if no strongly emerging changes in time can be observed by DNA-microarray experiments. We, however, allow for changes during the transition from one stable state to another. Those transitions are typically observed when cells are shifted from one metabolic state to another; e.g., the change from feeding on glucose to feeding on glutamate. In contrast, for example, tumor or apoptotic (dying) cells usually show a constant change in gene-expression. These changes are, from a biological point of view, unstable.

2. Time-Continuous and Discrete Models of Gene-Expression Data by Euler Discretization

Euler Discretization. In [12], by Euler's discretization applied on the time-continuous differential equations

(CE) E = M (E)E, (2.1)

which will be the outcome of our mathematical modeling and describing the metabolic process (cf. Paragraph 2), we get the vectors Ek (k £ IN0) of length n. By these Ek coming from time-discretization, we approximate the values E(tk) of the solution trajectories E(•) at times tk. These times may, but need not, be chosen as the discrete times ik (k £ {0,1,...,£}) of biological measurements (see Paragraph 2). Our discrete dynamics can also directly be obtained by inference (statistical learning [1, 14]).

Let E = E(t) be gene-expression patterns at different times t and (CE) our model describing the continuous process of the metabolism of a healthy cell. By using the following algorithm we obtain parameters, for which our system behaves in a stable way. By Euler discretization, for all k £ 1N0 we get:

Ek+1 — Ek

hk

M(Ek)Ek, (2.2)

delivering the sequence

Ek+1 = (I + hk M (Ek ))Ek, (2.3)

where hk = tk+1 — tk, and tk means the k-th time. Let us define

Mk := I + hkM(Ek) (2.4)

so that we obtain the following time-discrete equation and dynamics:

(DE) Ek+i = Mk Ek (k e No). (2.5)

We will analyze whether a (large approximate) set M = {M0, M1,..., Mm-1} of m matrices is stable or not. Our system's dynamics is strongly related to both the structure of the polytopes and the relevant matrices. We present a significant link between the geometry of polytopes and the theory of dynamical systems. We employ a slightly corrected version of the algorithm of R.K. Brayton and C.H. Tong [5]. A first application of this was given by S.W. Pickl [22] on environmental protection. For related attempts cf. [1, 7-10, 19, 24].

Underlying Mathematical Modeling. In [12], we justified our mathematical model (CE) IE = M(E)E. Our data at times tj do not only consist of the measurements Ek, but also of the tendencies of increase or decrease given, e.g., by the difference quotients

Ek

Efc+1 — Ek

tfc+1 — tk

Ef —

U — U-i

if k e {0,...,i- 1},

(2.6)

if k = i.

To point out the basic idea, we refer to the easiest but roughest case given by polynomial

regression. Let E be a vector of the type E = (E1,E2)T. Now, based on a diagonal form of

M(E) = Mai.-.-.aa (E), our ansatz looks, e.g., bi,...,b 6

aE2 + a2E1E2 + asE2 + aAEi + a5E2 + a6 0 1 (2 7)

0 bE + 62E1E2 + 63E2 + 64E1 + 65E2 + b6 I' ( 7)

Now, let us look at some values aj = a*, bj = b*, j e {4, 5, 6} as parameters reserved for stability analysis. Then, the discrete approximation problem takes the form [11, 12]:

... £

minimize ^^ .. . .^ ^ . Il2 .

(p) a a a b b b Ell(M-i-2-3-J-S-S (Ek)Ek — Ek)|2 , (2.8)

a1 , a2, a3 , b1 , b2, b3 k=0 b1,b2,b3,bJ,bJ,bJ

with || • || being Euclidean norm. The solutions generically are locally defined functions (aj(a4, a* , a* , b* , b* , b£) (j e {1 ,2 , 3}) , bj(a4, a5 , a* , b4 , b5 , b6) (j e {1 ,2 ,3})). Inserting them into Ma1,a2,a3,aj,a5,aj (E) gives M(E), still parametrically depending on a*, b*, j e {4, 5,6}. Herewith,

bi,b2,b3,b|,b*,bS

we interpret our entire problem as a two-stage problem [15]. Both classes of parameters can be interpreted in terms of statistical learning [14] as training set and testing set, respectively.

Stability of Matrices in Terms of Polytopes

Definition 2.1. [5] A given (n x n)-matrix M is called stable, if there exists a K £ R+ such that for all j £ IN, ||Mj || < K (|| • || being some natural matrix norm).

Thus, the matrix M is stable if and only if |A(M)| < 1 and in the case of |A(M)| = 1 the algebraic multiplicity is equal to the geometric multiplicity.

A set M of (n x n)-matrices is stable, if for every neighborhood of the origin U C C" there exists another neighborhood of the origin U such that, for each M £ M' it holds MU C U (M' denoting the semi-group generated by all finite products of element of M).

In the pioneering work [12], this stability condition was versified or falsified, respectively, with the help of the algorithm of Brayton and Tong [5], modified and implemented by Pickl [22]. In fact, this algorithm became applied on a finite set M of approximative matrices. A central idea of the underlying construction princle consists in applying any finite number of these matrices on a polyhedral neighbourhood of the origin 0 £ R" or, equivalently, by observing a sequence of stepwise defined polyhedra. The boundedness or unboundedness of the sequence of the latter ones just characterizes the stability or instability of the set M and, herewith, (in)stability of our time-discrete dynamical system (DE).

For numerical calculations based on all the preceding explanations we refer to [11, 12, 30]. In case of Euler discretization, both the mathematical modeling and the stability investigation have turned out to be very elegant and implementable by means of our contruction principle and algorithm. A structural frontier, however, consists in the fact that under Euler's discretization principle, stability of the time-continuous system (CE) and of the time-discrete system (DE) are not always equivalent [18]. We overcome this problem by obtaining (DE) directly by statistical learning from the experimental data, herewith avoiding a model by (DE) [30], or, alternatively, by turning to another discretization principle. This will be initiated and discussed by us in the following.

3. Time-Continuous and Discrete Models of Gene-Expression Data by Trapezoidal Rule

In [12], the mathematical model of metabolic process is described by the continuous differential equation

(CE) E = M(E)E, (3.1)

where E is vector of length n and M(E) is an (nx n)-matrix. This matrix M(E) is constructed by firstly discretizing the right-hand side of (CE) according to Euler:

Ek+! - Ek = M(Ek)Ek (3.2)

hk

and, then, the regarded parameters in M(E) are estimated using least squares approximation so that this discrete model best fits the data Ek (k £ {0,1,...,£}). Lateron, it becomes requested that the time-disrete system behaves stable with respect to the remaining so-called em (expression-metabolism) parameters. Indeed, stability is requested for the discrete equation

(DE) Ek+i = Mk Ek,

(3.3)

where Mk := I + hkM(Ek). The considered set M of matrices Mk is stable if and only if for all k it holds |A(Mk)| < 1, where in the case " = 1" additionally the algebraic and the geometric multiplicity must be the same. Here, A := A(Mk) = 1 + hk^(M(Ek)) stands for all eigenvalues of Mk, while ^ := ^(M(Ek)) represents all the eigenvalues of M(Ek). Within of our spectral way of stability characterization, let us concentrate on the generic case " < 1". This condition is satisfied if hk^ lies in a unit circle of the complex plane C , centered at the point (-1, 0). Since the order of Euler's method is one and its stability region is restricted, we can apply Runge — Kutta methods to construct a time-discrete mathematical model. We will illustrate the idea using Trapezoidal rule which is a second order 2-stage Runge — Kutta method and has stronger stability properties (for general information cf. [13]).

Discretization of equation (CE) by Trapezoidal rule yields

(TR) Efc+1 = Ek + 1 hfc(M(Efc+1)Efc+1 + M(Efc)Efc) (k G Nq). (3.4)

Therefore, the matrix M(E) which best fits the measurements Ek (k G {0,1,... ,£}) will be

constructed by minimizing

'-1 ~ 1 ^ ^ ^

J] IE + ^hk(M(Ek+1)Ek+1 + M(Ek)Ek) - Ek+1)||2 (3.5)

k=0

with respect to all considered parameters (|| • || denoting Euclidean norm; please cf. also (P)).

It can be easily seen that for the linear model E = ME, where the constant matrix M may be the Jacobian coming from a linearization of equation (CE), our equation (TR) in a generical sense reduces to the formula of a time-discrete system

(DE) Ek+1 = Mk Ek (3.6)

with the system matrices given by

1, .A-1 M 1

Mk = - 2hkMJ + 2hkMJ . (3.7)

Since for these matrices the eigenvalues are

1 +1 hk mm )

A(Mk) =-1-, (3.8)

1 - 2 hk MM)

the stabilty condition |A(Mk)| < 1 is satisfied with no restriction on hkMM), i.e., the system is unconditionally stable.

The stability investigation of (TR) for the nonlinear case of a nonconstant function M(E) and more general versions of Runge — Kutta methods as well are in the future research plan on genetical processes which we have together with our colleagues. This plan also includes the more general structure given by an additional shift vector B on the right-hand side of (CE) as introduced and successfully studied by F.B. Yilmaz [30], and using another methods of statistical learning, too. We conclude with a further encouraging remark.

Remark 3.1. Runge — Kutta discretization of model equation (CE) generates a nonlinear discrete equation for parameters, even if the model equation is linear. If we use implicit Runge —

Kutta methods, depending on the base functions used in construction of the matrix M(E), derivation of discrete equation may not be possible. We illustrate the nonlinearity by applying 2-stage explicit Runge — Kutta method to the linear model equation E = ME, where there is a constant matrix (to be determined). Discretization gives

where a + a2 obtain

which means

(RK)

1, ßa

Ek+i 1 2

Ek + hk(aiKi,k + a2K2,fc) (k £ No), (3.9)

and K1;k := MEk, K2,k := M(Ek + ßhkK1;k). Therefore, we

E

fc+i

(DE )2

Ek + hk MEk + hk M 2Ek,

E

k+i

2

Mk Ek

(3.10)

(3.11)

h2

with Mk := I + hkM + -2M2. Since the entries of M are our parameters in the mathematical

modeling problem, the term M2 represents the parametrical nonlinearity (quadratic regression). Please note that, nevertheless, for the time-discrete dynamics (DE)2 with its system matrices Mk the stability analysis of [12] including its algorithm can be carried over and applied!

Conclusion

This contribution introduces a new method for the dynamical analysis of the transcriptome which plays a crucial role in orchestrating the cell+s metabolism. The analysis of gene-expression opens new horizons, since the knowledge on regulatory non-coding RNA grows very fastly. In contrast to related approaches, we continue analyzing the system for stability. The innovation, however, consists in the utilization of Runge — Kutta methods. We discussed their advantages in terms of stability, their analyical structural frontiers and numerical challenges which are a topic of future research. Finally, we hope that our widened methodic concept will become enriching and giving new insights into cell metabolism and its regulation.

Acknowledgement: We express our gratitude to our colleagues and friends Prof. Dr. Marat Akhmet, Dipl.-Math. Jutta Gebert, Dipl.-Math. Martin Latsch, Prof. Dr. Bulent Kara-sozen, Dr. Hakan Oktem, Dr. Stefan W. Pickl, Prof. Dr. Karl Roesner, Dr. Nina Shokina, Prof. Dr. Yurii Shokin, Dr. Robbe Wunschiers and MSc F. Bilge Yilmaz for valuable discussions and encouragement.

References

[1] Akhmet M.U., Gebert J., Öktem H., Pickl S.W., Weber G.-W. An improved algorithm for analytical modeling and anticipation of gene expression patterns // J. Computational Technologies. 2004 (in press).

[2] Amann H. Gewöhnliche Differentialgleichungen. Berlin, N.Y.: Walther de Gruyter. 1983.

[3] Aster A., Borchers B., Thurber C. Parameter Estimation and Inverse Problems. N.Y.: Acad. Press, 2004.

[4] Bank B., Guddat J., Klatte D. et al. Non-Linear Parametric Optimization. Berlin: Akad. Publ. House, 1982.

[5] Brayton R.K, Tong C.H. Stability of dynamical systems: A constructive approach // IEEE Trans. on Circuits and Systems. 1979. Vol. 26, N 4. P. 224-234.

[6] Brayton R.K, Tong C.H. Constructive stability and asymptotic stability of dynamical systems // IEEE Trans. on Circuits and Systems. 1980. Vol. 27, N 11. P. 1121-1130.

[7] Chen T., He H.L, Church G.M. Modeling gene expression with differential equations // Proc. Pacific Symp. on Biocomputing. 1999. P. 29-40.

[8] D 'Haeseleer P. Reconstructing gene networks from large scale gene expression data. Univ. of New Mexico/USA, 2000.

[9] De Hoon M., Imoto S., Miyano S. Inferring gene regulatory networks from time-ordered gene expression data using differential equations // Lecture Notes in Computer Sci. 2002. N 2534. P. 267-274.

[10] De Hoon M., Imoto S. Kobayashi K. et al. Inferring gene regulatory networks from time-ordered gene expression data of Bacillus subtilis using differential equations // Proc. Pacific Symp. on Biocomputing. 2003. P. 17-28.

[11] Gebert J., LAtsch M., Pickl S.W., Weber G.-W., WUnschiers R. Genetic networks and anticipation of gene expression patterns, invited paper, Computing Anticipatory Systems: CASYS'03 // Sixth Intern. Conf., AIP Conf. Proc. 2004. N 718. P. 474-485.

[12] Gebert J., LAtsch M., Pickl S.W., Weber G.-W., WUnschiers R. An algorithm to analyze stability of gene-expression pattern, to appear in the special issue // Discrete Mathematics and Data Mining of Discrete Applied Mathematics / E. Boros, P.L. Hammer, A. Kogan (Eds). 2004.

[13] Hairer E., Lubich C., Wanner G. Geometric Numerical Integration. N.Y., Berlin: SpringerVerlag, 2002.

[14] Hastie T., Tibshirani R., Friedman J. The Elements of Statistical Learning. N.Y., Berlin: Springer-Verlag, 2001.

[15] Jongen H.Th., Weber G.-W. On parametric nonlinear programming // Annals of Operations Research. 1990. Vol. 27. P. 253-284.

[16] Kolodner R.D., Putnam C.D, Myung K. Maintainance of genome stability in Saccharomyces cerevisiae // Science. 2002. Vol. 297. P. 552-557.

[17] Kot M. Elements of Mathematical Ecology. Cambridge Univ. Press, 2001.

[18] Krabs W. Dynamische Systeme: Steuerbarkeit und chaotisches Verhalten. Teubner, Stuttgart, Leipzig, 1998.

[19] Mjolsness E., Mann T., Castano R., Wold B. From coexpression to coregulation: An approach to inferring transcriptional regulation among gene classes from large-scale expression data // Advances in Neural Information Proc. Systems. 2000. Vol. 12. P. 928-934.

[20] Monod J. Chance and Necessity. Random House, 1972.

[21] Palis J., de Melo W. Geometric Theory of Dynamical Systems. N.Y., Berlin: Springer-Verlag, 1982.

[22] PlCKL S. The t-value als Kontrollparameter. Modellierung und Analyse eines Joint-Implementation Programmes Mithilfe der Kooperativen Dynamischen Spieltheorie und der Diskreten Optimierung. Aachen: Shaker Verlag, 1999.

[23] RoCKAFELLAR R.T. Convex Analysis. Princeton Univ. Press, 1970.

[24] Sakamoto E., Iba H. Inferring a system of differential equations for a gene regulatory network by using genetic programming // Proc. Congress on Evolutionary Computation. 2001. P. 720-726.

[25] SCHENA M. DNA Microarrays. Oxford Univ., 1999.

[26] SCHENA M. Microarray Analysis. John Wiley & Sons, 2003.

[27] Schilling C.H., Letscher D, Palsson B. Theory for the systematic definition of metabolic pathways and their use in interpreting metabolic function from a pathway-oriented perspective // J. Theor. Biol. 2000. Vol. 203. P. 229-248.

[28] Schrödinger E. What is life? Cambridge Univ. Press, 1944.

[29] Weber G.-W. Generalized Semi-Infinite Optimization and Related Topics. Heldermann Verlag, 2003.

[30] Yilmaz F.B. A Mathematical Modeling and Approximation of Gene Expression Patterns by Linear and Quadratic Regulatory Relations, MSc Thesis at Institute of Applied Mathematics, METU, Ankara, 2004.

Received for publication October 18, 2004