Computing equilibrium concentrations for large heterodimerization networks (1110.2616v1)

Abstract: 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.

• The paper develops an iterative fixed-point algorithm to compute equilibrium concentrations in networks of reversible heterodimerization reactions.
• It reformulates the mass action kinetics equilibrium as a contraction mapping in a specially designed metric space, with convergence rigorously proven via the Banach fixed-point theorem.
• Application to mRNA hybridization in transcriptomics demonstrates rapid convergence and efficiency in handling networks with millions of species.

This paper (1110.2616) presents an efficient iterative algorithm for computing the equilibrium concentrations in large chemical reaction networks involving only reversible heterodimerization reactions of the form $A_i + A_j \rightleftharpoons A_{ij}$. Traditional numerical integration of the ordinary differential equations (ODEs) governing mass action kinetics can be computationally prohibitive for networks with a large number of species ($N$).

The core idea is to reformulate the equilibrium problem as a fixed-point iteration. Under mass action kinetics and assuming constant total concentrations $\overline{c}_i$ for each species $A_i$ (due to conservation laws in a closed system without production/degradation), the equilibrium concentrations $c_i$ of the free species $A_i$ satisfy the set of non-linear equations:

$c_i = \frac{\overline{c}_i}{1 + \sum_j K_{ij} c_j}$

where $K_{ij}$ are the equilibrium constants for the dimerization reactions $A_i + A_j \rightleftharpoons A_{ij}$. This equation defines a map $\vec{c} = \vec{T}(\vec{c})$.

The proposed algorithm is a simple iterative scheme: starting with an initial guess $\vec{c}^{(0)}$, subsequent iterates are computed as $\vec{c}^{(k+1)} = \vec{T}(\vec{c}^{(k)})$.

The paper provides a rigorous proof that this iterative scheme always converges to a unique fixed point, which corresponds to the true equilibrium concentrations. This convergence proof is based on the Banach fixed-point theorem. The key challenge in applying the theorem is finding a suitable metric space and demonstrating that the map $\vec{T}$ is a contraction mapping within that space. The authors construct a specific metric $$d(\vec{c}, \vec{c}') = \max_{1 \leq i \leq N} \frac{|c_i - c_i'|}{c_i + c_i'}$$ and prove that the map $\vec{T}$ is a contraction under this metric for species concentrations within a defined compact space $X$ (where $\varepsilon_i \leq c_i \leq \overline{c}_i$ for all $i$). The existence of a Lipschitz constant $q < 1$ is explicitly shown, guaranteeing convergence to a unique fixed point.

As a practical application, the authors apply the algorithm to a large-scale problem in computational biology: computing the equilibrium concentrations of messenger RNA (mRNA) fragments in solution due to their mutual hybridization.

• Problem Context: In experiments like DNA microarrays, mRNA fragments extracted from cells are in solution and can hybridize with each other if their sequences are partially complementary. This "solution hybridization" competes with the hybridization to probes on the microarray surface and can lead to underestimation of transcript abundance.
• Network Scale: The paper considered a network of $N \sim 3 \times 10^6$ different 48-nucleotide mRNA fragments derived from human transcriptome data.
• Implementation Details:
  • Initial concentrations $\overline{c}_i$ were based on typical mRNA expression levels from microarray experiments.
  • Equilibrium constants $K_{ij}$ were calculated from standard nearest-neighbor model parameters for RNA/RNA hybridization free energies.
  • A key computational optimization was implemented: since hybridization strength ($K_{ij}$) decays rapidly with decreasing complementarity, the sum $\sum_j K_{ij} c_j$ is dominated by interactions with only a few partners for each species. They approximated the sum by keeping only the top $\sim 10$ strongest interaction terms for each species. An efficient method was developed to quickly find potential hybridization partners by indexing sequences based on short (8-nucleotide) complementary stretches. This reduces the computational cost from $O(N^2)$ to something closer to $O(Nk)$ where $k$ is the number of dominant interactions considered per species.
• Results: The iterative algorithm was used to calculate the equilibrium concentrations $c_i^*$ of free mRNA fragments. The results showed that fragments from highly expressed transcripts were generally less affected by solution hybridization, while fragments from low-concentration transcripts, particularly those with strongly complementary partners, could be significantly depleted from the solution (their free concentration $c_i^*$ much lower than total $\overline{c}_i$). This underscores the importance of accounting for solution hybridization when interpreting microarray data and suggests avoiding highly interactive regions when designing microarray probes.

The paper also analyzes the convergence rate of the algorithm. While the Banach theorem guarantees convergence, it doesn't always provide a tight bound for practical purposes. Linear stability analysis around the fixed point shows that the convergence rate is related to the largest eigenvalue of a modified Jacobian matrix. They derive a practical upper bound for this eigenvalue $r < 1$. Empirical results showed that for the transcriptome example, most species converged rapidly within $\sim 50$ iterations, although some species exhibited slower convergence requiring a few hundred iterations. Analysis revealed that slow convergence is likely to occur in sub-networks where two or more species have similar high total concentrations and strongly interact with each other. Asymmetric strong interactions lead to faster convergence. For practical implementation, they suggest detecting such problematic sub-networks and potentially solving them analytically before running the iterative scheme on the full network to accelerate overall convergence.

The computational performance was demonstrated to be very efficient; the calculation for the large human transcriptome example ($N \sim 10^6$) converged rapidly and completed in only a few minutes on a standard PC, making it practical for large-scale analysis. The authors suggest that this iterative approach could potentially be extended to other types of chemical reaction networks beyond simple heterodimerization.