Short relaxation times but long transient times in both simple and complex reaction networks

When relaxation towards an equilibrium or steady state is exponential at large times, one usually considers that the associated relaxation time τ, i.e. the inverse of the decay rate, is the longest characteristic time in the system. However, that need not be true, other times such as the lifetime of an infinitesimal perturbation can be much longer. In the present work, we demonstrate that this paradoxical property can arise even in quite simple systems such as a linear chain of reactions obeying mass action (MA) kinetics. By mathematical analysis of simple reaction networks, we pin-point the reason why the standard relaxation time does not provide relevant information on the potentially long transient times of typical infinitesimal perturbations. Overall, we consider four characteristic times and study their behaviour in both simple linear chains and in more complex reaction networks taken from the publicly available database ‘Biomodels’. In all these systems, whether involving MA rates, Michaelis–Menten reversible kinetics, or phenomenological laws for reaction rates, we find that the characteristic times corresponding to lifetimes of tracers and of concentration perturbations can be significantly longer than τ.


Analytical derivation of the relaxation time
First, let us derive the general solution of the differential equation governing a set of concentrations x m (t) (m = 1, ..., N ): Separating the dependence on time and position, we set x m (t) = T (t)X m . Then The left hand side is independent of m while the right hand side is independent of t, so both must be equal to the same constant which we denote by −λ. The dependence on time is immediate: T (t) = T (0) exp(−λt). The equations eq.2 for the X m correspond to finding an eigenvector of a matrix H of eigenvalue −λ. The diagonal elements of H are −(A+B) while the off diagonal parts are restricted to nearest neighbors and are again independent of m. Given this translational invariance, the eigenvectors can be taken to be of the form X m = X 0 e (γ±jω)m . Plugging this in, one has For the studied system, both A and B are real. Let us start by searching for the real eigen values. We can solve for the real and imaginary parts: A general solution of eq.2 is a combination of the two complex conjugates X + m and X − m : Whith boundary conditions X 0 = X N +1 = 0: The search for real eigen values provides us with N real eigen values which is the rank of the matrix H, thus all the eigen values are real.
A relaxation time is associated to each λ κ : In practice, one refers to the relaxation time as the longest such times, i.e., τ = .
Case of a symmetric system: For a system that is left-right symmetric (A = B), τ ≈ (N +1) 2 Aπ 2 for large N which is characteristic of diffusing systems. (Clearly, if A = B, the equations describe a diffusing particle that has no bias towards left or right.) Case of an asymmetric system: If the system is not left-right symmetric (A = B), at large N the relaxation time become independent of N (and in fact at fixed κ all τ κ become independent of N : A crossover size for N between the two regimes τ ∝ N 2 and τ = cst: If the difference between A and B is sufficiently small, the relaxation time will behave as the case A = B, that is τ will scale as N 2 until N is large enough for the regime of constant τ to set in. To understand the scale in N where this crossover behavior arises, let us note A = B + . Then for large N This last equation shows that one has a crossover between the two regimes when The transition between a diffusing regime and a regime independent on the network size N thus occurs for a characteristic size 2 Generalisation of the concept of tracer to complex reactions

Kinetic for a tracer in a linear chain
The case of the linear chain provides a ideal case since every reaction has one substrate and one product. In such a case it is relatively easy to derive the reactions for tracers. Let us consider the reaction: with the kinetic for the conventional direction from left to right: with Where a and b are the concentrations of A and B and α, β and the Ks are the parameters of the equations. Note that the rate law is a summary of the difference between the forward(v + ) and backward(v − ) fluxes. When tracer molecules are introduced downtream in a chain of reactions, some of them will move against the current although the global flux of the reactions indicate the opposite direction. The correct expression for the tracer rate is The superscript t and tot stand respectively for tracer and total. Interestingly, the eq.7 is still true at steady state, then the dynamics of the tracers depends on their fraction a t a tot or b t b tot in the network.

Kinetic for a tracer in a realistic network
In a more complex reaction such as the dynamics of the tracer is not as straight forward as for the linear chain. Now a labelled carbon can go from A to C but it can also go from A to D; its fate depends on the chemical mechanism of the reaction. To investigate the dynamics of a tracer in a whole network, one has to determine the exact mechanism of every reaction, this task is tedious and can becomes for large networks. Thus we decided to add a simplification to the dynamics, during a reaction all the carbon of the tracer substrate are mixed and shuffled with the carbons of the of the non-tracer substrate. It results that every carbon coming from a labelled A can go either to the pool of labelled C or to the pool of labelled D.
Let us exemplify this on the tracer flow from A to C and D: a SS and a t are the total concentration of A at steady state and the tracer concentration of A. The chance for A t to take part in the reaction is proportional to its fraction. The flow produces C and D in proportions that depends on the number of carbons in each molecules and in the stoichiometry of the reaction. To prevent carbon creation, the following balance has to be satisfied: where the νs are the stoichiometric coefficients and the ns the number of carbons contains in one molecule.
Finally one obtains the evolution The get the full dynamics obtained by completting the flows F B→C+D , F C→A+B and F D→A+B .

Computing the different characteristic times
Our characteristic times are defined using infinitesimal perturbations of the steady state. Such a choice has the advantage of facilitating the treatment of the time dependence of the perturbation since one can exploit the linearity of the dynamics. Denoting by δC the infinitesimal vector giving the concentration deviation from the steady state and by J the jacobian matrix associated with the linearized dynamics, one has: C is the total concentration and C SS is the steady state concentration.
We showed in the first section that the eigenvalues of J are all real and negative in the case of a homogeneous linear chain as might be expected in a dissipative system. Our matrix formulation allows us to conveniently define the "lifetime" associated with a perturbation δC(t) introduced at time t = 0 via Because of the linearity, the characteristic time T does not depend on the amplitude of the initial perturbation, | C(t = 0)|. In addition, we shall impose the initial perturbation to be localized to just one site of the metabolic network. In consequence, and without any loss of generality, we can take the initial perturbation to be positive, and then it is easy to see that all of the components of C(t) remain positive at all times. using this property, we can interchange the integration and the taking of the norm if we work with the L 1 norm: The lifetime T as specified is dependent on the position where the perturbation was introduced at t = 0. To overcome this drawback, we define the characteristic time of the system as the maximal lifetime when considering all possible positions of the initial perturbation.
Our global algorithm to compute our characteristic time is then as follows: Data: vector_param: best_vector,conc_ss: concentration at steady state, F: function returning the derivative the system Result: lifetime of the system J = linearize F(. . . ,param) near conc_ss; max_time = 0; for i=1 to i=conc_ss do if time>max_time then max_time = time; end end return max_time; Algorithm 1: Algorithm to compute the lifetime Note that that all the previous derivation is done for the lifetime of a concentration perturbation but it remains valid for the lifetime of a tracer when taking the steaty state vector 0 and replacing the perturbation δC by the concentration of tracer.