Katja Rateitschak
email: katja.rateitschak@uni-rostock.de

Olaf Wolkenhauer
email: olaf.wolkenhauer@uni-rostock.dewww.sbi.uni-rostock.de

Thresholds in transient dynamics of signal transduction pathways

Keywords: signal transduction, transient dynamics, Lyapunov exponents, metabolic control coefficient, Hill equation, relative instability exponent

Thresholds in transient dynamics of signal transduction pathways

Katja Rateitschak, Olaf

A c c e p t e d m a n u s c r i p t 1 Introduction

The behavior of cellular processes can be classified into steady state and transient dynamics, whereby we include biochemical oscillations with constant amplitude as examples of steady state processes. While steady state processes are most important in metabolic systems [START_REF] Heinrich | The Regulation of Cellular Systems[END_REF], transient dynamics play an important role in cell communication [START_REF] Aldridge | Direct Lyapunov exponent analysis enables parametric study of transient signaling governing cell behavior[END_REF][START_REF] Bolouri | Transcriptional regulatory cascades in development: initial rates, not steady state, determine network kinetics[END_REF][START_REF] Kao | Identification of the mechanisms regulating the differential activation of the mapk cascade by epidermal growth factor and nerve growth factor in PC12 cells[END_REF][START_REF] Cheong | Transient IKK activity mediates NF-κB temporal dynamics in response to a wide range of TNFα doses[END_REF][START_REF] Hao | Regulators of G protein signaling and transient activation of signaling: experimental and computational analysis reveals negative and positive feedback controls on G protein activity[END_REF][START_REF] Sasagawa | Prediction and validation of the distinct dynamics of transient and sustained ERK activation[END_REF].

Transient dynamics describe the transition from some initial state of the system into a steady state. For an understanding of signaling pathways and their targets transient dynamics are often more important than steady state dynamics as demonstrated by the following examples: Aldridge et al. [START_REF] Aldridge | Direct Lyapunov exponent analysis enables parametric study of transient signaling governing cell behavior[END_REF] have established a mathematical model which describes the regulation of active caspase 3 dynamics in apoptosis. The temporal duration depends on the initial condition (concentration) of the protein xiap. Bolouri and Davidson [START_REF] Bolouri | Transcriptional regulatory cascades in development: initial rates, not steady state, determine network kinetics[END_REF] have shown by mathematical modelling that genes are activated successively in a regulatory network cascade, long before steady states are attained. The EGFR stimulated MAPK cascade in the PC12 cell line shows a transient activation of ERK with a peak at 5 min and a return to its basal level after 30 min, which leads to cellular proliferation [START_REF] Kao | Identification of the mechanisms regulating the differential activation of the mapk cascade by epidermal growth factor and nerve growth factor in PC12 cells[END_REF]. Using mathematical modelling Cheong et al. [START_REF] Cheong | Transient IKK activity mediates NF-κB temporal dynamics in response to a wide range of TNFα doses[END_REF] have predicted that the dynamic profile of the IKK signal must transiently peak at all TNFα doses in order to generate the observed NFκB activity which they could experimentally validate.

Aldridge et al. have employed direct finite-time Lyapunov exponents to identify domains in the transient dynamics of high sensitivity to initial conditions [START_REF] Aldridge | Direct Lyapunov exponent analysis enables parametric study of transient signaling governing cell behavior[END_REF]. These separatrices delineate regions with different transient dynamics and can be considered as a threshold. Furthermore, for the activation of a signalling pathway often threshold concentrations of proteins are required to prevent the unintentional activation of signalling pathways through random fluctuations.

The relevance of thresholds in biology is demonstrated by a large number of publications about ultrasensitive responses ranging from Goldbeter [START_REF] Goldbeter | An amplified sensitivity arising from covalent modification in biological systems[END_REF] and Ferrell [START_REF] Ferrell | Tripping the switch fantastic: how a protein kinase cascade can convert graded input into switch-like outputs[END_REF] to the recent work of Buchler [START_REF] Buchler | Protein sequestation generates a flexible ultrasensitive response in a genetic network[END_REF]. Ferrell has concluded, on the basis of experimental data, that the MAPK cascade is optimized to convert a graded input into a switch-like output [START_REF] Ferrell | Tripping the switch fantastic: how a protein kinase cascade can convert graded input into switch-like outputs[END_REF]. Ultrasensitive responses leading to sigmoidal stimulusresponse curves have been related to activation thresholds of biochemical networks in the literature [START_REF] Goldbeter | An amplified sensitivity arising from covalent modification in biological systems[END_REF][START_REF] Ferrell | The biochemical basis of an all-or-none cell fate switch in Xenopus oocytes[END_REF][START_REF] Nash | Multisite phosphorylation of a CDK inhibitor sets a threshold for the onset of DNA replication[END_REF][START_REF] Bhalla | MAP kinase phosphatase as a locus of flexibility in a mitogen activated protein kinase signaling network[END_REF][START_REF] Huang | Ultrasensitivity in the mitogen-activated protein kinase cascade[END_REF][START_REF] Bentele | Mathematical modeling reveals threshold mechanism in CD95-induced apoptosis[END_REF][START_REF] Gunarwardena | Multisite protein phosphorylation makes a good threshold but can be a poor switch[END_REF][START_REF] Salazar | Versatile regulation of multisite protein phosphorylation by the order of phosphate processing and protein-protein interactions[END_REF]. Experimental measurements in networks with multiple phosphorylations have shown that the stimulus has to exceed a threshold concentration to activate downstream events [START_REF] Buchler | Protein sequestation generates a flexible ultrasensitive response in a genetic network[END_REF][START_REF] Ferrell | The biochemical basis of an all-or-none cell fate switch in Xenopus oocytes[END_REF][START_REF] Nash | Multisite phosphorylation of a CDK inhibitor sets a threshold for the onset of DNA replication[END_REF][START_REF] Bhalla | MAP kinase phosphatase as a locus of flexibility in a mitogen activated protein kinase signaling network[END_REF]. Mathematical models for the MAPK cascade with dual phosphorylation [START_REF] Huang | Ultrasensitivity in the mitogen-activated protein kinase cascade[END_REF] and general models for multisite phosphorylations [START_REF] Salazar | Versatile regulation of multisite protein phosphorylation by the order of phosphate processing and protein-protein interactions[END_REF] have shown that multisite phosphorylation give rise to more threshold like responses that single site phosphorylation. So far, however, mathematical modelling has focused on network properties leading to sigmoidal stimulus-response curves in the steady state dynamics.

A c c e p t e d m a n u s c r i p t

In small networks, stimulus-response curves of the transient and steady dynamics can be studied in detail, leading to the identification of putative thresholds. In larger biochemical networks the identification of putative thresholds would become very time consuming because the network components could show sigmoidalstimulus response curves during different time intervals. Therefore quantitative measures to identify putative thresholds in transient dynamics are required. They should be designed such that an inspection of trajectories is not necessary, i. e. they should deliver initial conditions (concentrations) of proteins and respective time intervals where a response shows a threshold. Recent examples of quantitative methods to characterize transient signals include transient metabolic control coefficients (MCC) [START_REF] Ingalls | Sensitivity analysis of stoichiometric networks: an extension of metabolic control analysis to non-steady state trajectories[END_REF][START_REF] Hornberg | Principles behind the multifarious control of signal transduction ERK phosphorylation and kinase/phosphatase control[END_REF] and finite time Lyapunov exponents (FTL) [START_REF] Aldridge | Direct Lyapunov exponent analysis enables parametric study of transient signaling governing cell behavior[END_REF]. We here focus on system structures that generate thresholds in transient dynamics and investigate quantitative measures to identify such thresholds.

The outline of this paper is as follows. Following an introduction to MCCs and FTLs, we define thresholds in cellular signaling, discuss their properties and argue for the need of quantitative measures to identify them in transient dynamics. We then study two model structures that can generate thresholds in their transient response and we compare quantitative measures to identify thresholds in the transient dynamics. We first extend the analysis of thresholds by FTLs in an apoptosis decision network [START_REF] Aldridge | Direct Lyapunov exponent analysis enables parametric study of transient signaling governing cell behavior[END_REF] by calculating the three largest FTLs and the MCCs. Next, we study a gene transcription network as an alternative model structure to generate thresholds in transient dynamics. Our analysis of the gene transcription network has motivated us to analytically quantify the threshold of a Hill equation and to discuss an alternative measure which combines properties of the FTLs and the MCCs.

Finally the results of our comparative study will lead us to an evaluation of the investigated quantitative measures regarding their ability to identify thresholds in transient signaling.

Quantitative measures to characterize transient dynamics

Metabolic control coefficients (MCC) measure the relative response of a state variable x i , with respect to the relative perturbation by state variable x j . They are a standard quantitative measure in sensitivity analysis [START_REF] Heinrich | The Regulation of Cellular Systems[END_REF]. MCCs can also be calculated for finite times as described in [START_REF] Ingalls | Sensitivity analysis of stoichiometric networks: an extension of metabolic control analysis to non-steady state trajectories[END_REF][START_REF] Hornberg | Principles behind the multifarious control of signal transduction ERK phosphorylation and kinase/phosphatase control[END_REF]. They are defined by

C x i (t), x j (0) = ∂ log x i (t) ∂ log x j (0) . (1)
Comparing the influence of a perturbation at t = 0 on different state variables

A c c e p t e d m a n u s c r i p t

at t > 0 it is useful to consider relative changes of state variables because the values of different state variables can have different orders of magnitude. Thus control coefficients for different pairs (i, j) can be compared. In general one uses the same percentaged perturbation for all x j (0) resulting in the same denominator for all C x i (t), x j (0) . To perturb simultaneously initial conditions of several state variables, global approaches have to be applied [START_REF] Zhang | Comparative study of parameter sensitivity analyses of the TCR-activated Erk-MAPK signaling pathway[END_REF], which is not the focus of this work.

The MCC in Eq. ( 1) is also called "concentration control coefficient". Replacing in the denominator the initial condition x j (0) by a parameter, for example a reaction constant, leads to another type of concentration control coefficient. The concentration control coefficients can be calculated for the transient dynamics as well as in the steady state. In addition, flux control coefficients have been introduced to quantify the control of the flux of a metabolic system in the steady state [START_REF] Heinrich | The Regulation of Cellular Systems[END_REF].

Lyapunov exponents measure the exponential divergence between a reference trajectory and d orthogonal perturbations to the trajectory, where d is the dimension of the related mathematical model which is equal to the number of state variables [START_REF] Eckmann | Ergodic theory of chaos and strange attractors[END_REF]. Finite time Lyapunov exponents (FTL) have been introduced to quantify dynamical instabilities over a finite interval of time [START_REF] Grassberger | Scaling laws for invariant measures on hyperbolic and nonhyperbolic attractors[END_REF][START_REF] Ott | Chaos in dynamical systems[END_REF][START_REF] Abarbanel | Analysis of observed chaotic data[END_REF][START_REF] Haller | Distinguished material surfaces and coherent structures in threedimensional fluid flows[END_REF]. They depend on time and on the initial conditions of the dynamical system. The FTLs λ i (t, x(0)), i = 1, . . . , d can be calculated from the numerical solution of the ordinary differential equations (ODE) for a finite time t and initial state x(0) at t = 0:

