<figure><img src="image_1.png" /></figure>

HAL
open science

# Conservation Analysis and Discrete Probabilistic Approximations for Parameter Estimation of Biochemical Networks

Olivier Bouët-Willaumez, Adrien Le Coënt, Benoît Barbot, Nihal Pekergin

<figure><img src="image_2.png" /></figure>

## To cite this version:

Olivier Bouët-Willaumez, Adrien Le Coënt, Benoît Barbot, Nihal Pekergin. Conservation Analysis and Discrete Probabilistic Approximations for Parameter Estimation of Biochemical Networks. 2025. hal-05104966

HAL Id: hal-05104966

https://hal.science/hal-05104966v1

Preprint submitted on 10 Jun 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Conservation Analysis and Discrete Probabilistic
Approximations for Parameter Estimation of
Biochemical Networks

Olivier Bouët-Willaumez¹, Adrien Le Coënt¹, Benoît Barbot¹, and Nihal Pekergin¹

Université Paris-Est Créteil, LACL,
F-94010 Créteil, France

{olivier.bouet-willaumez, adrien.le-coent, benoit.barbot, nihal.pekergin}@u-pec.fr

**Abstract.** In this paper, we aim to efficiently estimate parameters of biological systems. Biological systems, and more specifically biochemi-cal networks, are usually modeled by Ordinary Differential Equations (ODEs). Instead of working directly on the ODE models, our approach relies on approximating biological systems with discrete probabilistic models, more specifically Dynamic Bayesian Networks (DBNs). The dis-crete approximation is used to efficiently estimate the parameters of the system by relying on Bayesian inference instead of classical optimiza-tion methods that fail due to the high dimensionality and number of parameters encountered in this type of systems. While more efficient than classical approaches, our method is also limited by dimensionality. To mitigate the dimensionality issue as much as possible, we tackle the problem of model reduction of biological systems, and the implications it has on our parameter estimation method. More specifically, we discuss an exact model reduction method that allows us to simplify the system without losing information, and we adapt the discrete models and pa-rameter estimation methods to this reduction method. We illustrate our methods on some concrete case studies.

**Keywords:** Parameter estimation · Markov chains · Bayesian networks · Biochemical networks · Model reduction

# 1 Introduction

Ordinary Differential Equations (ODEs) are a widely used formalism for mod-eling dynamical systems, including biological and biochemical networks. In such systems, estimating parameters from experimental data is a central challenge, with applications in drug discovery, clinical diagnostics, and system biology.

Several ODE-based approaches have been proposed for analyzing such sys-tems, each with specific strengths and limitations, as discussed in [16]. Clas-sical estimation methods rely on repeatedly simulating ODEs and minimiz-ing a cost function, but this becomes computationally prohibitive for stiff or
high-dimensional models. Probabilistic approaches, particularly those based on Markov Chain Monte Carlo (MCMC), have been explored to tackle uncertainty and improve inference robustness [47], but they also suffer from high computational costs when the parameter space is large or when model stiffness is significant.

In prior work [26], we proposed a method for parameter estimation that transforms ODE models into discrete probabilistic models, specifically, Dynamic Bayesian Networks (DBNs), to enable faster cost evaluation and avoid the need for repeated ODE simulations. This approach has been shown to outperform classical ODE-based optimization for moderately sized systems (approximatively 10 variables) [27]. However, it still faces scalability limitations for larger models (approximatively 100 variables) due to the exponential growth of the state space and the increasing complexity of computing and storing transition probabilities within the DBN framework.

To address this challenge, we explore the use of model reduction techniques, which aim to simplify the system while preserving key dynamical properties. A comprehensive overview of such methods for biochemical networks is provided in [45]. These approaches can generally be divided into two categories. The first includes **structural simplification** methods, which reformulate the model to reveal intrinsic properties such as conservation laws. In particular, conservation analysis [40,48] exploits algebraic relationships among species to reduce the system's dimensionality without altering its dynamic behavior. The second category consists of **component reduction** methods, which simplify the system by eliminating or aggregating species or reactions. This includes classical lumping techniques [13] as well as more recent exact approaches based on symbolic computation [4,21,32], species lumping in stochastic models [9], and equivalence-based reductions using backward bisimulation principles [10, 11]. Additionally, timescale separation methods [23, 41] fall under this category, treating fast and slow subsystems separately to facilitate simplification. While such methods are promising for handling large models, their integration into parameter estimation pipelines remains nontrivial and often requires careful adaptation.

In this paper, we extend our previous discrete probabilistic framework for parameter estimation [26] by integrating conservation analysis as a model reduction technique.

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Workflow of the method. Square nodes represent data structures, while rounded nodes indicate processing phases. Nodes with a gray background and double border correspond to novel contributions introduced in this paper.</figcaption></figure>
The goal is to reduce the dimensionality of the biochemical system prior to discretization, thereby improving scalability and computational efficiency. To better illustrate the overall approach, Fig. 1 summarizes the proposed parameter estimation pipeline. The procedure, including the discretization and Bayesian inference steps, will be detailed in Section 2.

Once the reduced system is obtained through conservation analysis, we propose several alternative structures of DBNs to represent its dynamics. These different structures reflect distinct interpretations of the reduced system in terms of node selection and dependency modeling, offering various trade-offs between accuracy and computational cost. We compare their performance in terms of parameter estimation quality and runtime on two biochemical case studies.

The paper is structured as follows: Section 2 introduces the problem, the case studies, and the parameter estimation framework. Section 3 presents the conservation analysis method. Section 4 details the integration of this last reduction step into the parameter estimation pipeline. Section 5 presents the numerical results. We conclude in Section 6.

# 2 Problem Setting

## 2.1 Biochemical Networks

We consider the following set of differential equations as the base model for biochemical networks:

$$ \dot{x}_i(t) = f_i(x(t), k), \quad (1) $$

the variables $x(t) = (x_1(t), \dots, x_n(t))$ are real-valued concentrations of species at time $t$. The real-valued vector $k = (k_1, \dots, k_m)$ contains parameter values, which are supposed to be constant over time, but unknown. The ODE can be written in matrix form as follows:

$$ \dot{x}(t) = N \cdot v(x(t), k) \quad (2) $$

where $v(x(t), k) = (v_1(x(t), k), \dots, v_r(x(t), k))^T$ is the rate vector that contains the rate equations of the $r$ reactions. Note that $v_1(x(t), k), \dots, v_r(x(t), k)$ are rational functions of the state variables $x_i$ and parameters $k_j$. The Stoichiometric matrix $N$ is $n \times r$ dimensional and is relating the species to the reactions they participate in [16,28]. We now present the running examples used in this paper.

**Enzyme-catalyzed Reaction Network** is a typical biochemical enzyme-catalyzed reaction with $n = 4$ variables and $m = 3$ parameters. Its continuous dynamics is given by the system of equations (3).

Let us write the rate vector $v$ with $r = 3$ reactions: $v_1(x(t), k) = k_1x_1x_2$, $v_2(x(t), k) = k_2x_3$, $v_3(x(t), k) = k_3x_3$. Then the system of equations (3) can be written in matrix form (2) with the stoichiometric matrix $N$ as defined in (4).
$$ \begin{cases} \dot{x}_1 = -k_1 x_1 x_2 + k_2 x_3 \\ \dot{x}_2 = -k_1 x_1 x_2 + (k_2 + k_3)x_3 \\ \dot{x}_3 = k_1 x_1 x_2 - (k_2 + k_3)x_3 \\ \dot{x}_4 = k_3 x_3 \end{cases} \quad (3) \quad \mathcal{N} = \begin{bmatrix} -1 & 1 & 0 \\ -1 & 1 & 1 \\ 1 & -1 & -1 \\ 0 & 0 & 1 \end{bmatrix} \quad (4) $$

The EGF-NGF signaling pathway that is considered in this paper as a case study is a biochemical network which comprises two interacting sub-networks, namely the Epidermal Growth Factor (EGF) signaling pathway and the Nerve Growth Factor (NGF) signaling pathway. Both of EGF and NGF, and their receptors were discovered by Nobel Laureates Rita Levi-Montalcini and Stanley Cohen [12] through the intense experimental research on the regulation of chemical reactions related to cell growth. The EGF-NGF signaling pathway model was extensively studied in [5, 27] to describe the effects of EGF and NGF in rat pheochromocytoma (PC12) cells. The system's dynamics are described by the ODEs (31) in Appendix, and the network structure is depicted in Fig. 5. The original model has 32 variables, 26 reactions and 48 parameters, 28 of which are known, the remaining ones needing to be learned.

## 2.2 Discrete Approximations of Biochemical Networks

To enable the dynamic Bayesian network formalizing, we assume that both the time and state space are discretized. We now present the Discrete-Time Markov Chain (DTMC) and Bayesian Network (BN) models used to construct approximate probabilistic representations of biochemical networks.

**DTMC Approximation** Let us suppose that the ODE state intervals $\mathcal{X} = [x_1^{min}, x_1^{max}] \times \dots \times [x_n^{min}, x_n^{max}]$ are divided into $L_i$ sub-intervals: $[x_i^{min}, x_i^{max}]$, $[x_i^{1}, x_i^{2})$, \dots, $[x_i^{L_i-1}, x_i^{max}]$, and that the parameter intervals $\mathcal{K} = [k_1^{min}, k_1^{max}] \times \dots \times [k_m^{min}, k_m^{max}]$ are divided into $M_j$ sub-intervals: $[k_j^{min}, k_j^{max}]$, $[k_j^{1}, k_j^{2})$, \dots, $[k_j^{M_j-1}, k_j^{max}]$. Suppose furthermore that time is discretized with time-step $\delta$ to specify an increasing sequence of time instants $(t_j)_{j=0}^T$, where $t_0 = 0$ and $\delta \cdot T = \tau$ is the observation horizon. The concentration level of each species is given by a vector of random variables:

$$ \mathbf{X}^t = (X_1^t, X_2^t, \dots, X_n^t) $$

