We consider a chemical reaction network governed by mass action kinetics and composed of N different species which can reversibly form heterodimers. A fast iterative algorithm is introduced to compute the equilibrium concentrations of such networks. We show that the convergence is guaranteed by the Banach fixed point theorem. As a practical example, of relevance for a quantitative analysis of microarray data, we consider a reaction network formed by N~10^6 mutually hybridizing different mRNA sequences. We show that, despite the large number of species involved, the convergence to equilibrium is very rapid for most species. The origin of slow convergence for some specific subnetworks is discussed. This provides some insights for improving the performance of the algorithm.
Systems of coupled chemical reactions involving many different species, i.e. reaction networks, have been intensively studied in chemistry, physics, mathematics, and engineering sciences (see e.g. [1][2][3][4][5][6][7][8][9][10]). If diffusion is fast enough and the number of molecules is sufficiently large, so that stochastic effects can be neglected, these systems can be described by a set of coupled first order ordinary differential equations (ODE), which govern the time evolution of the concentrations of each species. In the ODE description, the rates of production and consumption of the chemical species are given in term of mass action, Michaelis-Menten, or other cooperative-type kinetics [2]. In such systems different types of behavior are possible, as for instance relaxation to a unique stationary point, oscillations or multistability.
Usually the time evolution of the system can be computed through numerical integration of the ODE. However, this method can become very slow for large reaction networks. In addition, the main interest is typically the long time behavior of the system, which in absence of oscillations boils down to finding the stationary (equilibrium) concentrations of each of the chemical species.
In the present work, we will describe an efficient method to find equilibrium concentrations for a class of networks which we will refer to as hetero-dimerization networks. In these networks the species associate to form dimers, which eventually break apart giving back the single species. The method is based on an iterative scheme, of which we can rigorously prove the convergence. The proof relies on the Banach fixed point theorem.
The proposed method is very efficient for large reaction networks. As an example to show that convergence is fast, even for systems with ∼ 10 6 species, we consider the hybridization of RNA strands. This example is of relevance, for instance, for a better quantitative understanding of the reactions underlying the functioning of DNA microarrays [11][12][13][14]. It shows that some sequences tend to get effectively depleted from the solution because of partial complementarity with other strands. This brings some consequences for the design and interpretation of microarray experiments. This paper is organized as follows. In Section II we introduce the iterative algorithm and prove its convergence to the fixed point, irrespective of the initial condition. In Section III we show a concrete calculation for a network composed of a large number of hybridizing mRNA strands. Section IV discusses the convergence rate of the algorithm and Section V concludes the paper.
We consider a set of different chemical species A i (i = 1, 2 . . . N ) undergoing reversible association/dissociation reactions of the type:
where K ij and K ij are the forward and reverse rates.
The ratio of the rates must satisfy the detailed balance condition
where ∆G ij is the free energy of formation of the complex A ij (the free energy difference between the bound and unbound state).
The system is considered to be well-mixed, i.e. diffusion of all species is assumed to occur on a fast time scale compared to reaction time scales. Furthermore, production and degradation are assumed to be absent. A configuration can then be characterised solely by the concentrations of all species c i and of all dimers c ij at a given time. In addition, the following conservation laws hold for every species A i :
with c i the constant total concentration of a species.
We consider mass-action kinetics, so the concentrations evolve in time according to
and we are interested in the stationary solution dc i /dt = 0. A general theorem on Mass Action Reaction Networks, the so-called Deficiency Zero Theorem [1], guarantees that the system (1) has a unique stationary point, irrespective of the values of K ij and K ij . We omit the proof of this statement, it amounts to calculating the deficiency of the reaction network, which is an integer easily obtainable from the network topology (for more details see [1]). Given a set of rates K ij and K ij and some initial concentrations, it is always possible to compute the relaxation to equilibrium by solving the ODE in Eq. ( 4) numerically, for instance through discretization in small time steps ∆t. However, this is a very costly procedure and thus very unpractical for large networks. In addition the discretization brings some errors scaling as powers of ∆t, which accumulate during the calculation. We show here that the problem of finding the stationary point can be reformulated as an iterative problem, which is much more efficient and does not involve discretization approximations.
The detailed balance condition (Eq. ( 2)) provides some freedom in choosing the forward and reverse rates. Only their ratio needs to be fixed to guarantee convergence to the stationary point. One particularly interesting choice is K ij = 1, and thus K ij = e ∆Gij/RT . Substituting these values in Eq.(
This content is AI-processed based on open access ArXiv data.