ẋ = f (x) u = df dx u λ i t, x(0) = 1 2t log Λ i (u T • u) (2)
where ẋ denotes a system of ODEs, and u a matrix differential equation for d initially orthonormal perturbation vectors which are the columns of the matrix u, and u T being the transpose of u and Λ i (u T • u) denoting the eigenvalues of u T • u.

Alternatively, λ i can be calculated directly from differences between trajectories and initially perturbed trajectories [START_REF] Aldridge | Direct Lyapunov exponent analysis enables parametric study of transient signaling governing cell behavior[END_REF]. The number of Lyapunov exponents is equal to d and the whole set of them is also called the Lyapunov spectrum. The largest Lyapunov exponent quantifies the exponential divergence of the most unstable direction and the lower Lyapunov exponents quantify the exponential divergence in the d -1 orthogonal directions.

Comparing the definitions of MCC and FTL the differences can be summarized as follows: The FTLs identify the most unstable direction of the state space and the stability in all orthogonal directions on the basis of absolute distances between the trajectory and perturbed trajectories. The MCCs measure the response-

A c c e p t e d m a n u s c r i p t

perturbation ratio on the basis of relative changes of state variables.

Recent experimental results show that biological responses can be absolute or relative: In EGF stimulated H1299 cells the absolute change of ERK2 response varies in different cells but the relative response is the same in different cells [START_REF] Cohen-Saidon | Dynamics and Variability of ERK2 Response to EGF in individual Living Cells[END_REF]. An absolute response mechanism seems to occur in some bacterial systems [START_REF] Shinar | Input output robustness in simple bacterial signaling systems[END_REF].

In our study we have numerically calculated the MCC according to Eq. ( 1) and the FTL according to Eq. ( 2). We have coded the calculations in Matlab [28]. For the eigenvalues in Eq. ( 2) we have also coded an RQ decomposition but our results do not depend on it. The Matlab code can be obtained from the authors upon request.

Thresholds in transient dynamics

Cells respond to changes in protein concentrations, which suggests that concentration changes should be the basis for a definition of thresholds in cellular signaling.

As it has been discussed in the introduction sigmoidal stimulus-response curves have been related to activation thresholds of biochemical networks [START_REF] Goldbeter | An amplified sensitivity arising from covalent modification in biological systems[END_REF][START_REF] Ferrell | The biochemical basis of an all-or-none cell fate switch in Xenopus oocytes[END_REF][START_REF] Nash | Multisite phosphorylation of a CDK inhibitor sets a threshold for the onset of DNA replication[END_REF][START_REF] Bhalla | MAP kinase phosphatase as a locus of flexibility in a mitogen activated protein kinase signaling network[END_REF][START_REF] Huang | Ultrasensitivity in the mitogen-activated protein kinase cascade[END_REF][START_REF] Bentele | Mathematical modeling reveals threshold mechanism in CD95-induced apoptosis[END_REF][START_REF] Gunarwardena | Multisite protein phosphorylation makes a good threshold but can be a poor switch[END_REF][START_REF] Salazar | Versatile regulation of multisite protein phosphorylation by the order of phosphate processing and protein-protein interactions[END_REF]. So far, experimental and modeling efforts have focused on sigmoidal stimulus-response curves in the steady state.

We suggest that a threshold value of a sigmoidal stimulus response curve can be defined by the stimulus that corresponds to the inflection point of the response curve. An alternative definition could be based on the highest curvature of the stimulus-response curve. Due to the fact that an interval of initial conditions can have the same curvature, we can only determine a threshold interval instead of a threshold value. The two suggested definitions of a threshold contain only information about the level of a stimulus at which a threshold exists but they do not contain information about its shape. The threshold encoded by the shape of a sigmoidal stimulus-response curve can vary from very steep thresholds, approaching a switch, to smoother thresholds with a lower steepness.

A typical sigmoidal stimulus-response curve, found in many references, is shown in the first column of Fig. 1. In the second column a sigmoidal stimulus-response curve is shown, where the response increases without saturation. Another possibility is a hyperbolic stimulus-response curve, where the response increases little below the interval of the highest curvature, while above, it increases without saturation. We call this interval of the highest curvature a "hyperbolic threshold", illustrated in the third column of Fig. 1. These definitions can be applied to steady states as well as to transient dynamics.

We here consider a stimulus as either some external input, say extracellular ligands, as well as initial condition (concentrations) of intracellular signaling proteins [START_REF] Bentele | Mathematical modeling reveals threshold mechanism in CD95-induced apoptosis[END_REF]. An example for an extracellular stimulus is the gene transcription network studied in this work. The transient intracellular response can also sigmoidally depend on the initial condition (concentration) of an intracellular signaling protein as in the apoptosis decision network investigated further below. Both types of "stimuli" as thresholds have been investigated in [START_REF] Bentele | Mathematical modeling reveals threshold mechanism in CD95-induced apoptosis[END_REF]. If there is no specific interest in certain network components, one should vary the initial conditions of all variables to identify thresholds. In our two studied examples we have restricted ourselves to one or two initial conditions because our objective was to compare FTL and MCC and not a detailed analysis of the models.

In transient signaling, threshold values can depend on time, which has to be considered in the discussion about the biological relevance for a predicted threshold. The amplitude of a stimulus-response curve is also essential for this evaluation. If the amplitude is too small, then the threshold effect can be lost by noise. As raised in the introduction, in larger biochemical networks the identification of putative thresholds would become very time consuming because the network components could show sigmoidal-stimulus response curves during different time intervals. The activation of a signaling pathway at the membrane-bound receptor requires a particular threshold concentration of the stimulus. For downstream components their threshold value may be higher than further up in the network or a threshold values has to be passed to switch off an inhibitor of the pathway, as shown in the right-most column of Fig. 1 (see the related discussion in [START_REF] Legewie | Quantitative analysis of ultrasensitive response[END_REF]). Therefore quantitative measures to identify putative thresholds in transient dynamics and concepts to predict their downstream consequences are required.

Quantitative measures with an extreme value at the inflection point, or during the interval with the highest curvature, would allow a fast identification of putative

A c c e p t e d m a n u s c r i p t

thresholds in larger biochemical networks. These predictions would then have to be confirmed in experiments. In the case of different threshold values the analysis of the related downstream networks will answer the question whether the threshold to activate a downstream event is an integrated response.

Summarizing our discussion, we seek quantitative measures that identify the inflection point or the interval of the highest curvature of a stimulus-response curve. To evaluate whether FTLs or MCCs are appropriate we have chosen three examples, which are simple enough to study in detail the dynamics of the trajectories and stimulus-response curves but which at the same time are complex enough to show thresholds in their transient dynamics. In our work we focus on the identification of thresholds, while in other works quantitative measures to analyze the steepness of sigmoidal stimulus-response curves have been investigated [START_REF] Gunarwardena | Multisite protein phosphorylation makes a good threshold but can be a poor switch[END_REF][START_REF] Legewie | Quantitative analysis of ultrasensitive response[END_REF].

Apoptosis decision network

A mathematical model describing the kinetics of an apoptosis decision network has been presented in [START_REF] Aldridge | Direct Lyapunov exponent analysis enables parametric study of transient signaling governing cell behavior[END_REF]. The respective reaction network is shown in Fig. 2 and the variables are summarized in Table 1. Table 1: Assignment of names for variables.

c3* ub. c8 c3* [c8*,c3] [c8,c3*] [c3*,xiap] c3 c8* xiap
state analysis does not reveal the critical initial concentrations for a survival/death decision. The authors have identified a threshold for the survival/death decision by a maximum of the largest FTL for a critical value of xiap(0). They reported that the final decision between cell survival and controlled cell death critically depends on the temporal duration of the activity of c3 * in living cells. A short duration leads to survival, while a long duration leads to cell death. We here extend the analysis of this mathematical model by calculating the three largest Lyapunov exponents and the MCCs and by comparing these quantities with the dynamics of this network. We begin with a reduction of the original eight ODEs, Eqs. (1-8) of [START_REF] Aldridge | Direct Lyapunov exponent analysis enables parametric study of transient signaling governing cell behavior[END_REF] to six ODEs, by applying conservation equations for c8 and xiap. This leads to the following set of equations: where the conservation equations are given by

d dt c8 * = k d2 • [c8, c3 * ] + (k d3 + k d4 )[c8 * , c3] -k 3 • c3 • c8 * d dt c3 = k d3 • [c8 * , c3] -k 3 • c3 • c8 * d dt c3 * = k d4 • [c8 * , c3] + (k d1 + k d2 ) • [c8, c3 * ] + k d5 • [c3 * , xiap] -k 1 • c3 * • c8 -k 5 • c3 * • xiap d dt [c8, c3 * ] = k 1 • c3 * • c8 -(k d1 + k d2 ) • [c8, c3 * ] d dt [c8 * , c3] = k 3 • c3 • c8 * -(k d3 + k d4 ) • [c8 * , c3] d dt xiap = (k d5 + k d6 )[c3 * , xiap] -k 5 • c3 * • xiap (3)
c 8 = c8 + c8 * + [c8, c3 * ] + [c8 * , c3] xiap = xiap + [c3 * , xiap] . (4)
The parameters and the initial values of the variables, including respective units are the same as in [START_REF] Aldridge | Direct Lyapunov exponent analysis enables parametric study of transient signaling governing cell behavior[END_REF] (See Table 2).

We study the dynamics of the network components for different initial conditions of xiap, the inhibitor of c3 * , and thus a critical regulator of the survival/death decision [START_REF] Aldridge | Direct Lyapunov exponent analysis enables parametric study of transient signaling governing cell behavior[END_REF]. The trajectories of all network components for eleven different xiap(0) are shown in Fig. 3. They show a more or less pronounced clustering into two groups for each network component, thus indicating thresholds. Groups with fully overlapping trajectories are annotated by " * ". On the basis of the trajectories we have chosen time points between t = 6h and t = 18h to study transient dynamics of the system.