where $X_i^t$ denotes the the sub-interval (discrete state of) the $i^{th}$ species at time $t$. We assume that the model satisfies the Markov property:

$$ p(\mathbf{X}^{t+1} | \mathbf{X}^t, \mathbf{X}^{t-1}, \dots, \mathbf{X}^0) = p(\mathbf{X}^{t+1} | \mathbf{X}^t) $$

Thus, the joint probability distribution from time step 0 to T can be computed as follows:

$$ p(\mathbf{X}^{0:T}) = p(\mathbf{X}^0)p(\mathbf{X}^1|\mathbf{X}^0)\cdots p(\mathbf{X}^T|\mathbf{X}^{T-1}) = p(\mathbf{X}^0) \prod_{t=1}^{T} p(\mathbf{X}^t|\mathbf{X}^{t-1}) $$
When the ODE is time-homogeneous (time-invariant), the transition probabilities do not depend on the time step and remain constant, i.e. $\forall t \ge 1$:

$$p(\mathbf{X}^t | \mathbf{X}^{t-1}) = p(\mathbf{X}' | \mathbf{X})$$

This model is thus a (time-homogeneous) DTMC. The transition probabilities of a DTMC are typically stored in a transition probability matrix, $P$ such that

$$P[i, j] = p(\mathbf{X}^t = j | \mathbf{X}^{t-1} = i) = p(\mathbf{X}' = j | \mathbf{X} = i)$$

where $i, j$ are any two (vector) states of the underlying stochastic process.

However, the state space of the DTMC approximations grows exponentially with the dimension of the original system, causing the transition matrix $P$ to become prohibitively large to store explicitly for high-dimensional models. To address this limitation, we employ (dynamic) Bayesian networks, which provide a compact representation of transition probabilities.

**BN Approximation of the DTMC** In essence, a Dynamic Bayesian Network (DBN) [25] allows to model conditional (in)dependencies in a probabilistic transition system. In our case, it is indeed reasonable to assume that the evolution of a given species does not depend on the entire system (at least in the short term), but rather on a subset of the species. One could for example consider that a species' dynamic depends only on those that appear on the right hand side of the ODE. While in practice this would lead to reasonable numerical results, there are some subtleties associated with these independencies hypotheses. These are beyond the scope of this paper but are discussed in detail in [26]. To summarize, we can always make some independence hypotheses that lead to significant transition probability matrix size reduction. However, the validity and impact of these assumptions depend on how the underlying probabilistic models (DTMCs or BNs) are constructed. The formal definition a (non dynamic) Bayesian network is the following.

**Definition 1.** A Bayesian network is a pair $\mathcal{B} = (\mathcal{G}, P_B)$. $\mathcal{G}$ is a Directed, Acyclic Graph (DAG) whose nodes are random variables $\mathbf{X} = X_1, \dots, X_n$. The graph structure represent the dependency relations among random variables: there is an incoming arc to node $X_i$ from node $X_j$, if and only if the value of $X_i$ depends on $X_j$. Let $Pa_{\mathbf{X}_i}^{\mathcal{G}}$ denote the set of parents of $X_i$ in $\mathcal{G}$. The distribution over $\mathbf{X} = X_1, \dots, X_n$ can be factorized as a product:

$$p(\mathbf{X}) = \prod_{i=1}^{n} p(X_i | Pa_{\mathbf{X}_i}^{\mathcal{G}})$$

$p(X_i|Pa_{\mathbf{X}_i}^{\mathcal{G}})$ is called the Conditional Probability Distribution for node $X_i$ specified as a Conditional Probability Table (CPT) at node $X_i$. $P_B$ is the set of CPTs associated with $\mathcal{G}$'s nodes.
The independence assumptions encoded by a DAG for Bayesian networks allows to specify the joint distribution with fewer parameters compared to unstructured models, especially when the DAG structure is sparse.

Bayesian networks are generalized to temporal models. Under Markovian and time-homogeneity assumptions, the probability distribution over infinite trajectories can be specified compactly by the initial state distribution and the conditional probabilities for the transition model in one time step, $p(\mathbf{X}'|\mathbf{X})$.

**Definition 2.** A 2-Time-slice Bayesian Network (2-TBN) for a process over $\mathbf{X}$ is a Bayesian network specifying the evolution of random variables at time $t+1$, $\mathbf{X}^{t+1}$, given the random variables at time $t$, $\mathbf{X}^t$. $\mathbf{X}_I \subseteq \mathbf{X}$ is a set of input variables having no parents in $\mathcal{G}$. $\mathbf{X}'$ is the set of output random variables. The transition model, for any time step $t$, is given by conditional probabilities $p(\mathbf{X}'|\mathbf{X}_I)$. Therefore, for each output random variable $X_i$, the CPT $p(X_i|\mathbf{Pa}_{\mathbf{X}_i}^{\mathcal{G}})$ is given once to specify the transition probabilities from $\mathbf{X}^t$ at time $t$ to $\mathbf{X}^{t+1}$ at time $t+1$ with $\mathbf{X}^t \subseteq \mathbf{X}_I$ for all $t$.

Dynamic Bayesian Networks, represented as 2-TBNs, are well suited to the time-homogeneous models we consider. For this reason, we adopt them as approximate representations of such ODE systems. An example of a 2-TBN graph is shown in Fig. 2 (a), illustrating the enzyme-catalyzed reaction case. Before proceeding with parameter estimation, the CPTs must be computed. Since exact computation is infeasible, we approximate them using the following procedure.

## 2.3 Computation of a BN Approximation and Parameter Estimation of ODE-based Models

To approximate a given system with a DBN model, the transition probabilities need to be estimated for the said system. To do so, we perform a number of ODE simulations, and, for each simulation, we keep count of the visited discrete states. By doing that, we can estimate the probability of transitioning from a starting state to an arrival state as the sum of the number of trajectories that have led to transition from the starting state to the arrival state, divided by the sum of the number of trajectories that have departed from the starting state. This counting procedure also needs to keep track of the variables and parameters on which a variable depends. The details of the procedure and the associated formal definitions are given in [26].

The main task of parameter estimation consists in finding the parameter values of a system for which we have some measurements. The measurements are possibly incomplete, in the sense that they are scattered within time, or that some species are not observed, but we assume that at least some initial conditions are given. Once a DBN is computed, parameter estimation is performed using the fact that efficient inference algorithms allow to compute distributions of species concentrations over time for given parameters. More precisely, for given parameter values, and initial conditions (which can be given as distributions), the distributions of the species concentration over time are obtained by "unrolling
the DBN" (inferring the distributions at a time step given the distributions at the previous time step). From this distribution, we can reconstruct some point values that can be compared to the measurements. Then, an exploration of the parameter space allows to find the parameter values that best reproduce the measurements.

Formally, given a DBN approximating a biochemical network, and a set of measurements $(t^k, x^k)_{1 \le k \le M}$ taken from an instance of this network. We assume that $\forall k \in \{1, \dots, M\}$, there exists $j \in \{1, \dots, T\}$ such that $t^k = \delta j$. For a given parameter set $\mathbf{K}$ and initial condition $\mathbf{X}^0$, species concentrations over time are computed from the posterior by unrolling the DBN:

$$p(\mathbf{X}^{0:T} | \mathbf{K}, \mathbf{X}^0) \qquad (5)$$

Recall that the concentration of species $i$, $X_i$ is a discrete value representing an interval: $X_i = l$, $1 \le l \le L_i$ (see Subsection 2.2). Let $\bar{x}_i^l$ denote the midpoint of interval $l$: $\bar{x}_i^l = \frac{(x_i^{l-1}+x_i^l)}{2}$. For each time instant $j \in \{1, \dots, T\}$, the real-valued concentration of species $i$, $\tilde{x}_i^j$, is derived through the posterior distribution as follows:

$$\tilde{x}_i^j = \sum_{l \in \{1, \dots, L_i\}} p(X_i^j = l \mid \mathbf{K}, \mathbf{X}^0) \bar{x}_i^l \qquad (6)$$

The parameters are then estimated by minimizing the sum of squared errors:

$$\varepsilon(\mathbf{K}) = \sum_{\{j | \exists k \text{ s.t. } t^k = \delta j\}} (x^k - \tilde{x}^j)^2 \qquad (7)$$

The optimization problem: argmin $\varepsilon(\mathbf{K})$ can be solved by brute-force if the search space is small enough (e.g. for the enzyme-catalyzed reaction example with moderate state and parameter discretization), or with more advanced strategies when it cannot be explored exhaustively (e.g. for the EGF-NGF example). Many hybrid optimization strategies combining global exploration and local refinement have been proposed for discrete problems lacking derivative information [6, 20, 24]. When brute-force is infeasible, we use a two-phase approach: (1) global search via Halton sequence sampling [17], where the objective is evaluated at each point, and (2) local optimization of the top $n_{hj}$ candidates using a Hooke and Jeeves pattern search [18], with up to $n_{iter}$ iterations per run. This strategy follows similar hybrid methods found in the literature [22, 37].

**Limitations** The approach we propose for parameter estimation works well in practice for moderately sized case studies. It is more efficient than standard ODE-based optimization methods whenever the dimension and number of parameters gets too large. However, we are also limited by the dimension of the system since the size of the discrete state space grows exponentially with the dimension: $(n+m)^L$ if each species and parameter state is divided in $L$ subintervals. And ideally, one would want to explore the discrete state space as
much as possible when building the DBN approximation. From the numerous experiments we carried out, we concluded that a reasonable trade-off between accuracy and computation time regardless of the dimension is obtained with $L \approx 5$. Consequently, in our framework, dimensionality is the primary focus of system simplification methods. Another limiting factor is the CPTs sizes. They do not depend directly on the dimension, but on the number of parents a node depends on (referred to as the fan-in $F$ of a node, the CPT size of a node with fan-in $F$ is $(F+1)^L$ if every node is divided in $L$ sub-intervals).

A common drawback that we identified in usual biological systems reduction methods, such as classical lumping methods [45], is the increase in the variables dependencies in the system. Differential equivalence methods [8, 9], and more specifically backward differential equivalence, are another type of lumping reduction which is less subject to the increase in variables dependencies, but they still require keeping track of initial conditions. This leads to an unavoidable increase in the fan-in, and consequently a significant increase in the CPTs sizes for the reduced system, canceling any computation time gain made with the dimension reduction. Conservation analysis allows to remove variables that can be expressed as linear combinations of others. When implemented correctly, i.e., using a reduced system DBN graph structure that avoids a drastic increase in fan-in, it is possible to significantly reduce the system's dimensionality while maintaining a manageable fan-in. This, in turn, results in substantial computational time savings for both DBN construction and parameter estimation.

## 3 Model Reduction by Conservation Analysis

In the original system (2), species dynamics are determined by reaction rates. However, certain network structures induce intrinsic conservation relations [40, 45, 48], leading to linear dependencies among species and among the rows of the stoichiometric matrix $N$. Consequently, the rank of $N$ is reduced compared to the case without conservation relations. We denote this rank by $n_0$.

The conservation analysis method aims to partition species into independent and dependent subsets [48]. This is achieved by separating linearly independent and dependent rows of the stoichiometric matrix $N$. Let $n_0$ denote the number of independent species $x_i(t)$, with the remaining $n - n_0$ being dependent species $x_d(t)$. By rearranging rows, $N$ can be written as:

$$N = \begin{bmatrix} N_I \\ N_D \end{bmatrix}, \qquad (8)$$

where $N_I \in \mathbb{R}^{n_0 \times r}$ is full rank and contains the independent rows of $N$, and $N_D \in \mathbb{R}^{(n-n_0) \times r}$ contains the dependent rows of $N$.

Note that in certain networks, such a partition of the stoichiometric matrix $N$ may not be possible [40], as all species may be independent, i.e., $n = n_0$ and $N = N_I$. We propose a theoretical upper bound on $n_0$ in Appendix 7.1.
## 3.1 Link-zero Matrix $L_0$ and System Partitioning

We begin by analyzing the relationship between the submatrix of dependent species, $N_D$ and the submatrix of independent species, $N_I$.

**Proposition 1 (Link-Zero Matrix [36,40,48])**. Let $N_D$ and $N_I$ denote the dependent and independent row blocks of the stoichiometric matrix $N$, respectively. If the rows of $N_D$ are linear combinations of those in $N_I$, then there exists a matrix $L_0$ such that

$$N_D = L_0 \cdot N_I \quad (9)$$

This result enables the reformulation of the system dynamics in terms of independent and dependent species via the partitioned stoichiometric matrix.

**Corollary 1 (Partitioning of the System and Species Dynamics).** Given the stoichiometric matrix partition (8), the dynamics of the full system (2) can be rewritten as

$$\dot{x}(t) = \begin{bmatrix} N_I \\ L_0 \cdot N_I \end{bmatrix} \cdot v(t) \quad (10)$$

This form separates the dynamics of the independent and dependent species. More explicitly, the system can be expressed as:

$$\dot{x}_i(t) = N_I \cdot v(t) \quad (11)$$

$$\dot{x}_d(t) = L_0 \cdot N_I \cdot v(t) \quad (12)$$

where $x_i(t)$ and $x_d(t)$ denote the concentrations of independent and dependent species, respectively.

Building on this formulation, we next simplify the dynamics of the dependent species by expressing them directly in terms of the independent species:

**Corollary 2 (Simplified Dependent Species Dynamics).** From the partitioned equations (11) and (12), the dynamics of the dependent species reduce to

$$\dot{x}_d(t) = L_0 \cdot \dot{x}_i(t) \quad (13)$$

Integrating this equation gives the dependent species' concentrations as a function of the independent species:

$$x_d(t) = x_d(0) + L_0 \cdot [x_i(t) - x_i(0)] \quad (14)$$

where $x_d(0)$ is the initial concentration of dependent species, and $x_i(t)$ and $x_i(0)$ are the concentrations of the independent species at time $t$ and at $t = 0$, respectively.
## 3.2 Computing the Link-zero Matrix $L_0$

A variety of approaches were proposed to find the link-zero matrix $L_0$, depending on the size of the network considered as well as the nature of the equations [40,45]. These include computing the null space of the stoichiometric matrix $N$, using Gaussian elimination or SVD [45, 48]. However, such methods are often numerically unstable for very large systems. In [48], Vallabhajosyula et al. uses full pivoting QR decomposition with Householder reflections, a more numerically stable method, even for systems with large condition number. The definition of the full pivoting QR Decomposition is the following:

**Definition 3 (Full Pivoting QR Decomposition [48]).** Let $A \in \mathbb{R}^{n \times r}$ be a non-symmetric matrix. A full pivoting QR decomposition of $A$ is a factorization of the form

$$A \cdot P = Q \cdot R \qquad (15)$$

where $Q \in \mathbb{R}^{n \times n}$ is an orthogonal matrix ($Q^T Q = I_n$), $P \in \mathbb{R}^{r \times r}$ is a permutation matrix representing column exchanges in $A$, and $R \in \mathbb{R}^{n \times r}$ is an upper trapezoidal matrix (i.e., lower triangular elements are zero). The permutation matrix $P$ enhances numerical stability by reordering the columns of $A$ during the QR factorization process, typically using Householder reflections [19], as refined by Golub [15]. Full pivoting improves the decomposition's robustness by selecting optimal pivots at each step [46].

**Application to $N^T$ Matrix** To compute the link-zero matrix $L_0$, Vallabhajosyula et al. [48] use the transpose of the stoichiometric matrix $N^T$ and apply a full pivoting QR decomposition. This is because the permutation matrix $P$ reorders the columns of $N^T$ (equivalent to the rows of $N$), which is essential for generating $Q$ and $R$. The full pivoting QR decomposition of $N^T$ is:

$$N^T \cdot P_* = Q_* \cdot R_* \qquad (16)$$

where the index $*$ denotes matrices derived from the QR decomposition of $N^T$ [48]. Multiplying both sides of (16) by $Q_*^T$ yields:

$$Q_*^T \cdot N^T \cdot P_* = R_* \qquad (17)$$

Since $R_*$ is upper trapezoidal, if the rows of $N$ are linearly dependent, $R_*$ can be partitioned into 4 sub-matrices [48]. After performing a Gaussian elimination on $R_*$, each non-zero row will have a unity on the main diagonal and zeros above the diagonal, leading to:

$$R_* = \begin{bmatrix} I & M \\ 0 & 0 \end{bmatrix} \qquad (18)$$

Finally, the link-zero matrix is given by $L_0 = M^T$ [36, 40, 48].
### 3.3 Application on the Running Examples

For the enzyme-catalyzed reaction, the method reduces the original system (3) of 4 species to a subsystem with only 2 species:

$$ \begin{cases} \dot{x}_2 = -k_1x_1x_2 + (k_2 + k_3)(-x_2 + x_3(0) + x_2(0)) \\ \dot{x}_1 = -k_1x_1x_2 + k_2(-x_2 + x_3(0) + x_2(0)) \end{cases} \qquad (19) $$

The full reduction process and resulting system are detailed in Appendix 7.2.

In the case of the EGF-NGF model, the original system comprising 32 species is reduced to a subsystem with only 13 independent species, as shown in Appendix (32). This reduction is especially beneficial, as it significantly lowers the computational cost of numerical integration. For instance, running 100,000 simulations with 11 time steps using SciPy's `solve_ivp` LSODA integrator [33, 43] took 561 seconds for the original system, whereas the same setup applied to the reduced system required only 394 seconds. This numerical integration phase is approximately reduced by a factor of 1.42.

## 4 BN-based Parameter Estimation Relying on Conservation Analysis

The expression of the dependent species in (14) can be substituted into the independent species equations (11), allowing simulation of the system using only the independent variables. This reduced system is used for simulations which significantly decreases the computation time spent in the ODE solver during DBN construction as shown in Subsection 3.3. However, since our parameter estimation procedure relies on measurements that may include dependent species, we must still reconstruct the full set of species concentrations over time although they are rarely required in the parameter estimation procedure.

**DBN Construction** A key challenge is defining an appropriate DBN structure for the reduced system. As shown in Equation (14), species' initial conditions can be substituted into the system dynamics and retained in the graph. We propose four specific methods to incorporate this information. It is important to note that explicitly representing the independent variables while approximating the dependent variables through injected initial conditions may result in some loss of information.

In the *no-var* approach, the species' initial conditions are entirely omitted from the DBN, resulting in the most compact graph in terms of CPT size but at the cost of significant information loss. In contrast, the *unstacked* method retains each relevant initial condition as a distinct parameter node in the graph, preserving the full dynamics of the original system but increasing the *fan-in*. The *stacked* approach mitigates this increase by aggregating frequently recurring algebraic expressions of initial conditions (often sums, as seen in (14) and (19)) into a single parameter node. While this reduces complexity, it introduces information loss due to discretizing only the aggregated value. Finally, in the
merged method, suitable when the DBN simulation horizon is limited to a single time step (i.e. $\tau = \delta$), the initial condition of a dependent species is identical to its parent node's value, and only the latter is retained. This ensures good coverage of the discrete state space under tight simulation budgets. For longer time horizons, this method requires reconstructing intermediate values, and we propose an *ad hoc* procedure to reconstruct the initial condition distribution during unrolled inference.

In Fig. 2, we illustrate these proposed methods on the enzyme-catalyzed reaction for which the corresponding reduced system (19) is given in Subsection 3.3.

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** 2-TBN structures for the original and reduced models in the enzyme-catalyzed case. (a) Original model as depicted in [27]. (b) New parameters are treated as constants. (c) New parameters are grouped (Stacked) or treated individually (Unstacked). (d) Dependent species are computed via conservation laws at each step.</figcaption></figure>

**Reconstruction Procedure for Merged Graph Inference** We discuss here how the inference may be performed from the proposed DBN structures letting us take into account the species initial conditions. For the *stacked* and *unstacked* graphs, unrolling inference is straightforward, as only the distribution of the species at time 0 is required to estimate the distributions of the independent species at any time $t$. However, *merged* graphs require a distribution at time $t$ to estimate the distribution at time $t + 1$.

More precisely, let us denote by $X_i$ the random variables associated to the independent variables, and by $X_d$ the random variables associated to the dependent variables, we add an exponent $t$ to denote their values at (discrete) time $t$. In the *stacked* and *unstacked* graphs, we compute the distributions of species over time using:

$$p(\mathbf{X}_i^{t+1} | \mathbf{X}_i^t, \mathbf{K}, \mathbf{X}^0)$$
where the values of the species at time 0 can be viewed as parameters. We can see in Fig. 2 that $X_2^0$ and $X_3^0$ are injected as parameters $K_4$ and $K_5$.

In the *merged* graph, we estimate the distributions of species by iterating:

$$p(\mathbf{X}_i^{t+1} | \mathbf{X}_i^t, \mathbf{X}_d^t, \mathbf{K}, \mathbf{X}^0).$$

Since CPTs are computed on the interval $[0, \delta]$, $\mathbf{X}_d^t = \mathbf{X}_d^0$, this is thus the same as computing

$$p(\mathbf{X}_i^{t+1} | \mathbf{X}^t, \mathbf{K})$$

For a one step computation, there is no approximation. But to compute $\mathbf{X}_i^{t+1}$ for $t \ge 2$, we need to reconstruct a value for $\mathbf{X}_d^t$, we do this with a reconstruction function $\mathbf{X}_d^{t+1} = h(\mathbf{X}^t, \mathbf{X}_i^{t+1})$ which implements a heuristics that computes a distribution for $\mathbf{X}_d^{t+1}$ by first computing a continuous value from $\mathbf{X}^t$, $\mathbf{X}_i^{t+1}$ and then computing a random variable with realizations that fit this continuous value. The reconstruction is illustrated by the red dotted lines in Fig. 2.

**Fan-in Analysis** The maximum fan-in is generally the most critical, bottleneck inducing value, as it reflects the size of the largest CPT in the DBN. This section outlines how different reduced model graphs affect the fan-in of a node.

The *no-var* graph always reduces fan-in when conservation analysis identifies dependent species, as their removal leads to fewer incoming arcs to the remaining nodes. The *unstacked* graph significantly increases *fan-in*, often counter-balancing any dimensionality reduction. Although dependent variables are removed, a separate parameter node is introduced per species' initial condition, adding numerous new arcs. The *stacked* graph leads to a null or moderate increase in *fan-in*. While dependent variables are removed, they are replaced by a new node representing the species' initial conditions sums (see (14)). As a result, some independent species may depend on more variables than in the original model. The *merged* graph preserves the original *fan-in*, and all the information about species' initial conditions is kept for inference. Overall, the *stacked* graph is the most broadly applicable, as it requires no hypothesis on the DBN construction, whereas the *merged* graph offers the best balance between reduction and inference accuracy.

# 5 Numerical Results

**Implementation** The approach was implemented in Python and is available online as a prototype called BayeSBML [2]. Biochemical networks were specified using SBML files via SimpleSBML [7, 39]. CPTs were constructed with pgmpy [34], and Bayesian network modeling and inference were performed using pyAgrum [35]. To improve efficiency, we applied Just-In-Time (JIT) compilation with Numba [30], achieving near-C performance.
All computations were performed on a standard laptop<fnref n="1" />, except for CPTs generation in the EGF-NGF case, which required a more powerful machine<fnref n="2" />. The code and data enabling the reproduction of the results presented in this paper are available as a permanent artifact online [3].

We now compare the reduced models presented in Section 4 with the original model using our running examples.

**Enzyme-catalyzed Reaction Network** The CPTs are generated with 78,125 ODE simulations via SciPy's `solve_ivp` LSODA integrator [33,43], over 11 time steps.

<table><caption>**Table 1.** Comparison of the original and reduced models for the enzyme-catalyzed reaction network example.</caption><thead><tr><th>Model</th><th>Original</th><th>No-var</th><th>Stacked</th><th>Unstacked</th><th>Merged</th></tr></thead><tbody><tr><td>argmin<sub>K</sub> ε(K) := K̃</td><td>(3,4,1)</td><td>(3,2,2)</td><td>(3,2,1)</td><td>(1,1,1)</td><td>(2,2,1)</td></tr><tr><td>min<sub>K</sub> ε(K) := ε(K̃)</td><td>0.0538</td><td>0.0604</td><td>0.0502</td><td>0.0431</td><td>0.0481</td></tr><tr><td>State-space Size</td><td>78125</td><td>3125</td><td>31250</td><td>78125</td><td>15625</td></tr><tr><td>Total CPTs Rows</td><td>172500</td><td>18750</td><td>187500</td><td>468750</td><td>93750</td></tr><tr><td>fan-in (Cum., Max)</td><td>(20,6)</td><td>(9,5)</td><td>(11,6)</td><td>(13,7)</td><td>(11,6)</td></tr></tbody></table>

Due to the system's low dimensionality, the increased *fan-in* in the *unstacked* model is not substantial. As reported in Table 1, this is the only model with a higher maximum *fan-in*, despite a reduced number of species. All other reduced models achieve lower *fan-in* values, both cumulative and maximal, compared to the original. As expected, the reduction in state space size is most significant for the *no-var* model, with an approximate 20-fold decrease. The state space is reduced by a factor of 5 in the *merged* model and by a factor of 2.5 in the *stacked* model, while no reduction is observed in the *unstacked* model.

We also observe that the total CPTs size (summed across all nodes) increases significantly for the *unstacked* model. In the *no-var* and *merged* models, some nodes are entirely removed, resulting in smaller CPTs. For the *stacked* model, although two variables are eliminated, a new variable representing species' initial conditions sum is introduced. This variable, discretized 10 times (5 per original variable), appears in all reduced equations, slightly increasing the overall CPTs size. A similar mechanism occurs in the *unstacked* model, but since a separate variable is introduced for each initial condition, the increase in CPTs size is more pronounced.

<footnote n="1">¹ Asus Zenbook 15 equipped with an 8-Core AMD Ryzen™ 7 7735U processor and 32 GB of RAM</footnote>

<footnote n="2">² 16-Core AMD EPYC 7281 processor with 128 GB of RAM</footnote>
To evaluate the predictive quality of the models and the accuracy of their DBN approximations, we carry out a parameter estimation procedure identical to that described in Subsection 2.3. The continuous parameter space is discretized into 5 intervals per parameter as described in Subsection 2.2, resulting in a grid of $5^3 = 125$ possible configurations. Next, we generate a reference trajectory with a known initial condition $x_0 = (5, 1, 7.5, 1.5)$ and a continuous parameter vector $k = (0.3, 0.3, 0.3)$, $K = (1, 1, 1)$ if discretized, that we aim to recover. A brute-force search is then performed over this discrete space to identify the optimal parameter configuration, $\tilde{K}$ that minimizes the objective function defined in (7). Once the optimal discrete combination $\tilde{K}$ is found, it is mapped back to a continuous estimate $\tilde{k}$ by taking the midpoints of the corresponding intervals. This estimated parameter vector is then used to simulate the model and assess its predictive performance.

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Species concentrations in the enzyme-catalyzed reaction example. ODE$_k$: simulation with nominal parameters; ODE$_{\tilde{k}}$: simulation with learned parameters through LSODA scheme and 5 discrete states; BN$_K$: the expected values for the proposed BN-based approximation with learned parameters.</figcaption></figure>

Fig. 3 compares the original ODE trajectory (ODE$_k$), the ODE trajectory using the estimated parameters (ODE$_{\tilde{k}}$), and the BN-derived trajectory based on $\tilde{K}$. The latter is indeed the expected value of the state at each time step using the midpoint-weighted distributions as described in (6).

As shown in Table 1, all reduced models, except *no-var*, match or outperform the original in predictive accuracy. However, the estimated parameter configurations $\tilde{K}$ differ notably from the reference $K = (1, 1, 1)$ for most models, highlighting a well-known issue in biochemical modeling: *parameter identifiability* [14, 44, 49, 50]. That is, distinct parameter sets can yield similar system behavior, complicating reliable parameter inference.
**The EGF-NGF Model** The CPTs are generated using 1,000,000 ODE simulations over 2 time steps. The high dimensionality makes the *unstacked* version impractical, with a state space size of $5^{84}$ (vs. $5^{52}$ for the original), resulting in massive CPTs of size 5.9 Go (vs 10.6 Mo for the original) and infeasible inference times, due to the addition of a parameter node per species' initial condition.

<table><caption>**Table 2.** Comparison of the original and reduced models for the EGF-NGF case study.</caption><thead><tr><th>Model</th><th>Original</th><th>No-var</th><th>Stacked</th><th>Unstacked</th><th>Merged</th></tr></thead><tbody><tr><td>Total CPTs - Rows</td><td>1.01e6</td><td>3.31e4</td><td>3.96e7</td><td>5.69e8</td><td>4.66e5</td></tr><tr><td>Total CPTs - Size (Mo)</td><td>10.62</td><td>0.35</td><td>415.68</td><td>5967.24</td><td>4.88</td></tr><tr><td>fan-in (Cum., Max)</td><td>(206,11)</td><td>(77,10)</td><td>(108,13)</td><td>(131,16)</td><td>(94,11)</td></tr></tbody></table>

The *merged* and *no-var* models are efficient with significantly smaller CPTs, as shown in Table 2. The *stacked* model is more complex but still usable. We focus on the *original*, *merged*, *stacked*, and *no-var* models, excluding the *unstacked* model due to its prohibitive computational cost.

Regarding structural complexity, the cumulative *fan-in* in the reduced models is up to half that of the original (see Table 2). The *stacked* model has a higher maximum *fan-in* than the original, leading to prohibitive inference costs, while the *merged* and *no-var* models maintain a maximum *fan-in* equal to or smaller than the original.

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** boundNGFR($x_6$), inactive Akt($x_{24}$) and inactive PI3K($x_{22}$) concentrations in the EGF-NGF reaction example. ODE$_k$ simulation with nominal parameters; ODE$_{\bar{k}}$: simulation with learned parameters through LSODA scheme and 5 discrete states; BN$_K$: the expected values for the proposed BN-based approximation with learned parameters.</figcaption></figure>
For parameter estimation, we generated a reference trajectory with $\boldsymbol{x}_0$ and $\boldsymbol{k}$ specified in the Appendix Tables 3 and 4. The continuous parameter space is discretized into 5 intervals, resulting in a grid of $5^{20}$ possible configurations. Due to the large size, we employ the exploration/exploitation method outlined in Section 2.3, using 250 Halton sequence samples, $n_{hj} = 3$, and $n_{iter} = 25$.

We selected key signaling species from the MAPK/ERK and PI3K-Akt pathways, which are the primary cascades activated by EGF and NGF stimulation [27, 31, 38, 42]. The analysis focuses on the dynamics of boundNGFR ($x_6$), and the inactive forms of PI3K ($x_{22}$) and Akt ($x_{24}$), due to their central roles in signal transduction and feedback regulation.

Similarly to the enzyme-catalyzed case, Fig. 4 compares the original ODE trajectory ($ODE_k$), the ODE trajectory using the estimated parameters ($ODE_{\bar{k}}$), and the BN-derived trajectory based on $\bar{K}$.

For both running examples, we observe comparable behaviors, except for the *stacked* model that struggles due to its large state space in the EGF-NGF case (see Akt* in Fig. 4).

# 6 Conclusion

The proposed approach leverages a model simplification technique that directly impacts the structure of our probabilistic approximation. Though previously unexplored, it yields significant performance gains, enhancing the method's applicability to high-dimensional case studies.

In future work, we aim to address the parameter identifiability problem to ensure robust and reliable estimation. We also plan to collaborate with biologists to apply our methods to real-world data, further validating their practical applicability. Other possible directions include parameter identification of more complex models, such as partial differential equations, and better sampling techniques.

Regarding implementation, we intend to integrate efficient CPU parallelization and batching into our existing framework BayeSBML to accelerate numerical simulations (and thus reduce the CPTs construction time) as well as the parameter optimization procedure.
References

1. Banerjee, S., Roy, A.: Linear Algebra and Matrix Analysis for Statistics, pp. 129–137. Chapman and Hall/CRC (2014)

2. bayesbml, https://git.lacl.fr/barbot/bayesbml.git

3. Bouët-Willaumez, O., Barbot, B., Le Coënt, A., Pekergin, N.: Artifact: Conservation analysis and discrete probabilistic approximations for parameter estimation of biochemical networks (2025). https://doi.org/10.5281/zenodo.15284362

4. Brenig, L.: Reducing nonlinear dynamical systems to canonical forms. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 376(2124), 20170384 (2018)

5. Brown, K.S., Hill, C.C., Calero, G.A., Myers, C.R., Lee, K.H., Sethna, J.P., Cerione, R.A.: The statistical mechanics of complex signaling networks: nerve growth factor signaling. Physical biology 1(3), 184 (2004)

6. Caflisch, R.E., Morokoff, W., Owen, A.B.: Valuation of mortgage-backed securities using brownian bridges to reduce effective dimension. Journal of Computational Finance 1(1), 27–46 (1997)

7. Cannistra, C., Medley, K., Sauro, H.M.: Simplesbml: A python package for creating and editing sbml models. bioRxiv (2015), https://apisemanticscholar.org/CorpusID:61794957

8. Cardelli, L., Tribastone, M., Tschaikowski, M., Vandin, A.: Symbolic computation of differential equivalences. ACM SIGPLAN Notices 51(1), 137–150 (2016)

9. Cardelli, L., Tribastone, M., Tschaikowski, M., Vandin, A.: Erode: a tool for the evaluation and reduction of ordinary differential equations. In: Tools and Algorithms for the Construction and Analysis of Systems: 23rd International Conference, TACAS 2017, Held as Part of the European Joint Conferences on Theory and Practice of Software, ETAPS 2017, Uppsala, Sweden, April 22-29, 2017, Proceedings, Part II 23. pp. 310–328. Springer (2017)

10. Cardelli, L., Tribastone, M., Tschaikowski, M., Vandin, A.: Maximal aggregation of polynomial dynamical systems. Proceedings of the National Academy of Sciences 114(38), 10029–10034 (2017). https://doi.org/10.1073/pnas.1702697114, https://www.pnas.org/doi/abs/10.1073/pnas.1702697114

11. Cardelli, L., Tribastone, M., Tschaikowski, M., Vandin, A.: Symbolic computation of differential equivalences. Theoretical Computer Science 777, 132–154 (2019). https://doi.org/10.1016/j.tcs.2019.03.018, https://www.sciencedirect.com/science/article/pii/S0304397519301719

12. Cohen, S.: Origins of growth factors: NGF and EGF. Journal of Biological Chemistry 283(49), 33793–33797 (2008). https://doi.org/10.1074/jbc.X800008200

13. Dokoumetzidis, A., Aarons, L.: Proper lumping in systems biology models. IET systems biology 3(1), 40–51 (2009)

14. Eisenberg, M.C., Hayashi, M.A.: Determining identifiable parameter combinations using subset profiling. Mathematical Biosciences 256, 116–126 (2014). https://doi.org/10.1016/j.mbs.2014.08.008, https://www.sciencedirect.com/science/article/pii/S0025556414001631

15. Golub, G.H., Van Loan, C.F.: Matrix Computations. The Johns Hopkins University Press (1996)

16. Gratie, D.E., Iancu, B., Petre, I.: ODE Analysis of Biological Systems, pp. 29–62. Springer Berlin Heidelberg (2013)

17. HALTON, J.: On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals. Numerische Mathematik 2, 84–90 (1960), http://eudml.org/doc/131448
18. Hooke, R., Jeeves, T.A.: Direct search solution of numerical and statistical problems. Journal of the ACM (JACM) **8**(2), 212–229 (1961)

19. Householder, A.S.: Unitary triangularization of a nonsymmetric matrix. J. ACM **5**(4), 339–342 (1958)

20. Huyer, W., Neumaier, A.: Global optimization by multilevel coordinate search. Journal of Global Optimization **14**(4), 331–355 (1999)

21. Jiménez-Pastor, A., Jacob, J.P., Pogudin, G.: Exact linear reduction for rational dynamical systems. In: International Conference on Computational Methods in Systems Biology. pp. 198–216. Springer (2022)

22. Kang, F., Li, J., Li, H.: Artificial bee colony algorithm and pattern search hybridized for global optimization. Applied Soft Computing **13**(4), 1781–1791 (2013). https://doi.org/10.1016/j.asoc.2012.12.025, https://www.sciencedirect.com/science/article/pii/S1568494612005698

23. Klonowski, W.: Simplifying principles for chemical and enzyme reaction kinetics. Biophysical chemistry **18**(2), 73–87 (1983)

24. Kolda, T.G., Lewis, R.M., Torczon, V.: Optimization by direct search: New perspectives on some classical and modern methods. SIAM review **45**(3), 385–482 (2003)

25. Koller, D., Friedman, N.: Probabilistic Graphical Models - Principles and Techniques. MIT Press (2009)

26. Le Coënt, A., Barbot, B., Pekergin, N.: Efficient probabilistic inference in biochemical networks. Computers in Biology and Medicine **183**, 109280 (2024). https://doi.org/10.1016/j.compbiomed.2024.109280, https://www.sciencedirect.com/science/article/pii/S0010482524013659

27. Liu, B., Hsu, D., Thiagarajan, P.S.: Probabilistic approximations of odes-based bio-pathway dynamics. Theoretical Computer Science **412**(21), 2188–2206 (2011)

28. Maarleveld, T., Khandelwal, R., Olivier, B., Teusink, B., Bruggeman, F.: Basic concepts and principles of stoichiometric modeling of metabolic networks. Biotechnology Journal **8**(9), 997–1008 (2013)

29. Malik-Sheriff, R.S., Glont, M., et al.: Biomodels—15 years of sharing computational models in life science. Nucleic Acids Research **48**(D1), D407–D415 (2019)

30. numba, https://numba.pydata.org/

31. Okada-Hatakeyama, M., et al.: A computational model on the modulation of mitogen-activated protein kinase (mapk) and akt pathways in heregulin-induced erbb signalling. The Biochemical journal **373**, 451–63 (08 2003). https://doi.org/10.1042/BJ20021824

32. Ovchinnikov, A., Pérez Verona, I., Pogudin, G., Tribastone, M.: Clue: exact maximal reduction of kinetic models by constrained lumping of differential equations. Bioinformatics **37**(12), 1732–1738 (2021)

33. Petzold, L.: Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations. SIAM Journal on Scientific and Statistical Computing **4** (1983)

34. pgmpy, https://pgmpy.org/

35. pyagrum, https://agrum.gitlab.io/

36. Reder, C.: Metabolic control theory: A structural approach. Journal of Theoretical Biology **135**(2), 175–201 (1988)

37. Rios-Coelho, A., Sacco, W., Henderson, N.: A metropolis algorithm combined with hooke-jeeves local search method applied to global optimization. Applied Mathematics and Computation **217**(2), 843–853 (2010). https://doi.org/10.1016/j.amc.2010.06.027, https://www.sciencedirect.com/science/article/pii/S0096300310007125
38. Santos, S.D.M., Verveer, P.J., Bastiaens, P.I.H.: Growth factor-induced mapk network topology shapes erk response determining pc-12 cell fate. Nature Cell Biology 9(3), 324-330 (2007)

39. Sauro, H.M.: Simplesbml: A python package for creating, editing, and interrogating sbml models: Version 2.0 (2021), https://arxiv.org/abs/2009.01969

40. Sauro, H.M., Ingalls, B.: Conservation analysis in biochemical networks: Computational issues for software writers. Biophysical Chemistry 109(1), 1-15 (2004)

41. Schauer, M., Heinrich, R.: Quasi-steady-state approximation in the mathematical modeling of biochemical reaction networks. Mathematical biosciences 65(2), 155-170 (1983)

42. Schoeberl, B., Eichler-Jonsson, C., Gilles, E., Müller, G.: Computational modeling of the dynamics of the map kinase cascade activated by surface and internalized egf receptors. Nature biotechnology 20, 370-5 (05 2002). https://doi.org/10.1038/nbt0402-370

43. scipy.integrate, https://docs.scipy.org/doc/scipy/reference/integrate.html

44. Simpson, M.J., Baker, R.E.: Parameter identifiability, parameter estimation and model prediction for differential equation models (2025), https://arxiv.org/abs/2405.08177

45. Snowden, T.J., van der Graaf, P.H., Tindall, M.J.: Methods of model reduction for large-scale biological systems: A survey of current methods and trends. Bulletin of Mathematical Biology 79(7), 1449-1486 (2017)

46. Strang, G.: Linear Algebra and Learning from Data. Cambridge University Press (2019)

47. Valderrama-Bahamóndez, G.I., Fröhlich, H.: Mcmc techniques for parameter estimation of ode based models in systems biology. Frontiers in Applied Mathematics and Statistics 5, 55 (2019)

48. Vallabhajosyula, R.R., Chickarmane, V., Sauro, H.M.: Conservation analysis of large biochemical networks. Bioinformatics 22(3), 346-353 (2006)

49. Villaverde, A., Barreiro, A.: Identifiability of large nonlinear biochemical networks. MATCH Communications in Mathematical and in Computer Chemistry 76, 259-296 (06 2016)

50. Wieland, F.G., Hauber, A.L., Rosenblatt, M., Tönsing, C., Timmer, J.: On structural and practical identifiability. Current Opinion in Systems Biology 25, 60-69 (Mar 2021). https://doi.org/10.1016/j.coisb.2021.03.005, http://dx.doi.org/10.1016/j.coisb.2021.03.005
# 7 Appendix

## 7.1 The Rank Theorem and the Upper-bound on $n_0$ Value

**Theorem 1 (Fundamental Theorem of Ranks [1]).** Let $A \in \mathbb{R}^{n \times r}$ be a non-symmetric matrix. Then, the rank of $A$ satisfies the following equality:

$$ \mathrm{rank}(A) = \mathrm{rank}(A^T A) = \mathrm{rank}(AA^T) = \mathrm{rank}(A^T) \quad (20) $$

*Proof.* Refer to [1], Theorem 5.2, page 132. □

**Corollary 3.** As a direct consequence of Theorem 1, the rank of a matrix $A$ is equal to the number of linearly independent rows or columns of $A$ [1]. Furthermore, since the number of linearly independent rows (resp. columns) cannot exceed $n$ (resp. $r$), we obtain the bound:

$$ \mathrm{rank}(A) \leq \min(n, r) \quad (21) $$

Applying this result to the stoichiometric matrix of a biochemical network, we can establish an upper bound on its rank. Specifically, this leads to the following theorem:

**Theorem 2.** Let $N \in \mathbb{R}^{n \times r}$ be the stoichiometric matrix associated with a biochemical network. If the number of reactions $r$ in the network is strictly less than the number of species $n$, then:

$$ \mathrm{rank}(N) = n_0 \leq r < n \quad (22) $$

This implies that biochemical systems with significantly more species than reactions can be efficiently reduced using the conservation analysis method.

*Proof.* Suppose that $r < n$, then $\min(n, r) = r$. Using Corollary 3, it comes that $\mathrm{rank}(N) = n_0 \leq \min(n, r) = r < n$. □

## 7.2 Application of the Conservation Analysis Method to the Enzyme-catalyzed Reaction

For this example, performing a full pivoting Householder QR decomposition of $N^T$, followed by a Gaussian elimination, allows us to obtain the following $R_*$ and $P_*$ matrices:

$$ R_* = \begin{bmatrix} 1 & 0 & -1 & 1 \\ 0 & 1 & 0 & -1 \\ 0 & 0 & 0 & 0 \end{bmatrix} = \left[ \frac{I}{0} \middle| \frac{M}{0} \right], \quad P_* = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \quad (23) $$
In this case, the rank of $R_*$ is 2, so it means $n_0 = 2$. Using (18), it follows upon that

$$L_0 = M^T = \begin{bmatrix} -1 & 0 \\ 1 & -1 \end{bmatrix} \qquad (24)$$

Now, we make use of the permutation matrix $P_*$ (23) to obtain the new reordered species list and stoichiometry matrix. Specifically, the reordered species list is the following

$$x_{\text{new}}(t) = (x_1(t), x_2(t), x_3(t), x_4(t)) \cdot P_* = (x_2(t), x_1(t), x_3(t), x_4(t)) \quad (25)$$

It comes that the first $n_0$ species are the independent ones, and the last $n - n_0$ are the dependent ones, i.e.

$$x_i(t) = (x_2(t), x_1(t)), \quad x_d(t) = (x_3(t), x_4(t)) \qquad (26)$$

Similarly, we extract the matrices $N_I$ and $N_D$ as sub-matrices of the reordered stoichiometry matrix (first $n_0$ and last $n - n_0$ rows), such as

$$N_I = \begin{bmatrix} -1 & 1 & 1 \\ -1 & 1 & 0 \end{bmatrix}, \quad N_D = \begin{bmatrix} 1 & -1 & -1 \\ 0 & 0 & 1 \end{bmatrix} \qquad (27)$$

Subsequently, using Corollary 1, we compute the independent species $x_i$ sub-system as follows

$$\dot{x}_i(t) = N_I \cdot v(t) \rightarrow \begin{cases} \dot{x}_2 = -k_1 x_1 x_2 + (k_2 + k_3)x_3 \\ \dot{x}_1 = -k_1 x_1 x_2 + k_2 x_3 \end{cases} \qquad (28)$$

Therefore, by leveraging the link-zero matrix $L_0$ along with Corollary 2, we can obtain an explicit expression for the dependent species $x_d$ in terms of the independent species $x_i$ and the species initial concentrations:

$$x_d(t) = \begin{bmatrix} x_3(t) \\ x_4(t) \end{bmatrix} = \begin{bmatrix} x_3(0) \\ x_4(0) \end{bmatrix} + L_0 \cdot \begin{bmatrix} x_2(t) - x_2(0) \\ x_1(t) - x_1(0) \end{bmatrix} \qquad (29)$$

$$= \begin{bmatrix} -x_2(t) \\ x_2(t) - x_1(t) \end{bmatrix} + \begin{bmatrix} x_3(0) + x_2(0) \\ x_1(0) - x_2(0) + x_4(0) \end{bmatrix} \qquad (30)$$

Finally, using the expressions of the dependent species we just obtained, we can simplify the sub-system (28) so that it only depends on the independent species $x_i$, the parameters, and the species initial concentrations, which leads to the sub-system (19), i.e.:

$$\begin{cases} \dot{x}_2 = -k_1 x_1 x_2 + (k_2 + k_3) * (-x_2 + x_3(0) + x_2(0)) \\ \dot{x}_1 = -k_1 x_1 x_2 + k_2 * (-x_2 + x_3(0) + x_2(0)) \end{cases}$$
$$ \dot{x}_1 = -k_1 * x_1 * x_3 + k_2 * x_4 $$

$$ \dot{x}_2 = -k_3 * x_2 * x_5 + k_4 * x_6 $$

$$ \dot{x}_3 = -k_1 * x_1 * x_3 + k_2 * x_4 $$

$$ \dot{x}_4 = k_1 * x_1 * x_3 - k_2 * x_4 $$

$$ \dot{x}_5 = -k_3 * x_2 * x_5 + k_4 * x_6 $$

$$ \dot{x}_6 = k_3 * x_2 * x_5 - k_4 * x_6 $$

$$ \dot{x}_7 = k_9 * x_{10} * x_8/(x_8 + k_{10}) - k_5 * x_4 * x_7/(x_7 + k_6) - k_7 * x_6 * x_7/(x_7 + k_8) $$

$$ \dot{x}_8 = -k_9 * x_{10} * x_8/(x_8 + k_{10}) + k_5 * x_4 * x_7/(x_7 + k_6) + k_7 * x_6 * x_7/(x_7 + k_8) $$

$$ \dot{x}_9 = -k_{27} * x_{21} * x_9/(x_9 + k_{28}) $$

$$ \dot{x}_{10} = k_{27} * x_{21} * x_9/(x_9 + k_{28}) $$

$$ \dot{x}_{11} = -k_{11} * x_8 * x_{11}/(x_{11} + k_{12}) + k_{13} * x_{13} * x_{12}/(x_{12} + k_{14}) $$

$$ \dot{x}_{12} = k_{11} * x_8 * x_{11}/(x_{11} + k_{12}) - k_{13} * x_{13} * x_{12}/(x_{12} + k_{14}) $$

$$ \dot{x}_{13} = 0 $$

$$ \dot{x}_{14} = -k_{15} * x_{12} * x_{14}/(x_{14} + k_{16}) + k_{45} * x_{32} * x_{15}/(x_{15} + k_{46}) + k_{35} * x_{25} * x_{15}/(x_{15} + k_{36}) $$

$$ \dot{x}_{15} = k_{15} * x_{12} * x_{14}/(x_{14} + k_{16}) - k_{45} * x_{32} * x_{15}/(x_{15} + k_{46}) - k_{35} * x_{25} * x_{15}/(x_{15} + k_{36}) $$

$$ \dot{x}_{16} = -k_{43} * x_{29} * x_{16}/(x_{16} + k_{44}) + k_{47} * x_{32} * x_{17}/(x_{17} + k_{20}) $$

$$ \dot{x}_{17} = k_{43} * x_{29} * x_{16}/(x_{16} + k_{44}) - k_{47} * x_{32} * x_{17}/(x_{17} + k_{20}) $$

$$ \dot{x}_{18} = -k_{17} * x_{15} * x_{18}/(x_{18} + k_{18}) - k_{19} * x_{17} * x_{18}/(x_{18} + k_{48}) + k_{21} * x_{31} * x_{19}/(x_{19} + k_{22}) $$

$$ \dot{x}_{19} = k_{17} * x_{15} * x_{18}/(x_{18} + k_{18}) + k_{19} * x_{17} * x_{18}/(x_{18} + k_{48}) - k_{21} * x_{31} * x_{19}/(x_{19} + k_{22}) $$

$$ \dot{x}_{20} = -k_{23} * x_{19} * x_{20}/(x_{20} + k_{24}) + k_{25} * x_{31} * x_{21}/(x_{21} + k_{26}) $$

$$ \dot{x}_{21} = k_{23} * x_{19} * x_{20}/(x_{20} + k_{24}) - k_{25} * x_{31} * x_{21}/(x_{21} + k_{26}) $$

$$ \dot{x}_{22} = -k_{29} * x_4 * x_{22}/(x_{22} + k_{30}) - k_{31} * x_{12} * x_{22}/(x_{22} + k_{32}) $$

$$ \dot{x}_{23} = k_{29} * x_4 * x_{22}/(x_{22} + k_{30}) + k_{31} * x_{12} * x_{22}/(x_{22} + k_{32}) $$

$$ \dot{x}_{24} = -k_{33} * x_{23} * x_{24}/(x_{24} + k_{34}) $$

$$ \dot{x}_{25} = k_{33} * x_{23} * x_{24}/(x_{24} + k_{34}) $$

$$ \dot{x}_{26} = -k_{37} * x_6 * x_{26}/(x_{26} + k_{38}) $$

$$ \dot{x}_{27} = k_{37} * x_6 * x_{26}/(x_{26} + k_{38}) $$

$$ \dot{x}_{28} = -k_{39} * x_{27} * x_{28}/(x_{28} + k_{40}) + k_{41} * x_{30} * x_{29}/(x_{29} + k_{42}) $$

$$ \dot{x}_{29} = k_{39} * x_{27} * x_{28}/(x_{28} + k_{40}) - k_{41} * x_{30} * x_{29}/(x_{29} + k_{42}) $$

$$ \dot{x}_{30} = 0 $$

$$ \dot{x}_{31} = 0 $$

$$ \dot{x}_{32} = 0 $$

(31)
<figure><img src="image_8.png" /><figcaption>Fig. 5. The EGF-NGF signaling pathway [27].</figcaption></figure>
$$ \dot{x}_7 = -(k_5 * x_4 * x_7/(x_7 + k_6)) - (k_7 * (x_6(0) + x_5(0) - x_5) * x_7/(x_7 + k_8)) + (k_9 * (x_{10}(0) + x_9(0) - x_9) * (x_8(0) + x_7(0) - x_7)/((x_8(0) + x_7(0) - x_7) + k_{10})) $$

$$ \begin{aligned} \dot{x}_{14} = & -(k_{15} * (x_{12}(0) + x_{11}(0) - x_{11}) * x_{14}/(x_{14} + k_{16})) \\ & + (k_{45} * x_{32} * (x_{15}(0) + x_{14} - x_{14})/((x_{15}(0) + x_{14} - x_{14}) + k_{46})) \\ & + (k_{35} * (x_{25}(0) + x_{24}(0) - x_{24}) * (x_{15}(0) + x_{14} - x_{14})/((x_{15}(0) + x_{14} - x_{14}) + k_{36})) \end{aligned} $$

$$ \begin{aligned} \dot{x}_{18} = & -(k_{17} * (x_{15}(0) + x_{14}(0) - x_{14}) * x_{18}/(x_{18} + k_{18})) \\ & - (k_{19} * (x_{17}(0) + x_{16}(0) - x_{16}) * x_{18}/(x_{18} + k_{20})) \\ & + (k_{21} * x_{31} * (x_{19}(0) + x_{18} - x_{18})/((x_{19}(0) + x_{18} - x_{18}) + k_{22})) \end{aligned} $$

$$ \dot{x}_4 = k_1 * (x_1(0) + x_4(0) - x_4) * (x_3(0) + x_4(0) - x_4) - (k_2 * x_4) $$

$$ \dot{x}_5 = -(k_3 * (x_2(0) - x_5(0) + x_5) * x_5) + (k_4 * (x_6(0) + x_5(0) - x_5)) $$

$$ \begin{aligned} \dot{x}_{11} = & -(k_{11} * (x_8(0) + x_7(0) - x_7) * x_{11}/(x_{11} + k_{12})) \\ & + (k_{13} * x_{13} * (x_{12}(0) + x_{11}(0) - x_{11})/((x_{12}(0) + x_{11}(0) - x_{11}) + k_{14})) \end{aligned} $$

$$ \begin{aligned} \dot{x}_{16} = & -(k_{43} * (x_{29}(0) + x_{28}(0) - x_{28}) * x_{16}/(x_{16} + k_{44})) \\ & + (k_{47} * x_{32} * (x_{17}(0) + x_{16} - x_{16})/((x_{17}(0) + x_{16} - x_{16}) + k_{48})) \end{aligned} $$

$$ \begin{aligned} \dot{x}_{20} = & -(k_{23} * (x_{19}(0) + x_{18}(0) - x_{18}) * x_{20}/(x_{20} + k_{24})) \\ & + (k_{25} * x_{31} * (x_{21}(0) + x_{20} - x_{20})/((x_{21}(0) + x_{20} - x_{20}) + k_{26})) \end{aligned} $$

$$ \dot{x}_{22} = -(k_{29} * x_4 * x_{22}/(x_{22} + k_{30})) - (k_{31} * (x_{12}(0) + x_{11}(0) - x_{11}) * x_{22}/(x_{22} + k_{32})) $$

$$ \begin{aligned} \dot{x}_{28} = & -(k_{39} * (x_{27}(0) + x_{26}(0) - x_{26}) * x_{28}/(x_{28} + k_{40})) \\ & + (k_{41} * x_{30} * (x_{29}(0) + x_{28} - x_{28})/((x_{29}(0) + x_{28} - x_{28}) + k_{42})) \end{aligned} $$

$$ \dot{x}_9 = -(k_{27} * (x_{21}(0) + x_{20}(0) - x_{20}) * x_9/(x_9 + k_{28})) $$

$$ \dot{x}_{24} = -(k_{33} * (x_{23}(0) + x_{22}(0) - x_{22}) * x_{24}/(x_{24} + k_{34})) $$

$$ \dot{x}_{26} = -(k_{37} * (x_6(0) + x_5(0) - x_5) * x_{26}/(x_{26} + k_{38})) $$

$$ x_{15} = x_{15}(0) + x_{14}(0) - x_{14} $$

$$ x_{21} = x_{21}(0) + x_{20}(0) - x_{20} $$

$$ x_{23} = x_{23}(0) + x_{22}(0) - x_{22} $$

$$ x_{17} = x_{17}(0) + x_{16}(0) - x_{16} $$

$$ x_{19} = x_{19}(0) + x_{18}(0) - x_{18} $$

$$ x_{29} = x_{29}(0) + x_{28}(0) - x_{28} $$

$$ x_8 = x_8(0) + x_7(0) - x_7 $$

$$ x_1 = x_1(0) + x_4(0) - x_4 $$

$$ x_{12} = x_{12}(0) + x_{11}(0) - x_{11} $$

$$ x_2 = x_2(0) - x_5(0) + x_5 $$

$$ x_3 = x_3(0) + x_4(0) - x_4 $$

$$ x_6 = x_6(0) + x_5(0) - x_5 $$

$$ x_{13} = x_{13}(0) $$

$$ x_{27} = x_{27}(0) + x_{26}(0) - x_{26} $$

$$ x_{10} = x_{10}(0) + x_9(0) - x_9 $$

$$ x_{25} = x_{25}(0) + x_{24}(0) - x_{24} $$

$$ x_{30} = x_{30}(0) $$

$$ x_{31} = x_{31}(0) $$

$$ x_{32} = x_{32}(0) $$
<table><caption>**Table 3.** Initial concentrations, names and bounds for species in the EGF-NGF signaling pathway model [5, 27, 29].</caption><thead><tr><th>Species</th><th>Name</th><th>Initial value</th><th>Range</th></tr></thead><tbody><tr><td>x<sub>1</sub></td><td>EGF</td><td>1.0002 × 10<sup>6</sup></td><td>[0, 1.1002 × 10<sup>6</sup>]</td></tr><tr><td>x<sub>2</sub></td><td>NGF</td><td>4.56 × 10<sup>5</sup></td><td>[0, 5.016 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>3</sub></td><td>freeEGFReceptor</td><td>8.0 × 10<sup>4</sup></td><td>[0, 8.8 × 10<sup>4</sup>]</td></tr><tr><td>x<sub>4</sub></td><td>boundEGFReceptor</td><td>0</td><td>[0, 8.8 × 10<sup>4</sup>]</td></tr><tr><td>x<sub>5</sub></td><td>freeNGFReceptor</td><td>1.0 × 10<sup>4</sup></td><td>[0, 1.1 × 10<sup>4</sup>]</td></tr><tr><td>x<sub>6</sub></td><td>boundNGFReceptor</td><td>0</td><td>[0, 1.1 × 10<sup>4</sup>]</td></tr><tr><td>x<sub>7</sub></td><td>SosInactive</td><td>1.2 × 10<sup>5</sup></td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>8</sub></td><td>SosActive</td><td>0</td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>9</sub></td><td>P90RskInactive</td><td>1.2 × 10<sup>5</sup></td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>10</sub></td><td>P90RskActive</td><td>0</td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>11</sub></td><td>RasInactive</td><td>1.2 × 10<sup>5</sup></td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>12</sub></td><td>RasActive</td><td>0</td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>13</sub></td><td>RasGapActive</td><td>1.2 × 10<sup>5</sup></td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>14</sub></td><td>Raf1Inactive</td><td>1.2 × 10<sup>5</sup></td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>15</sub></td><td>Raf1Active</td><td>0</td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>16</sub></td><td>BRafInactive</td><td>1.2 × 10<sup>5</sup></td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>17</sub></td><td>BRafActive</td><td>0</td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>18</sub></td><td>MekInactive</td><td>6.0 × 10<sup>5</sup></td><td>[0, 6.6 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>19</sub></td><td>MekActive</td><td>0</td><td>[0, 6.6 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>20</sub></td><td>ErkInactive</td><td>6.0 × 10<sup>5</sup></td><td>[0, 6.6 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>21</sub></td><td>ErkActive</td><td>0</td><td>[0, 6.6 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>22</sub></td><td>PI3KInactive</td><td>1.2 × 10<sup>5</sup></td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>23</sub></td><td>PI3KActive</td><td>0</td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>24</sub></td><td>AktInactive</td><td>1.2 × 10<sup>5</sup></td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>25</sub></td><td>AktActive</td><td>0</td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>26</sub></td><td>C3GInactive</td><td>1.2 × 10<sup>5</sup></td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>27</sub></td><td>C3GActive</td><td>0</td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>28</sub></td><td>Rap1Inactive</td><td>1.2 × 10<sup>5</sup></td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>29</sub></td><td>Rap1Active</td><td>0</td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>30</sub></td><td>RapGapActive</td><td>1.2 × 10<sup>5</sup></td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>31</sub></td><td>PP2AActive</td><td>1.2 × 10<sup>5</sup></td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr><tr><td>x<sub>32</sub></td><td>Raf1PPtase</td><td>1.2 × 10<sup>5</sup></td><td>[0, 1.32 × 10<sup>5</sup>]</td></tr></tbody></table>
<table><caption>**Table 4.** Initial concentrations, names and bounds for parameters in the EGF-NGF signaling pathway model [5, 27, 29]. Known parameters with no variation are indicated as constant.</caption><thead><tr><th>Parameter</th><th>Name</th><th>Initial value</th><th>Range</th></tr></thead><tbody><tr><td>k<sub>1</sub></td><td>krbEGF</td><td>2.19 × 10<sup>-5</sup></td><td>[0, 4.37 × 10<sup>-5</sup>]</td></tr><tr><td>k<sub>2</sub></td><td>kruEGF</td><td>1.21 × 10<sup>-2</sup></td><td>[0, 2.42 × 10<sup>-2</sup>]</td></tr><tr><td>k<sub>3</sub></td><td>krbNGF</td><td>1.38 × 10<sup>-7</sup></td><td>[0, 2.76 × 10<sup>-7</sup>]</td></tr><tr><td>k<sub>4</sub></td><td>kruNGF</td><td>7.24 × 10<sup>-3</sup></td><td>[0, 1.45 × 10<sup>-2</sup>]</td></tr><tr><td>k<sub>5</sub></td><td>kEGF</td><td>6.95 × 10<sup>2</sup></td><td>constant</td></tr><tr><td>k<sub>6</sub></td><td>KmEGF</td><td>6.09 × 10<sup>6</sup></td><td>constant</td></tr><tr><td>k<sub>7</sub></td><td>kNGF</td><td>3.89 × 10<sup>2</sup></td><td>constant</td></tr><tr><td>k<sub>8</sub></td><td>KmNGF</td><td>2.11 × 10<sup>3</sup></td><td>constant</td></tr><tr><td>k<sub>9</sub></td><td>kdSos</td><td>1.61 × 10<sup>3</sup></td><td>constant</td></tr><tr><td>k<sub>10</sub></td><td>KmdSos</td><td>8.97 × 10<sup>5</sup></td><td>constant</td></tr><tr><td>k<sub>11</sub></td><td>kSos</td><td>3.23 × 10<sup>1</sup></td><td>[0, 6.47 × 10<sup>1</sup>]</td></tr><tr><td>k<sub>12</sub></td><td>KmSos</td><td>3.59 × 10<sup>4</sup></td><td>[0, 7.19 × 10<sup>4</sup>]</td></tr><tr><td>k<sub>13</sub></td><td>kRasGap</td><td>1.51 × 10<sup>3</sup></td><td>constant</td></tr><tr><td>k<sub>14</sub></td><td>KmRasGap</td><td>1.43 × 10<sup>6</sup></td><td>constant</td></tr><tr><td>k<sub>15</sub></td><td>kRasToRaf1</td><td>8.84 × 10<sup>-1</sup></td><td>[0, 1.77]</td></tr><tr><td>k<sub>16</sub></td><td>KmRasToRaf1</td><td>6.25 × 10<sup>4</sup></td><td>constant</td></tr><tr><td>k<sub>17</sub></td><td>kpRaf1</td><td>1.86 × 10<sup>2</sup></td><td>[0, 3.72 × 10<sup>2</sup>]</td></tr><tr><td>k<sub>18</sub></td><td>KmpRaf1</td><td>4.77 × 10<sup>6</sup></td><td>constant</td></tr><tr><td>k<sub>19</sub></td><td>kpBRaf</td><td>1.25 × 10<sup>2</sup></td><td>constant</td></tr><tr><td>k<sub>20</sub></td><td>KmpBRaf</td><td>1.58 × 10<sup>5</sup></td><td>constant</td></tr><tr><td>k<sub>21</sub></td><td>kdMek</td><td>2.83</td><td>constant</td></tr><tr><td>k<sub>22</sub></td><td>KmdMek</td><td>5.19 × 10<sup>5</sup></td><td>constant</td></tr><tr><td>k<sub>23</sub></td><td>kpMekCytoplasmic</td><td>9.85</td><td>[0, 1.97 × 10<sup>1</sup>]</td></tr><tr><td>k<sub>24</sub></td><td>KmpMekCytoplasmic</td><td>1.01 × 10<sup>6</sup></td><td>constant</td></tr><tr><td>k<sub>25</sub></td><td>kdErk</td><td>8.89</td><td>constant</td></tr><tr><td>k<sub>26</sub></td><td>KmdErk</td><td>3.49 × 10<sup>6</sup></td><td>constant</td></tr><tr><td>k<sub>27</sub></td><td>kpP90Rsk</td><td>2.14 × 10<sup>-2</sup></td><td>[0, 4.27 × 10<sup>-2</sup>]</td></tr><tr><td>k<sub>28</sub></td><td>KmpP90Rsk</td><td>7.64 × 10<sup>5</sup></td><td>[0, 1.53 × 10<sup>6</sup>]</td></tr><tr><td>k<sub>29</sub></td><td>kPI3K</td><td>1.07 × 10<sup>1</sup></td><td>[0, 2.13 × 10<sup>1</sup>]</td></tr><tr><td>k<sub>30</sub></td><td>KmPI3K</td><td>1.85 × 10<sup>5</sup></td><td>constant</td></tr><tr><td>k<sub>31</sub></td><td>kPI3KRas</td><td>7.71 × 10<sup>-2</sup></td><td>constant</td></tr><tr><td>k<sub>32</sub></td><td>KmPI3KRas</td><td>2.72 × 10<sup>5</sup></td><td>constant</td></tr><tr><td>k<sub>33</sub></td><td>kAkt</td><td>5.66 × 10<sup>-2</sup></td><td>[0, 1.13 × 10<sup>-1</sup>]</td></tr><tr><td>k<sub>34</sub></td><td>KmAkt</td><td>6.54 × 10<sup>5</sup></td><td>[0, 1.31 × 10<sup>6</sup>]</td></tr><tr><td>k<sub>35</sub></td><td>kdRaf1ByAkt</td><td>1.51 × 10<sup>1</sup></td><td>constant</td></tr><tr><td>k<sub>36</sub></td><td>KmRaf1ByAkt</td><td>1.19 × 10<sup>5</sup></td><td>constant</td></tr><tr><td>k<sub>37</sub></td><td>kC3GNGF</td><td>1.47 × 10<sup>2</sup></td><td>[0, 2.94 × 10<sup>2</sup>]</td></tr><tr><td>k<sub>38</sub></td><td>KmC3GNGF</td><td>1.29 × 10<sup>4</sup></td><td>[0, 2.58 × 10<sup>4</sup>]</td></tr><tr><td>k<sub>39</sub></td><td>kC3G</td><td>1.40145</td><td>[0, 2.8029]</td></tr><tr><td>k<sub>40</sub></td><td>KmC3G</td><td>1.10 × 10<sup>4</sup></td><td>[0, 2.19 × 10<sup>4</sup>]</td></tr><tr><td>k<sub>41</sub></td><td>kRapGap</td><td>2.73 × 10<sup>1</sup></td><td>constant</td></tr><tr><td>k<sub>42</sub></td><td>KmRapGap</td><td>2.96 × 10<sup>5</sup></td><td>constant</td></tr><tr><td>k<sub>43</sub></td><td>kRap1ToBRaf</td><td>2.21</td><td>[0, 4.42]</td></tr><tr><td>k<sub>44</sub></td><td>KmRap1ToBRaf</td><td>1.03 × 10<sup>6</sup></td><td>[0, 2.05 × 10<sup>6</sup>]</td></tr><tr><td>k<sub>45</sub></td><td>kdRaf1</td><td>1.26 × 10<sup>-1</sup></td><td>constant</td></tr><tr><td>k<sub>46</sub></td><td>KmdRaf1</td><td>1.06 × 10<sup>3</sup></td><td>constant</td></tr><tr><td>k<sub>47</sub></td><td>kdBRaf</td><td>4.41 × 10<sup>2</sup></td><td>constant</td></tr><tr><td>k<sub>48</sub></td><td>KmdBRaf</td><td>1.08 × 10<sup>7</sup></td><td>constant</td></tr></tbody></table>