In Fig. 4 we compare the concentrations of the network components (first row) with the FTLs λ i (second row) and with the MCCs C(x i , x j ) (third row) as a function of xiap(0) (abscissa) for the three chosen time points (columns) during the transient signal. We call the curves in the first row of Fig. 4 stimulus-response curves, to use standard notation though the "stimulus" is the initial condition xiap(0) in our study. The proteins c8 * and c3 show a sigmoidal stimulus-response curve in agreement with the clusters of the related trajectories in Fig. 3 and thus each curve has an individual sigmoidal threshold. The intervals of the stimulus, which corresponds to the increasing part of the stimulus-response curve, are partially overlapping within a small xiap(0) interval whose position is independent of Our results for the FTLs are shown in the second row of Fig. 4, where only the three largest FTLs are shown 3 . The maximum of the largest FTL λ 1 is located at a xiap(0), which is close to the inflection points of the overlapping stimulusresponse curves for c8 * and c3. The stimulus-response curve of c8 * has the largest amplitude, while the steepness of the three curves are similar. Thus c8 * contributes most to the maximum of λ 1 , although the largest FTL has identified a sigmoidal threshold as an integrated response of c8 * , c3, and [c3, c8 * ]. The profile of λ 1 for t = 6h agrees with the result in Fig. 4 of [START_REF] Aldridge | Direct Lyapunov exponent analysis enables parametric study of transient signaling governing cell behavior[END_REF].

The second largest FTL has a maximum at a different critical xiap(0), which is found close to the thresholds for c3 * and xiap. The dynamics of the protein c3 * contribute probably only at t = 6h to the maximum of the second largest FTL because for the other two time points its stimulus-response curve is hyperbolic. The value of this critical xiap(0) decreases with time, in contrast to the case above. This means, the longer the damage persists, the higher becomes the chance for survival. This sounds unrealistic and thus it is questionable whether the xiap(0), leading to a maximum of the second largest Lyapunov exponent, is a true biological threshold. On the other hand, for the opposite case, it could be realistic that the chance for survival decreases with time.

The second largest FTL is negative and thus the respective direction in state space is stable. A maximum of a negative FTL indicates minimal stability. Related stimulus-response curves are sigmoidal and initial conditions with the point of inflection form a threshold according to our definition. These relations cannot be directly inferred from Fig. 4 due to several interacting proteins. In Appendix B we present a one-dimensional example which allows the reader to follow this fact.

Our results for the MCCs leading to the largest maxima and minima at critical xiap(0) are shown in the third row of Fig. 4. We have independently perturbed

Gene transcription network

The question arises whether there are alternative model structures whose transient dynamics are characterized by activation thresholds. Systems with cooperativity can show sigmoidal stimulus-response curves, which suggests them as promising candidates. Cooperativity is a property of many gene transcription networks [START_REF] Courey | Cooperativity in transcriptional control[END_REF]. We therefore study in this section a simple example for a gene transcription net-

ẋ1 = -x 1 ẋ2 = x 1 • (10 -x 2 ) -x 2 Model I ẋ3 = 10 • x 4 2 1 + x 4 2 -x 3 . (5)
The first equation describes the degradation of a stimulus x 1 , the second the activation of transcription factor x 2 by x 1 and the third gene transcription activated by x 2 , where we have chosen the value four for the Hill-coefficient, which is the cooperativity measure. All variables represent concentrations but without units. In contrast to the example in the previous section, the transient response is induced by a transient stimulus x 1 in this model. In addition we study two extended versions of this network by independently adding different feedback loops. In Model II the deactivation of x 2 is controlled by x 3 , while in Model III a positive feedback enhances the activation of x 2 by x 3 . The modified equations are

ẋ2 = x 1 • (10 -x 2 ) -x 2 • x 3 Model II ẋ2 = (x 1 + x 3 ) • (10 -x 2 ) -x 2 Model III . (6)
The initial condition of the stimulus x 1 (0) will be varied to identify a critical value which could form an activation threshold for the output x 3 (t). For the other initial conditions we have chosen x 2 (0) = x 3 (0) = 0. The trajectories for x 2 and x 3 of the three models for different x 1 (0) are shown in Fig. 6. The trajectories of x 3 for Models I and III show a clustering into two groups during the transients. It is difficult to detect a clustering for Model II but we will later see that the stimulus-response curve of x 3 is sigmoidal4 . We have chosen time point t = 2 to study the transient dynamics.

The related quantitative measures are presented in Fig. 7, where we compare the stimulus-response curves with the spectrum of the FTLs and with the MCCs. The stimulus-response curves are sigmoidal for x 3 in all three models and for x 2 in Model III.

The largest FTL has a maximum at the inflection point of the stimulus-response curves of x 3 in each model. It therefore identifies putative activation thresholds in the transient dynamics for these examples too, though a threshold was hard to detect in the trajectories of Model II. Note, that there is in general no one-to-one x 1 (0)

x 1 (0) x 1 (0)

x 1 (0)

x 1 (0)

A c c e p t e d m a n u s c r i p t

correspondence between the FTLs and the variables. For the model without feedback, the dynamics of x 3 is also reflected in λ 2 and λ 3 . This can easily be verified by removing the equation for x 3 from Eqs. ( 5). The resulting two-dimensional model leads to two monotonously decreasing FTLs. The maximum of the largest FTL is preserved over a large time interval during the transient signal as shown for the model without feedback in Fig. 8. Thus the threshold is a temporal robust property of the system.

For the two feedback models the maximal λ 1 is shifted to lower values of x 1 (0) as shown in the second and third column of Fig. 7. For Model II, the second Lyapunov exponent λ 2 shows a maximum at a higher critical x 1 (0) than the respective value for the maximum of λ 1 . The stimulus-response curve of x 2 has the inflection point of a sigmoidal decreasing part at this value which could be related to a stabilization of the trajectory induced by the controlled deactivation of x 2 by x 3 , acting for x 3 > 1 as a negative feedback, while for x 3 < 1 as positive feedback, according to Eq. ( 6). The other variables could also influence the profile of λ 2 . This one can also see for Model III, where only λ 1 shows a maximum and the stimulus-response curves of x 2 and x 3 are sigmoidal with their inflection points at the same x 1 (0).

For Models I and II, the curve of the MCC C(x 3 , x 1 ) is monotonously decreasing over the whole x 1 (0) interval, as shown in the last row of Fig. 7. For these two models C(x 3 , x 1 ) is not sensitive enough to identify thresholds in the transient dynamics. The positive feedback regulation (Model III) enhances the strength of the activation threshold by showing maxima in C(x 2 , x 1 ) and C(x 3 , x 1 ). For Model II, C(x 2 , x 1 ) has a negative minimal value at the inflection point of the sigmoidal decreasing part of the respective stimulus-response curve.

So far we have studied this gene transcription network for the initial conditions x 2 (0) = 0 and x 3 (0) = 0. Due to activation by other enzymes or transcription factors it could be possible that initial values of proteins or mRNA are greater than zero. The results for Model I, with two positive x 3 (0), are shown in Fig. 9. The profile of the stimulus-response curve and of the largest FTL are robust against different x 3 (0) but the profiles of the MCCs C(x 3 , x 1 ) show a maximum whose position x 1 (0) increases with increasing x 3 (0) 5 .

Our investigation of the gene transcription network has thus shown qualitative differences between FTLs and MCCs regarding the identification of thresholds in the transient dynamics. The FTLs have identified the inflection point of the sigmoidal stimulus-response curves as a threshold. The MCC C(x 3 , x 1 ) with x 3 (0) = 0 is monotonous for Models I and II. The results for the FTLs are robust

A c c e p t e d m a n u s c r i p t

for different initial x 3 (0) while the profile of the MCC shows a maximum whose position x 1 (0) depends on the initial value x 3 (0). Thus the MCCs do not reliably detect all thresholds and suffer from false positives.

The threshold of the Hill equation

Our analysis of the gene transcription network without feedback has shown that the MCC C(x 3 , x 1 ) is monotonously decreasing with increasing stimulus x 1 (0) for the default case with x 3 (0) = 0 but shows a maximum for a positive x 3 (0) whose position x 1 (0) is proportional to x 3 (0). This property motivated our analytical investigation of the MCC, the FTL and two other quantitative measures from the literature for a Hill curve with a positive response for a zero input. A positive response for zero input could have been induced by other transcription factors in case of gene transcription. The Hill equation is a standard example of a sigmoidal stimulus-response curve. With a positive initial response b, it is defined by the following equation

y = f (x) = b + x n /(1 + x n ) . (7)
Assuming that the stimulus-response curve of a one-dimensional mathematical model can be described by the Hill equation with a Hill coefficient n ≥ 2 at a time point during the transient dynamics or in the steady state. Then the FTL has its maximum at the stimulus with the inflection point of the Hill equation. Thus the position of the maximum of the FTL is independent of a positive initial response b. A right translation of the Hill curve would also change the position of the maximum of the FTL but not its value. The cooperativity index f -1 (0.9)/f -1 (0.1) is also independent of a positive initial response but its value depends on a right translation of the Hill curve [START_REF] Goldbeter | An amplified sensitivity arising from covalent modification in biological systems[END_REF].

Perturbing the stimulus x by h percent and inserting Eq. ( 7) in Eq. ( 1) results in the following equation for the MCC

C(f (x), x) = C(x) = (b + x n (1 + b)(1 + h) n )(1 + x n ) (1 + x n (1 + h) n )(b + x n (1 + b)) -1 1 h . ( 8
)
Next we calculate the stimulus x with the maximum of C(x):

d dx C(x) = n((1 + h) n -1)(b -(1 + b)(1 + h) n x 2n )x n-1 (b + x n (b(1 + h) n + 1 + b) + (1 + b)(1 + h) n x 2n ) 2 h
with the roots at

x 1 = 0 and x 2 = b (1 + b)(1 + h) n 1/2n . (9)

A c c e p t e d m a n u s c r i p t

Finally

d 2 dx 2 C = (b 2 (n -1) -bpx n -4bnqx 2n + pq(1 -2n)x 3n + q 2 (1 -n)x 4n )x n-2 (b + px n + qx 2n ) 3 h d 2 dx 2 C(x 2 ) < 0 and d 2 dx 2 C(x 1 ) = 0 with p = b(1 + h) n + b + 1 and q = (1 + b)(1 + h) n .
(10) Our result of Eq. [START_REF] Ferrell | Tripping the switch fantastic: how a protein kinase cascade can convert graded input into switch-like outputs[END_REF] shows that the MCC C(f (x), x) for the Hill curve without a positive initial response does not have any extreme value for x > 0. If the Hill curve has a positive initial response then the MCC has a maximum whose position x depends on the value of the initial response, as shown by Eqs. ( 9) and [START_REF] Buchler | Protein sequestation generates a flexible ultrasensitive response in a genetic network[END_REF]. Thus a maximum of the MCC neither detects the inflection point, nor the interval of the highest curvature hyperbolic threshold whose positions are independent of b.

Alternatively, the Hill equation can be defined with a basal activity which could correspond to a basal transcription rate or a basal rate for phosphorylation. The transcription rate or the phosphorylation rate have a maximum independent on the basal activity which is different to the first model where the maximal activation is proportional to the initial positive response. Introducing a basal activity leads to the following form of the Hill equation:

y = f (x) = b + (1 -b)x n /(1 + x n ) . (11)
For the MCC it follows that

C(f (x), x) = C(x) = (b + (1 + h) n x n )(1 + x n ) (1 + (1 + h) n x n )(b + x n ) -1 1 h d dx C(x) = (b -(1 + h) n x 2n )(b + (1 + h) n (1 -b) -1)nx n-1 (b + ((1 + h) n b + 1)x n + (1 + h) n x 2n ) 2 h
with the root leading to a maximum at

x = b (1 + h) n 1/2n . ( 12
)
The MCC has a maximum for a positive basal activity for this definition of a Hill equation too.

The Hill equation on its own is frequently used to quantify ultrasensitivity in biochemical networks by fitting the Hill equation to the stimulus-response curve of a network. However this approach is not appropriate if the response exhibits basal activation [START_REF] Legewie | Quantitative analysis of ultrasensitive response[END_REF]. An alternative approach to quantify sensitivities of networks with 7) in [START_REF] Legewie | Quantitative analysis of ultrasensitive response[END_REF])

C(f ) = n (1 -f )f (1 -b) b + f (1 -b) with f = x n 1 + x n . ( 13
)
The response as a function of the activated fraction shows a maximum at

f = 1 1 -b -b + b 2 + b(1 -b) . ( 14
)
The result is similar to the results for the MCC according to Eqs. ( 9) and ( 12):

The response coefficient of Eq. ( 13) has a maximum at a critical f whose position depends on the basal activity and the maximum vanishes for zero basal activity.

In experiments, cells are frequently starved to study their response to a specific stimulation profile, without the influence of basal activity or without the influence of a positive initial response. The reason for this is that a signaling pathway can be active without stimulation in a cell culture. Our analysis has shown that quantifying the response by the MCC, or by the response coefficient of Eq. ( 13), can lead to qualitatively different results for experiments with and without starvation. In contrast, a maximum of the FTL identifies the increasing part of a stimulusresponse curve and is thus robust against changes of basal activities or against a positive initial response, as our results for the gene transcription network of Eqs. [START_REF] Cheong | Transient IKK activity mediates NF-κB temporal dynamics in response to a wide range of TNFα doses[END_REF] have shown (Fig. 9).

Most unstable direction of relative changes

The investigated examples have shown differences of the two quantitative measures regarding the prediction of thresholds in transient dynamics, which could be related to the different properties of these two measures (Section 2): The FTLs identify the most unstable direction of the state space and the stability in all orthogonal directions on the basis of absolute distances between the trajectory and perturbed trajectories. The MCCs measure the response-perturbation ratio on the basis of relative changes of state variables.

In this section, we introduce an alternative measure which combines properties of the FTL and of the MCC: The most unstable direction of the state space is a useful property of biochemical networks with overlapping stimulus-response curves of different proteins as shown in Fig. 1. We start with the definition of the MCCs according to Eq. ( 1). The column vectors of Eq. ( 1) form an orthonormal system for t = 0. We calculate the spectral norm SN of Eq. (1) according to The different results for FTLs and MCCs are related to the different properties of these two measures. We introduced the relative instability exponent, called RIE. The RIE identifies the most unstable direction of relative changes. Our analysis led to the result that the RIE shares more properties with the MCCs than with the FTLs and is therefore not appropriate to identify thresholds.

SN(t, x 0 ) = Λ max (C T • C) , (15)
We investigated examples which are complex enough to show thresholds in transient signaling but at the same time allow us to compare in detail the dynamics of the trajectories with respect to the FTLs and MCCs. In large biochemical networks such a detailed comparison would be very time consuming. Therefore quantitative measures to identify putative thresholds in the transient dynamics are required. Our investigation suggests that they should identify the most unstable direction and that they should be defined on the basis of absolute differences between responses. The spectrum of the finite-time Lyapunov exponents unambigu-

A c c e p t e d m a n u s c r i p t

ously and reliably identifies putative thresholds as inflection points of sigmoidal stimulus-response curves in transient dynamics. Relative measures, like metabolic control coefficients, do not reliably detect all thresholds and suffer from false positives. A consequence for biochemical networks that show relative responses, is that FTLs are needed to identify thresholds, while MCCs should be applied to compare the amplitude of different thresholds.

Finally we discuss the application of FTLs and MCCs for the prediction of therapeutical interventions in biochemical networks [START_REF] Wolkenhauer | SysBioMed report: Advancing systems biology for medical application[END_REF]. For mitochondrial diseases it has been shown that the hyperbolic threshold for the inhibition of different mitochondrial complexes is a linear decreasing function of the respective flux control coefficient [32]. However, it has not been studied whether a threshold can be identified by extreme values of the MCCs.

For the administration of a drug, which binds as an extracellular ligand to a specific receptor on the cell membrane, the dose of the drug must exceed a threshold to activate the intracellular biochemical network. On the other hand, the dose of the drug should not be higher than necessary to avoid unwanted side effects. A maximum of an FTL can identify the minimal initial concentration of a drug that leads to the activation of a downstream network. FTL could also support the design of dose split experiments to avoid that the dose is decreasing below a threshold value. The threshold value could have changed over time due to changes of intracellular protein concentrations.

MCCs are applicable in medical systems biology for a comparison of the responses of different state variables in dependence on perturbations of different parameters or initial conditions of state variables [33]. Such perturbations could be induced by small molecules. In case that one is interested in a measure summarizing the responses of all variables the RIE could be useful. Global methods of sensitivity analysis need to be applied to study simultaneous perturbations of several initial conditions, [START_REF] Zhang | Comparative study of parameter sensitivity analyses of the TCR-activated Erk-MAPK signaling pathway[END_REF]. A comparison with the result for the RIE would be interesting.

FTLs and MCCs can also be applied in combination to support the design of experiments in systems biology: An analysis of MCC leads to the identification of drugs or small molecules which most effectively perturb a biochemical network. Then FTLs will determine the optimal stimulus profile of a drug or small molecule.

A c c e p t e d m a n u s c r i p t A Quantified threshold values

In the following two tables we summarize the values for xiap(0) leading to maxima of the largest and second largest FTL and to extreme values of the MCCs in Figs. 4 and5.

Row in Fig.

A c c e p t e d m a n u s c r i p t B Maximum of negative FTL

The dynamics of the following one-dimensional example ẋ = -x/(1 + x 6 )x leads to a maximum of a negative FTL . The trajectories and the FTL are shown in Fig.

Figure 1 :

 1 Figure1: Thresholds in stimulus-response curves. The possible thresholds are indicated by dashed lines. For the two sigmoidal curves the threshold can be defined by the inflection point or by the interval with the highest curvature. For the hyperbolic curve the threshold can be defined by the interval with the highest curvature. The right-most column shows a situation in a biochemical network with stimulus-response curves of different proteins.

Figure 2 :

 2 Figure 2: Apoptosis decision network. Figure adopted from [2].

  k 1 2.67•10 -9 cell • (s • molecules) -1 c8 3 .5 • 10 5 molecules/cell k d1 0.01 s -1 c8 * 5 • 10 4 or 10 5 molecules/cell k d2 0.008 s -1 c3 2 . 6 •10 5 molecules/cell k 3 6.8•10 -8 cell • (s • molecules) -1 c3 * 100 molecules/cell k d3 0.05 s -1 [c8 * , c3] 100 molecules/cell k d4 0.001 s -1 [c8 * , c3] 100 molecules/cell k 5 0.00007 cell • (s • molecules) -1 xiap 10 3 to 10 5 molecules/cell k d5 0.0000167 s -1 [c3 * , xiap] 100 molecules/cell k d6 0.000167 s -1

Figure 3 :

 3 Figure 3: Apoptosis decision network. Trajectories of the network components for xiap(0) = 0 molecules/cell, 10 4 molecules/cell, 2 • 10 4 molecules/cell . . . , 10 5 molecules/cell. An arrow indicates the direction of increasing xiap(0). A "*" indicates a cluster of overlapping trajectories.

A c c e p t e d m a n u s c r i p t time 1 .

 1 Thus their inflection points and the intervals with the highest curvature, where the response raises, are in close vicinity. The three components are directly linked in the network of Fig.2.The curve of xiap is sigmoidal too and for t = 6h also the curve of c3 * . For t = 12h and t = 18h, c3 * has not a sigmoidal but a hyperbolic threshold. The sigmoidal threshold of xiap and the hyperbolic threshold of c3 * are in close vicinity 2 . Their positions depend on time and are located around different xiap(0), compared to the sigmoidal threshold of c8 * and c3 for t = 12h and t = 18h.

  ,c3) C(c3,c3) C(xiap,c3) C(c3,c3)

Figure 4 :

 4 Figure 4: Apoptosis decision network. Quantitative measures to characterize transient dynamics. c8 * (0) = 5 • 10 4 molecules/cell, abscissa: xiap(0). First row: stimulus-response curves. Second row: FTLs and third row: Profiles of MCCs which lead to largest maxima and minima at the two critical xiap(0) are shown. Columns: different time points. Dashed lines have been added for a better comparison of the maxima of FTLs, extreme values of MCCs and thresholds in stimulusresponse curves. Same colors are chosen for the same response variable in the first and third row. The quantified threshold values are summarized in Table 3 in Appendix A.

A c c e p t e d m a n u s c r i p t

  the initial values of all cellular components by one percent and then calculated the resulting responses for each xiap(0). Comparing the results for the MCCs with the results for the FTLs and the results for the stimulus-response curves, shows that the MCCs have extreme values for the two values of xiap(0) at which the first and second largest FTL have their maximum. Perturbations in c3(0) cause most sensitive responses. The global maximum is the MCC, C(c3 * , c3), found at the critical xiap(0) with the maximum of the second largest Lyapunov exponent. At this xiap(0) the stimulus-response curve of c3 * shows a sharp concentration change from an almost linear decrease to zero, which we have defined as a hyperbolic threshold. For an initial condition of c8 * (0) = 10 5 molecules/cell, the difference between the global extremum of the MCCs and the other extreme MCCs is much more pronounced, as shown in Fig. 5. A maximum of MCCs is visible only at the smaller critical xiap(0). For this initial value, c3 * has a hyperbolic threshold for all three time points. Comparing our results for the FTLs and MCCs, in order to identify thresholds in the transient dynamics of the apoptosis decision network, shows that the maxima of the two largest FTLs are found at two critical xiap(0) where also some MCCs have extreme values. The larger value of xiap(0) is located close to the inflection points of c8 * , c3 and [c8 * , c3]. The smaller xiap(0) is located close to the inflection point of xiap and for c8 * (0) = 5 • 10 4 molecules/cell, t = 6h it is also close to the inflection point of c3. Both measures have qualitatively identified the same two critical xiap(0) as thresholds in the transient dynamics. However, there is a quantitative difference: The largest FTL has identified a critical xiap(0), where some MCCs have extreme values too but they do not form the global extremum of the MCCs. In contrast, the maximum of the second largest FTL is found at the same critical xiap(0) as the global extremum of the MCCs. This difference drastically increases with increasing c8 * (0). The value of the critical xiap(0) with the global extremum of the MCCs depends on time such that the chance for survival increases with time which appears biologically unrealistic. The value of the critical xiap(0) with the maximum of the largest FTL is independent of time and is thus more likely to be a biological threshold.

Figure 5 :

 5 Figure 5: Apoptosis decision network. Quantitative measures to characterize transient dynamics. c8 * (0) = 10 5 molecules/cell. Abscissa: xiap(0). First row: stimulus-response curve, second row: FTLs and third row: Profiles of MCCs which lead to largest extreme values are shown. Columns: different time points. Dashed lines have been added for a better comparison of the maxima of FTLs, extreme values of MCCs and thresholds in stimulus-response curves. Same colors are chosen for the same response variable in the first and third row. The quantified threshold values are summarized inTable 4 in Appendix A.

A c c e p t e d m a n u s c r i p t

  work. The temporal concentration changes of the network components are described by the following set of ODE (Model I):

Figure 6 :Figure 7 :

 67 Figure6: Gene transcription network. Trajectories. First column: without feedback regulation, second column: control of deactivation of x 2 , third column: positive feedback. The different curves result from x 1 (0) = 0.1, 0.2, . . . , 1 for the first column and from x 1 (0) = 0.025, 0.05, . . . , 0.25 for the second and third column and with x 2 (0) = 0 and x 3 (0) = 0 for all columns. An arrow indicates the direction of increasing x 1 (0). A "*" indicates a cluster of overlapping trajectories. The variables are without unit.

1 Figure 8 : 7 .

 187 Figure 8: Gene transcription network, Model I, transient response. Largest FTL has a maximum for a broad time interval. Parameters as in the left column of Fig. 7.

Figure 9 :

 9 Figure 9: Gene transcription network, Model I. Transient response at t = 2 for different x 3 (0) in dependence on x 1 (0) (abscissa). First column: x 3 (0) = 0.01, second column: x 3 (0) = 1.

A c c e p t e d m a n u s c r i p t

  basal activity has been introduced in[START_REF] Legewie | Quantitative analysis of ultrasensitive response[END_REF]. The authors have defined the MCC as a function of the activated fraction f (Eq. (

A c c e p t e d m a n u s c r i p t

  stimulus-response curves are considered as typical thresholds in biochemical networks. We examined a mathematical model of an apoptosis decision network (mutual activation cycle) and of a gene transcription network (cooperative regulation) as well as the Hill equation that can show thresholds in their transient dynamics. We then compared finite-time Lyapunov exponents (FTL) with transient metabolic control coefficients (MCC) as quantitative measures to predict thresholds in the dynamics of these networks.We obtained the following results: A maximum of an FTL identifies the value of an external stimulus or an initial condition, leading to the inflection point of a sigmoidal stimulus-response curve. In case of overlapping sigmoidal stimulusresponse curves a FTL identifies a stimulus value close to the different inflection points of the responses. The largest FTL can identify a different threshold value than the second largest FTL, which are related to stimulus-response curves with nonoverlapping rising parts.For the gene transcription network and for the Hill equation our results demonstrate that at the stimulus value, where the respective stimulus-response curve has an inflection point, in several cases the MCCs have no extreme values. Furthermore, some MCCs have extreme values at stimuli where the respective stimulusresponse curve has no point of inflection and in two cases they show extreme values near the interval with the highest curvature.The FTLs are robust against variations of initial values of network components which are not the stimulus of the relevant stimulus-response curves. The profiles of the MCCs can depend on the initial values of other components, as shown for the gene transcription network and the Hill equation. These results are important for experimental studies of signalling pathways under starvation because mathematical models, based on data from experiments either with or without starvation, could show different sensitivities.

12

 12

Figure 10 :

 10 Figure 10: Relative instability exponent (RIE) for the apoptosis decision network.

Figure 11 :

 11 Figure 11: Relative instability exponent (RIE) for the gene transcription network without feedback at t = 2.

Figure 12 :

 12 Figure 12: Maximum of a negative FTL.

Table 2 :

 2 Apoptosis model: Parameter values and initial conditions.

Table 3 :

 3 Apoptosis decision network: xiap(0) in 10 4 molecules/cell leading to maximal FTLs and extreme values of MCCs in Fig. 4.

	4	component	6h 12h 18h
	inflection point	c8 *		45	44	44
		c3		46	46	46
	FTL	λ 1		45	44	44
	MCC	C(c8 * , c3), C(c3, c3)	46	45	45
	inflection point	c3 *		45
		xiap		46	34	23
	highest curvature c3 *			33	33
	FTL	λ 2		45	32	25
	MCC	C(xiap, c3), C(c3 * , c3) 46	33	22
	Row in Fig. 5		component	6h 8h 10h
	inflection point	c8 *	60 60	60
			c3	66 67	68
	FTL		λ 1	61 62	61
	MCC
	inflection point	c3 *	57
			xiap	60 48	39
	highest curvature c3 *	45	37
	FTL		λ 2	55 45	38
	MCC		C(c3 * , c3)	58	38
			C([c8, c3 * ], c3)	47

Table 4 :

 4 Apoptosis decision network: xiap(0) in 10 4 molecules/cell leading to maximal FTLs and extreme values of MCC in Fig. 5.

The sigmoidal curve of [c8 * , c3] has its inflection point within this interval too but it has a too small response amplitude and is thus not visible in Figure4.

The curve of [c8, c3 * ] has a maximum at the same xiap(0) but it is too small and is thus not visible in this figure.

The three smallest FTLs could not be determined due to numerical problems with orthonormalization of the respective tangent vectors. We have used the RQ decomposition for orthonormalization.

The clustering is more pronounced for a higher Hill-coefficient. In order to study whether quantitative measures can predict thresholds that are difficult to see by visual inspection, we preferred a moderate value, which is still greater than a typical case of a dimer. To see whether this holds for all cases requires confirmation by experiments.

In Fig.9, merging the two figures of each row into one figure, then no visible differences can be seen for x 3 and λ 1 . The maximum of the two curves for C(x 3 , x 1 ) are found at different positions, which are close to the interval with the highest curvature of the stimulus-response curve for x 3 .

Acknowledgements K.R. acknowledges support by the German Ministry for Education and Research (BMBF) through the FORSYS junior research group programme (grant number 0315255). O.W. received support through the Helmholtz Alliance on Systems Biology and acknowledges support by the Stellenbosch Institute for Advanced Study (STIAS), South Africa. We thank M. Harle, T. Millat, M. Swat, N. Blüthgen, J. Snoep, F. Winter and R. Klages for useful discussions.

A c c e p t e d m a n u s c r i p t

where Λ max (C T • C) denotes the largest eigenvalue of C T • C. The eigenvector associated with the largest eigenvalue has the maximal stretching. Thus the SN identifies the most unstable direction of the relative changes ∂ log x i (t)/∂ log x j (0). The value of the SN depends on the size of the percentaged perturbation of the x j (0). For perturbations of one percent or less the differences are negligible in our studied examples.

In analogy to the largest FTL we define the relative instability exponent RIE by

The motivation for this is as follows. While the spectral norm helps us to identify the maximal stretching, taking the logarithm and dividing by time ensures that the profiles obtained are comparable to the largest FTL.

The following example illustrates the differences between FTL and RIE: ẋ = cx. We obtain a FTL equal to c and RIE = log(1)/t = 0. This is due to the fact that the system variable and the perturbation evolve according to the same equation, which leads to a value for RIE equal to zero. Fig. 10 shows the RIE of the apoptosis decision network for the same parameters which have been used in Figs. 4 and5. The RIE shows maxima at the same critical xiap(0), where the FTLs and MCCs show extreme values too. For the two later time points the RIE at the sigmoidal threshold of xiap is greater than the RIE at the sigmoidal threshold of c8 * in three of four cases. Thus for this example the profile of the RIE is more similar to the profile of the MCCs than to the profile of the FTLs.

The RIE for the gene transcription network without feedback is shown in Fig. 11 for three different x 3 (0). The RIE sensitively depends on this initial condition. Its curve has a similar profile as the curves of the MCC C(x 3 , x 1 ) in Figs. 8 and9.

Our two examples have shown that the dependence of the RIE on the stimulus is more influenced by the relative changes than by the most unstable direction. From our discussion and comparison of quantitative measures it follows that the RIE is less appropriate to identify thresholds in transient signaling, suggesting that quantitative measures introduced to identify thresholds in transient signaling should be defined on the basis of absolute differences between trajectories and perturbed trajectories.

Summary and conclusions

We investigated transient dynamics of mathematical models describing biochemical networks and compared quantitative measures to identify thresholds in transient responses. Thresholds separate different behaviors of a system. Sigmoidal


Ahmed Ahmed

Gabriel Frey

Christian J Michel

Essential molecular functions associated with the circular code evolution

Keywords: circular code, evolution, evolutionary trinucleotides, circular code stability, trinucleotide stability, comparative genomics, database, molecular functions, essential genes

A circular code is a set of trinucleotides allowing the reading frames in genes to be retrieved locally, i.e. anywhere in genes and in particular without start codons, and automatically with a window of few nucleotides. In 1996, a common circular code, called , was identified in large populations of eukaryotic and prokaryotic genes.

Hence, it is believed to be an ancestral structural property of genes. A new computational approach based on comparative genomics is developed to identify essential molecular functions associated with circular codes. It is based on a quantitative and sensitive statistical method (FPTF) to identify 3 permuted trinucleotide sets in the 3 frames of genes, a flower automaton algorithm to determine if a trinucleotide set is a circular code or not, and an integrated Gene Ontology and Taxonomy (iGOT) database. By carrying out automatic circular code analyses on a huge number of gene populations where each population is associated with a particular molecular function, it identifies 266 gene populations having circular codes close to . Surprisingly, their molecular functions include 98% of those covered by the essential genes of the DEG database (Database of Essential Genes). Furthermore, 3 trinucleotides GTG, AAG and GCG, replacing 3 trinucleotides of the code and called "evolutionary" trinucleotides, significantly occur in these 266 gene populations. Finally, a new method developed to analyse and quantify the stability of a set of trinucleotides demonstrates that these evolutionary trinucleotides are associated with a significant increase of the stability of the common circular code . Indeed, its stability increases from the 1502th rank to the 16th rank after the replacement of the 3 evolutionary trinucleotides among 9920 possible trinucleotides replacement sets.

Introduction 1.Minimal gene set

The study of essential genes and the analysis of their features are obviously of great interest in basic and applied researches. We recall a common definition of a minimal gene set. A minimal gene set, or essential genes, is a smallest possible group of genes that would be sufficient to sustain a functioning cellular life form under the most favorable conditions imaginable, that is, in the presence of a full complement of A c c e p t e d m a n u s c r i p t 3 essential nutrients and in the absence of environmental stress [START_REF] Cho | Policy forum: genetics-ethical considerations in synthesizing a minimal genome[END_REF][START_REF] Hutchison | Global transposon mutagenesis and a minimal Mycoplasma genome[END_REF][START_REF] Mushegian | The minimal genome concept[END_REF][START_REF] Nelson | Evidence for lateral gene transfer between archaea and bacteria from genome sequence of Thermotoga maritima[END_REF]. The upper bound of a minimal gene set is the number of genes of the smallest known genome [START_REF] Koonin | Comparative genomics, minimal gene-sets and the last universal common ancestor[END_REF].

The sequencing of the Mycoplasma genitalium which is the smallest known genome until now with a size of 580 kb and 470 predicted genes, allowed an upper bound to be determined [START_REF] Fraser | The minimal gene complement of Mycoplasma genitalium[END_REF]. This bacterium is capable of living independent of its host.

The analysis of a minimal gene set aims to identify those genes that are essential to support the cell life using comparative and experimental methods. However, a minimal gene set depends on environmental conditions of the cell [START_REF] Koonin | How many genes can make a cell: the minimal-gene-set concept[END_REF][START_REF] Gerdes | Essential genes on metabolic maps[END_REF]. Some researchers defined a set of functions that should be covered by a minimal gene set to maintain the cell integrity, including translation, transcription, replication, membrane-transport and energy conversion [START_REF] Alberts | Molecular Biology of the Cell[END_REF][START_REF] Gil | Determination of the core of a minimal bacterial gene set[END_REF].

A minimal gene set is an essential factor for further experimental and theoretical researches as the full-synthesize and semi-synthesize of a functional cell [START_REF] Luisi | Approaches to semi-synthetic minimal cells[END_REF][START_REF] Forster | Towards synthesis of a minimal cell[END_REF][START_REF] Gibson | Complete chemical synthesis, assembly, and cloning of a Mycoplasma Genitalium genome[END_REF] or the reconstruction of the last universal common ancestor [START_REF] Lazcano | The molecular search for the last common ancestor[END_REF][START_REF] Koonin | Comparative genomics, minimal gene-sets and the last universal common ancestor[END_REF]. A specialized DEG database (Database of Essential Genes) is developed to gather all published essential genes. It organizes the genes according to its kingdom and allows search using gene criteria or BLAST (Zhang et al., 2004;Zhang and Lin, 2009).

Comparative genomic approach

The comparative genomic is a bioinformatics approach to identify essential genes. It consists of selecting those genes that are shared between distantly related organisms [START_REF] Mushegian | A minimal gene set for cellular life derived by comparison of complete bacterial genomes[END_REF]. Evidently, genes are not identically shared between genomes. Hence, the term "shared genes" is redefined by orthologue or homologue genes [START_REF] Fitch | Homology: a personal view on some of the problems[END_REF]. This comparative approach was used when the first 2 bacterial genomes sequences of H. influenzae and M. genitalium were completed [START_REF] Fraser | The minimal gene complement of Mycoplasma genitalium[END_REF]. This comparison revealed 241 direct orthologue genes but although this gene set cannot maintain a viable cell. For certain substantial essential molecular functions, different organisms do not use genes that are orthologue or even homologue. Such genes are not identified by orthologue comparative approach, rather they are detected by examining non-orthologue genes for missing A c c e p t e d m a n u s c r i p t 4 essential functions. Such genes between genomes are called non-orthologue gene displacement NOGD [START_REF] Koonin | Non-orthologous gene displacement[END_REF]. By adding these NODG genes of H.

influenzae and M. genitalium, this minimal gene set reaches 256 genes.

Although only 15 NODG genes were identified in this bacterial study, large genome comparisons shown that the NODG genes are associated with most essential genes including transcription, translation and replication [START_REF] Koonin | Sequence -Evolution -Function[END_REF]. These genetic findings show that the concept of essential gene functions is more appropriate than that of essential genes.

Experimental approach

The experimental approach to identify essential genes existed before the comparative genomic one. It is simply based on gene-knockouts to determine lethal genes. The first experimental study generated 79 random gene-knockouts in the bacterial genome B. subtilis and identified 73 genes which did not kill the cell and 6 lethal genes [START_REF] Itaya | An estimation of minimal genome size required for life[END_REF]. It did not list the essential genes but only their proportion regarding to the total number of genes. This proportion (6 essential genes out of 79 genes, i.e. about 8%) is very close to that obtained by comparative genomic approach (256 essential genes out of 4100 genes of the studied organism, i.e. about 6%).

Genetic methods used by this approach include transposon-insertion mutagenesis (Judson and [START_REF] Mekalanos | Transposon-based approaches to identify essential bacterial genes[END_REF], plasmid-insertion mutagenesis (Vagner et al., 1998) and gene inactivation using antisense RNAs [START_REF] Woodnutt | Identification of essential genes in Staphylococcus aureus using inducible antisense RNA[END_REF].

The study of full wide genomes using experimental approaches showed many limitations. A first problem appeared because about half of gene-knockouts tend to be interrupted during experiences, and hence, a full list of essential genes is difficult to obtain. A second problem is the difficulty to identify essential cell functions performed by multiple proteins. Indeed, the disruption of one gene may not be lethal as gene redundancy exists, even if the function is essential to the cell. Another limitation of this experimental approach is that it does not consider the notion of gene evolution at all, in contrast to the comparative approach which is based on common orthologue genes.

The common circular code which was identified in large populations of prokaryotic and eukaryotic genes, is considered as an ancestral structural property of genes, and hence, as an essential property. It satisfies the principle of the comparative genomic approach. In this paper, we study the essential molecular functions using a

A c c e p t e d m a n u s c r i p t

5 comparative genomic approach based on this common circular code and its evolution. Hence, essential molecular functions will be analysed according to this circular code property.

In the next section, we recall the definition of the common circular code identified in genes of eukaryotes and prokaryotes and its main properties which will be used in the different methods developed.

The common circular code

Identification

In 1996, a simple occurrence study of the 64 trinucleotides AAA,…,TTT in the 3 frames of genes showed that the trinucleotides are not uniformly distributed in these 3 frames. By excluding the 4 trinucleotides with identical nucleotides AAA,CCC,GGG,TTT and by assigning each trinucleotide to a preferential frame (frame of its highest occurrence frequency), the same 3 subsets , and of 20 trinucleotides are found in the frames 0, 1 and 2, respectively, of 2 large and different gene populations (protein coding regions) of eukaryotes (26,757 sequences, 11,397,678 trinucleotides) and prokaryotes (13,686 sequences, 4,709,758 trinucleotides) [START_REF] Arquès | A complementary circular code in the protein coding genes[END_REF] (Table 1). Note: A 2010 statistical study of archaeal genes (150 genomes, 85,804 sequences, 70,411 kb) shows that this set is only partially retrieved in these genes (AAT is replaced by ATA, CAG by GCA, CTG by GCT, GGT by GTG and TTC by CTT). By convention, the reading frame established by a start codon ATG,GTG,TTG is the frame 0, and the frames 1 and 2 are the reading frame 0 shifted by 1 and 2 nucleotides in the direction, respectively. These 3 trinucleotide subsets present several strong biomathematical properties, particularly the fact that they are circular codes.

Before to give a few basic definitions of circular codes, the principle and the biological importance of circular codes which could construct the coding sequences, are described. In present day-genes, the principle of decoding of a DNA sequence is "not complicated": by knowing the beginning of a sequence (a start codon), a complex translation apparatus (ribosome with ARNs et proteins) allows a reliable reading of nucleotides three by three, each trinucleotide being a word coding an amino acid (except the stop codon). This efficient translation process has no frameshifting (except for the particular cases of frameshift genes). In contrast, reading a DNA sequence from a random initial position in the sequence is a "more complicated" problem. Indeed, the correct reading frame must be retrieved among the 3 potential frames. Is the randomly selected nucleotide the 1st, the 2nd or the 3rd nucleotide of a codon ?

Circular codes are sets of words. Trinucleotide circular codes, i.e. circular codes which are subsets of the genetic code, have interesting synchronizing properties. Any DNA sequence composed of a concatenation of words of a circular code can be read (decoded) from any position randomly chosen in the sequence as well as automatically (no need for a translation apparatus). Indeed, circular words sufficiently long always synchronize (at least 13 nucleotides with the common circular code). In other words, the reading frame of the sequence can always be retrieved.

In addition to this random position property, circular words have polymorphic patterns. Indeed, a maximal trinucleotide circular code (see Definition 3 and Remark 3 below) leads to 20 n synchronizing (circular) words of codon length ( nucleotides). Note: it also exists short synchronizing words of codon length (not detailed).

Synchronizing words (or circular words) in present-day genes could be markers of specific DNA boxes. Their multiple localization in the DNA sequence, even if some of them were destroyed by random mutations, and their polymorphic patterns make their identification specific and quick according to a parallel process (no need for a sequential reading of the DNA sequence from a start codon). Thus, they could help the ribosome to maintain the reading frame during translation (by checking the position of the ribosome with respect to the correct frame). As there are numerous and polymorphic, they can potentially be superposed to the genetic code. Thus, genes could have (or had) 2 codes: the classical genetic code to code the amino acids and a circular code to retrieve the reading frames of genes.

Definitions and main properties of the common circular code

We recall the definitions and the main properties of this common circular code which will be involved in this paper. , the condition implies and for .

Remark 1. The set itself is a code, more precisely a uniform code. Consequently, its non empty subsets are codes. In this paper, we call them trinucleotide codes.

Definition 2. Trinucleotide circular code: A trinucleotide code is circular if for each , , , , the conditions and imply , (empty word) and for (Figure 1).

Remark 2. The set is obviously not a trinucleotide circular code.

A circular code allows the reading frames in genes to be retrieved. It is a set of words over an alphabet such that any word written on a circle (the next letter after the last letter of the word being the first letter) has a unique decomposition into words of the circular code. As an example, let the set be composed of the 6 following words:

AAT,ATG,CCT,CTA,GCC,GGC and the word , be a series of the 9 following letters: ATGGCCCTA. The word , written on a circle, can be factorized into words of according to 2 different ways: ATG,GCC,CTA and AAT,GGC,CCT, the commas showing the way of decomposition (Figure 2). Therefore, is not a circular code. In contrast, if the set obtained by replacing the word GGC of by GTC is considered, i.e. AAT,ATG,CCT,CTA,GCC, GTC , then there never exists an ambiguous word with , such as for , and then is a circular code.

The construction frame of a word generated by any concatenation of words of a circular code can be retrieved anywhere in the generated word after the reading of a certain number of nucleotides depending on the code. This series of nucleotides is called the window of the circular code. Then, the minimal window length to retrieve the construction frame is the size of the longest ambiguous word which can be read in at least 2 frames, plus one letter. Figure 3 shows an example with the word 1). Is the letter A the 1st, the 2nd or the 3rd nucleotide of a codon ? By trying the 3 possible factorizations (frames) , and of into words of , only is possible, i.e. the first letter A of is the 2nd letter of a codon. Indeed, all the codons NAG, GTA, ATT, ACC and AAN (N being any letter) belong to . The factorization is impossible as no word of starts with AG (Table 1). The factorization is also impossible as CAA is not a word of (Table 1). The chosen word AGGTAATTACCA ( without the last A) is one of the longest ambiguous words with the circular code (12 letters).

Definition 3. Maximal trinucleotide circular code:

A trinucleotide circular code is maximal if for each , is not a trinucleotide circular code.

Remark 3. Any trinucleotide circular code with 20 words is maximal. Result 1 [START_REF] Arquès | A complementary circular code in the protein coding genes[END_REF]; Table 1). The trinucleotide set coding the reading frames (frames 0) of eukaryotic and prokaryotic genes is a maximal selfcomplementary trinucleotide circular code. will be also noted and simply called common circular code.

Therefore, the common circular code and its 2 permuted circular codes and can exist in a DNA double helix simultaneously: in a given DNA strand can be paired with in the antiparallel complementary DNA (cDNA) strand, ( shifted by 1 nucleotide in the direction) in a given DNA strand can be paired with ( shifted by 2 nucleotides in the direction) in the cDNA strand and in a given DNA strand can similarly be paired with in the cDNA strand.

Furthermore, the code allows retrieval of any frame in genes, locally anywhere in the 3 frames and in particular without start codons in reading frames, and automatically with a series of a few nucleotides (13 nucleotides in each frame; Michel, 2008, property (via)).

A recent review of circular codes in genes details the research context, the history and the other properties of this circular code (rarity, largest window length, higher frequency of "misplaced" trinucleotides in shifted frames, flexibility) [START_REF] Michel | A 2006 review of circular codes in genes[END_REF].

Result 2 [START_REF] Ahmed | Plant microRNA detection using the circular code information[END_REF]. Circular codes exist not only in genes (protein coding regions) but also in plant microRNAs which inherit the same structure of their target genes.

Result 3 [START_REF] Ahmed | Frameshift signals in genes associated with the circular code[END_REF]. Circular codes are absent in the frameshift sites of genes where the ribosome shifts the reading frame to a shifted frame 1 or 2.

Method

We present here a new computational approach based on comparative genomics to identify essential molecular functions associated with circular codes. It will also allow an evolution study of the common circular code. It is based on a quantitative and sensitive statistical method (FPTF) to identify 3 permuted trinucleotide sets in the 3 frames of genes, a flower automaton algorithm to determine if a trinucleotide set is a circular code or not, and an integrated Gene Ontology and Taxonomy (iGOT)

database. Finally, a new method is developed to analyse and quantify the stability of a set of trinucleotides. It is then applied to the circular code .

Recall of the main steps of the FPTF (Frame Permuted Trinucleotide Frequency) method

The FPTF method published earlier is a quantitative and sensitive statistical analysis of the occurrence frequencies of the 3 permuted trinucleotides in their 3 frames which is developed in order to have an automatic approach for the trinucleotide assignment to a frame [START_REF] Frey | Identification of circular codes in bacterial genomes and their use in a factorization method for retrieving the reading frames of genes[END_REF].

The set having the highest value with the function among these 3 sets is defined by such that (5)

The FPTF method allows the identification of 20 preferential sets of 3 trinucleotides in a gene population such that in each set , 3 permuted trinucleotides are assigned to 3 different frames. Therefore, 3 sets , and of 20 trinucleotides can be associated with the frames 0, 1 and 2 of a gene population , respectively.

Application example of the FPTF method

As an application example to comment on the FPTF method, we consider the population of drug binding genes. Drug binding in the gene ontology terms is a molecular function interacting selectively with a drug, any naturally occurring or synthetic substance. This population contains 217 genes from different species. 3 shows the function (Formula 3) of 9 trinucleotides , , in the group associated with the set AAC,ACA,CAA of genes in order to evaluate a trinucleotide simultaneously compared to its 2 other frames and its 2 other permuted trinucleotides. For the rest of the application example, the computations will be shown only for this set , the other sets being computed in the same manner. A trinucleotide with the highest value occurs preferentially in the frame . Table 4 computes the function , (Formula 4), evaluating the set AAC,ACA,CAA of genes . As has the highest value (among the 3 values), then AAC 0 ,ACA 1 ,CAA 2 . Therefore, the preferential frames for the 3 trinucleotides AAC, ACA and CAA are the frames 0, 1 and 2, respectively. In the same way, the preferential frames for the other permuted trinucleotide sets are computed. Table 5 shows the 3 preferential sets of 20 trinucleotides , and in frames 0, 1 and 2, respectively, of genes .

Flower automaton algorithm

A flower automaton algorithm is developed to determine if a set of trinucleotides is a circular code or not (not detailed; for the mathematical concept see [START_REF] Berstel | Theory of Codes[END_REF].

The iGOT database (integrated Gene Ontology and Taxonomy)

An integrated Gene Ontology and Taxonomy (iGOT) database is developed and associated with the FPTF (Frame Permuted Trinucleotide Frequency) method (see above) to analyse massive gene populations in order to identify automatically circular codes close to the common circular code of eukaryotic and prokaryotic genes (Figure 4). It is freely available upon request. The iGOT database is based on the Gene Ontology (GO) database (The Gene Ontology Consortium, 2000Consortium, , 2008)).

The GO database describes genes and gene products in a species independent manner. It covers the molecular function of gene products, their roles in multi-step biological processes and their localization in cellular components. For example, the genes responsible for the drug binding are found under the binding function which The GO database focuses on the gene product itself rather than the nucleic sequence of genes, thus the genes are presented as a sequence of amino acids.

Therefore, we applied some modifications on the table structures and imported the nucleic sequences of their genes. Figure 4 shows the main entities of the iGOT database which inherits some entities of the GO database. iGOT has been updated

with the last GO database release.

In the GO role hierarchy, the first 4 levels cover 1183 populations of genes while the first 5 levels cover 5566 populations of genes. Due to the exponential increase in the data-processing computing time, the number of populations in our statistical analysis has been limited to the first 4 levels, i.e. 1183 gene populations representing in total about 300,000 genes (312,000 kb).

The TERM and TERM2TERM tables store the molecular functions as well as their hierarchical organizations. The GO terms are organized in a network like structure, i.e.

a term may be assigned to many parent term groups. The GENE_PRODUCT table stores the gene entity itself rather than its sequence which is stored in the SEQ table

for the amino acid sequences and in the SEQ_EMBL table for the nucleic sequences.

The SPECIES table stores the species and their organization. The DBXREF table stores all the external references. An extension to store the FPTF calculation results is developed (not shown in Figure 4).

Circular code stability method

Principle

We present here a new method to analyse and quantify the stability of a set of trinucleotides. This general concept applied to a circular code is related to the notion of ambiguity of circular words. Short sequences of circular words (necessarily shorter than the code window length of 13 nucleotides) can occur in non-reading frames (shifted frames) of words generated by a concatenation of circular words (e.g. the words GGT, AAT and TAC in the factorization of the example in Figure 3).

The stability function of a code will measure the proportion of words of which occur in shifted frames when concatenating words of . Then, circular codes . This class of codes with stronger synchronization properties are called comma-free codes [START_REF] Berstel | Codes and Automata[END_REF]. More generally, the more the value of a code is raised, the more the probability to retrieve the reading frame increases with shorter words of . We apply here this concept to the common circular code .

Definition

As the code occurs in the reading frame (frame 0) of genes, we denote by . This property exists with the circular codes but not with the comma-free codes (see Definitions 4 and 5 in [START_REF] Michel | A relation between trinucleotide comma-free codes and trinucleotide circular codes[END_REF]. For example, the concatenation of the trinucleotides TAC and CTC (Table 1), i.e. TACCTC, leads to the trinucleotide ACC (Table 1) which thus occurs in frame 0 but also in frame 1.

Our objective is the study of the stability of trinucleotides of in the reading frame after the mixing of all possible trinucleotides of (with equiprobability) or in other words, the tendency of trinucleotides of to occur in other frames than the reading frame. Let the indicator function be equal to 1 if a trinucleotide is The more the value is raised, the more the trinucleotide in the frame is stable.

Then, the stability function (in %) of the common circular code in a frame is naturally measured by the average values of ( 8)

Similarly, the more the value is raised, the more the common circular code in the frame is stable.

Results

Essential molecular functions associated with circular codes

The first 4 levels of molecular functions of the iGOT database are scanned and 1183 gene populations sharing the same molecular function are thus constructed. Then, the FPTF method and the flower automaton algorithm are applied on these gene populations in order to identify the gene populations which are maximal trinucleotide circular codes (Definitions 2, 3 and 8) close to circular code . These maximal trinucleotide circular codes close to are noted .

The common circular code is considered as a common ancestral property of genes, and hence, as an essential structural property of genes. However, in small and particular gene populations, such as sharing the same molecular function, the code and the codes derived by evolution of are considered. This code Furthermore, a statistical study is being performed to evaluate if the circular code property is affected by gene GC content. Thus, the identified circular codes are analysed according to the GC ratio of their 266 gene populations . Figure 6 shows that they are distributed on both sides of the random GC ratio 0.5, i.e. with values greater or less than 0.25. Therefore, according to this statistical study, there is nocorrelation between the circular code property and the gene GC content.

Evolutionary trinucleotides in the common circular code

For most of the 266 gene populations , the evolutionary trinucleotides are mainly the 3 trinucleotides GTG, AAG and GCG among the 40 possible trinucleotides occurring in frames 1 or 2. Belonging to frame 1 with (belonging to , Table 1), they occur in frame 0 (reading frame) with 81.3%, 40.0% and 24.6% of these 266 gene populations , respectively (Table 6). Remark: This statistical result explain the previous choice to set the word distance (see above).

Very surprisingly, these evolutionary trinucleotides may be explained by the stability of the common circular code . Table 7 shows the stability function (Formula 8;

in %) of the 20 trinucleotides of the common circular code in the 3 frames ordered by increasing values of stability.

Result 4. The most stable trinucleotides of are CAG, CTG, CTC and GAG with a value , i.e. with no occurrence in the shifted frames (Table 7). The least stable trinucleotides of are ACC and GGT with a value (Table 7).

Remark 5. and whatever (Result 1 and consequence of Definition 9).

Result 5. The stability of the common circular code has a value (Table 7).

Remark 6.

(Result 1 and consequence of Definition 9).

In order to determine if a relation between the evolutionary trinucleotides GTG, AAG and GCG and the stability of the common circular code exists, we develop an iterative procedure which replaces a trinucleotide of by another trinucleotide of the same permuted set , , (defined in Section 2.1) and computes the stability of this new set (Formula 8). For example, if the least stable trinucleotide GGT of which tends to occur in frame 2 (Table 7) is replaced by GTG in frame 0 (as GTG GGT) then the stability computation of this new set leads to a value which is unexpectedly greater than (Result 5).

A c c e p t e d m a n u s c r i p t 19

Therefore, this set is more stable than . As there are (mainly) 3 evolutionary trinucleotides, the procedure analyses all the trinucleotides sets until a maximum of 3 trinucleotide replacements in the common circular code . There are such sets . Figure 7 shows these 9920 sets sorted according to their decreasing values of stability. This function has a maximal value equals to 89.3%, a minimal value equals to 67.8% and an average value equal to 78.6 at rank 4768 (Figure 7). The common circular code occurs at the 1502th rank (Figure 7).

Notation 2. The maximal circular code is now defined precisely as being the common circular code (maximal self-complementary trinucleotide circular code) in which the 3 evolutionary trinucleotides GTG, AAG and GCG replace the 3 trinucleotides GGT, GAA and GGC, respectively (Table 8). Note: is not selfcomplementary.

Result 6. The stability of the circular code has a value and occurs at the 16th rank among the 9920 sets, i.e. among the first 0.15% sets (Figure 7).

This surprising result agrees with a relation between the evolutionary trinucleotides and the stability of the common circular code .

A comparative statistical study of the common circular code between genes coding membrane proteins (GO term 0016020: membrane, 46,182 sequences, 71,000 kb) and genes coding soluble proteins (GO term 0005625: soluble fraction, 567 sequences, 983 kb) is also performed. In genes of membrane proteins, the common circular code has 4 evolutionary trinucleotides (CAG is replaced by GCA, CTG by GCT, GGC by GCG and GGT by GTG). In contrast, it is more conserved in genes of soluble proteins with 2 evolutionary trinucleotides (GAA is replaced by AAG and GGT by GTG). Both classes of genes have the main evolutionary trinucleotide GTG (Table 6) and each class contains one to the 2 other identified evolutionary trinucleotide AAG or GCG (Table 6).

Conclusion

In 1996, a common circular code was identified in large populations of eukaryotic and prokaryotic genes (Result 1 and Table 1). It allows the reading frames in genes to be retrieved locally, i.e. anywhere in genes and in particular without start codons, and automatically with a window of a few nucleotides. It is believed to be an ancestral structural property of genes [START_REF] Arquès | A complementary circular code in the protein coding genes[END_REF] and their subsequent papers).

We developed here a new computational approach based on comparative genomics to identify molecular functions associated with this common circular code and its evolution in close circular codes (Definition 8). In order to carry out automatic circular code analyses on a huge number of gene populations, it is based on a quantitative and sensitive statistical method (FPTF) which allows 3 permuted trinucleotide sets in the 3 frames of genes to be identified (Section 2.1), a flower automaton algorithm to determine if a trinucleotide set is a circular code or not (Section 2.3), and an integrated Gene Ontology and Taxonomy database (iGOT) (Section 2.4 and Figure 4). It identifies 266 molecular functions (Figure 5) which surprisingly include 98% of the functions covered by the essential genes of the DEG database. This observation suggests that the common circular code and its evolution in close circular codes are associated with essential molecular functions.

This approach also identifies 3 (main) evolutionary trinucleotides in the common circular code : GTG, AAG and GCG (Table 6). Very surprisingly, these 3 evolutionary trinucleotides of may be explained by the stability of . Indeed, the trinucleotide GGT GTG of (31.0% in frame 2) with a value GGT is the least stable trinucleotide of (Table 7). Similarly, the trinucleotide GAA AAG of (23.1% in frame 2, i.e. the second highest value) is also one of the most unstable trinucleotides of (Table 7). Therefore, a new circular code is identified in gene populations.

is the common circular code in which the 3 evolutionary trinucleotides GTG, AAG and GCG replace the 3 trinucleotides GGT, GAA and GGC, respectively (Table 8). A new method is developed to analyse and quantify the stability of a circular code, precisely to study the tendency of trinucleotides of to and to be compared. The common circular code has a stability value of (Table 7) and occurs at the 1502th rank among the 9920 trinucleotide sets (Figure 7).

Surprisingly, the circular code has a significant stability increase with a value of associated with the 16th rank among these 9920 trinucleotide sets (Figure 7).

This stability increase of could be related to a better reading of gene frames from a point of view of accuracy and speed. anywhere in a sequence generated by a concatenation of words of the common circular code (Table 1). 7) and occurs at the 1502th rank (right square) among these 9920 sets. The circular code with the 3 evolutionary trinucleotides GTG, AAG and GCG has a stability equal to and occurs at the 16th rank (left square) among these 9920 sets.

  The letters (or nucleotides or bases) define the genetic alphabet A,C,G,T . The set of non-empty words (words resp.) over is denoted by ( resp.). Let be the concatenation of the words for .Definition 1. Code: A set of words is a code if for each ,

  AGGTAATTACCAA... located anywhere in a sequence generated by a concatenation of words of the common circular code (Table

  the 2 occurrence probabilities of its 2 permuted trinucleotides and in the same frame . Let be the observed occurrence probability of a trinucleotide in a frame of a population of genes.Then, in a group , the function of a a group , the function of a trinucleotide computes the average probability of the 3 permuted trinucleotides , and in the frame as follows(2)A trinucleotide occurring with the highest (or lowest) probability in a frame compared to the 2 other frames can have a probability lower (or higher) than the probabilities of its 2 permuted trinucleotides in this frame . In order to evaluate a trinucleotide simultaneously compared to its 2 other frames and its 2 other permuted trinucleotides, the function of a trinucleotide is defined as the mean of its weight simultaneously in its frame and in its permutation set. Therefore, a trinucleotide with the highest value occurs preferentially in the frame , i.e. does not occur preferentially in the 2 other frames and , and the 2 other permuted trinucleotides and do not occur preferentially in the frame .

  one relation between a trinucleotide and its frame allows the 2 others relations between the permuted trinucleotides and their frames to be deduced by permutation. In order to quantify a set , , the statistical function is defined as being the mean of the function with the 3 words , and

  functions and which is localized under the molecular functions. In this paper, we excluded the biological processes and the cellular components which are not directly related to a circular code function.

  ambiguous for the frame retrieval (codes with few words in shifted frames) have a higher stability function . A circular code with a stability value of only contains words of in the reading frame (no words of in the shifted frames)

,

  being a nucleotide A,C,G,T , a trinucleotide of in frame 0. The cardinality of is . Let the di-trinucleotide be a concatenation of the

  figure. Three values are presented in each cell: the size, the and values. The size

  frame or to expand in shifted frames (Section 2.5). Furthermore, an iterative procedure analysing the stability of the 9920 trinucleotide sets (until a maximum of 3 trinucleotide replacements in ) allows the stability of the 2 codes

  , Dervyn E., Ehrlich S.D. (1998). A vector for systematic gene inactivation in Bacillus subtilis. Microbiol. 144, 3097-3104. Zhang R., Ou H., Zhang C. (2004). DEG, a Database of Essential Genes. Nucleic Acids Res. 32, D271-D272. Zhang R., Lin Y. (2009). DEG 5.0, a database of essential genes in both prokaryotes and eukaryotes. Nucleic Acids Res. 37, D455-D458. ATC ATT CAG CTC CTG GAA GAC GAG GAT GCC GGC GGT GTA GTC GTT TAC TTC ACA ATA CCA TCA TTA AGC TCC TGC AAG ACG AGG ATG CCG GCG GTG TAG TCG TTG ACT TCT CAA TAA CAC CAT TAT GCA CCT GCT AGA CGA GGA TGA CGC CGG TGG AGT CGT TGT CTA CTT

  set AAC,ACA,CAA for the example of drug binding genes . As has the highest value (among the 3 values), then AAC 0 ,ACA 1 ,CAA 2 . Therefore, the preferential frames for the 3 trinucleotides AAC, ACA and CAA are the frames 0, 1 and 2, respectively. ACC GAC TAC CAG GAG GTA ATC GAT ATT GCC CTC GGC GTC GCT TTC GTG GTT ACA AAG ATA CCA ACG ACT AGC AGG TAG TCA ATG TTA CCG TCC GCG TCG CTG TCT TGG TTG CAA AGA TAA CAC CGA CTA GCA GGA AGT CAT TGA TAT CGC CCT CGG CGT TGC CTT GGT TGT

  Stability value (in %) of the 20 trinucleotides of the common circular code in the 3 frames ordered by increasing values of stability and stability function (in %) of the common circular code in the 3 frames (last row). AAC AAT ACC ATC ATT CAG CTC CTG GAA GAC GAG GAT GCC GGC GGT GTA GTC GTT TAC TTC AAC AAT ACC ATC ATT CAG CTC CTG AAG GAC GAG GAT GCC GCG GTG GTA GTC GTT TAC TTC

Figure 1 .

 1 Figure 1. A graphical representation of the circular code definition. A set of words is a circular code if any word generated by a concatenation of words of and written on a circle has a unique decomposition into words of .

Figure 2 .Figure 3 .

 23 Figure 2. The set AAT,ATG,CCT,CTA,GCC,GGC is not a circular code as the word ATGGCCCTA written on a circle can be factorized into words of according to 2 different ways: ATG,GCC,CTA and AAT,GGC,CCT.

Figure 4 .

 4 Figure 4. Main entities of the iGOT database (integrated Gene Ontology and Taxonomy).

Figure 5 .Figure 6 .

 56 Figure 5. The molecular functions encoded by the 266 gene populations (maximal

Figure 7 .

 7 Figure 7. Stability function (in %) of 9920 sets associated with a maximum of 3 trinucleotide replacements in the common circular code and sorted according to their decreasing values of stability. The common circular code has a stability

Table 2

 2

shows the occurrence probabilities of 9 trinucleotides , , in

Table 1 .

 1 The common circular code identified in both eukaryotic and prokaryotic genes: ,

	and	are the preferential sets of 20 trinucleotides in frames 0, 1 and 2 of genes,
	respectively.

Table 5 .

 5 The

	3 preferential sets of 20 trinucleotides	,	and	in frames 0, 1 and 2,
	respectively, of the example of drug binding genes .
	Evolutionary trinucleotides	Occurrence
	GTG	81.3%
	AAG	40.0%
	GCG	24.6%

Table 6 .

 6 Main evolutionary trinucleotides and their occurrence percentages in the 266 gene populations .

Table 8 .

 8 The common circular code identified in both eukaryotic and prokaryotic genes

(first row which is identical to the first row of Table

1

) and the circular code with the 3 evolutionary trinucleotides GTG, AAG and GCG (underlined) (last row).

Acknowledgments

We thank a reviewer for his advices and suggestions.