PCA-based synthetic sensitivity coefficients for chemical reaction network in cancer

Background: computation of CRNs equilibriumWe consider a CRN comprising r reactions involving n well-mixed chemical species and assume that the law of mass action kinetics holds. Throughout the paper, we shall denote with \(x_i(t)\), \(i=1, \dots , n\), the molar concentration (nM) at time t of the \(i-\)th species \(A_i\), and with \(k_j\), \(j=1, \dots , r\), the rate constant corresponding to the \(j-\)th reaction. The dynamic of the species concentration can be described by a system of n ordinary differential equations (ODEs) as in Eq. (1).We further assume the CRN to be weakly elemented4, that is \(p=n-{{\,\textrm{rank}}}({\varvec{S}})\) independent conservation laws exist described by the linear system$$\begin{aligned} {\varvec{N}}{\varvec{x}}= {\varvec{c}}, \end{aligned}$$
(3)
where \({\varvec{c}}\in R^p_+\) is the vector of conservations laws’ constants whose values can be defined as \({\varvec{c}}:= {\varvec{N}}{\varvec{x}}_0\), being \({\varvec{x}}_0\) the initial values of the species concentrations , and the matrix \({\varvec{N}}\in N^{p \times n}\) is determined by computing the kernel of \({\varvec{S}}^T\) and is assumed to contain a minor equal to the identity matrix of size p. Therefore, up to a change in the order of the components of \({\varvec{x}}\), we can consider the decompositions$$\begin{aligned} {\varvec{N}}= \left[ \begin{array}{cc} {\varvec{I}}_p&{\varvec{N}}_2^{p \times (n-p)} \end{array}\right] , \quad {\varvec{x}}= \left[ \begin{array}{c} {\varvec{x}}_1 \\ {\varvec{x}}_2 \end{array}\right] , \quad {\varvec{S}}= \left[ \begin{array}{c} {\varvec{S}}_1^{p \times r} \\ {\varvec{S}}_2^{(n-p) \times r} \end{array} \right] \ , \end{aligned}$$
(4)
where \({\varvec{x}}_1\) is a p-vector of the so called elemental species4, and \({\varvec{S}}_2\in R^{(n-p) \times r}\) has rank \((n-p)\) and is the submatrix of \({\varvec{S}}\) corresponding to the non-elemental species. An illustrative example on how to compute and interpret the conservation laws and the elemental species for a simpler network is provided in the supplementary material (see Supplementary Note S1).It follows from the definition of conservation laws that \({\varvec{N}}{\varvec{S}}=0\), whence \({\varvec{S}}_1 = – {\varvec{N}}_2 {\varvec{S}}_2\). It can be shown13,26 that the equilibrium (stationary) solution of the Cauchy problem is determined by the algebraic system$$\begin{aligned} \left\{ \begin{array}{c} {\varvec{S}}_2 {\varvec{v}}({\varvec{x}},{\varvec{k}})=0 \\ {\varvec{N}}{\varvec{x}}- {\varvec{c}}=0 \end{array} \right. \, . \end{aligned}$$
(5)
We observe that seeking the solution of this algebraic system is equivalent to the root-finding problem$$\begin{aligned} {\varvec{f}}({\varvec{x}},{\varvec{k}},{\varvec{c}}) =0, \end{aligned}$$
(6)
where \(f: R^n \times R^r \times R^p \longrightarrow R^n\) is a continuously differentiable function defined as$$\begin{aligned} {\varvec{f}}({\varvec{x}},{\varvec{k}},{\varvec{c}})= \left[ \begin{array}{c} {\varvec{S}}_2 {\varvec{v}}({\varvec{x}},{\varvec{k}})\\ {\varvec{N}}{\varvec{x}}- {\varvec{c}}\end{array}\right] \, . \end{aligned}$$
(7)
Since each reaction involves up to two reactants, from the law of mass action it follows that the components of \({\varvec{f}}\) are polynomials of degree two or lower.A CRN is said to satisfy the global stability condition4 if to any fixed \({\varvec{k}}^*\) and \({\varvec{c}}^*\) there corresponds a unique asymptotically stable equilibrium \({\varvec{x}}^*\), which may be determined by either a dynamic simulation4 or the solution of (5)13.It has been demonstrated that the action of two particular classes of mutations can be simulated by modifying the original Cauchy problem associated to the CRN through proper projectors acting on the initial conditions, and thus on the conservation laws’ constants \({\varvec{c}}\), or on the stoichiometric matrix \({\varvec{S}}\)4,5. In detail, a class of mutations resulting in the loss of function (LoF) of an elemental species can be implemented by setting to zero the constant of the conservation law where it is involved. It follows that the concentration of all the species involved in such a conservation law will remain null over time. A class of mutations resulting in the gain of function (GoF) of a species is implemented by removing from the network the reactions involved in its deactivation.Local sensitivity matrix and statistical sensitivity indicesIn this section we thoroughly describe the two steps that form the proposed algorithm for computing the SSIs.Step 1. Analytic computation of the sensitivity matrixWe consider a weakly elemented CRN and we assume that$$\begin{aligned} \det [(\nabla _{\varvec{x}}{\varvec{f}})({\varvec{x}}^*,{\varvec{h}}^*)] \ne 0\ \ . \end{aligned}$$
(8)
On account of the implicit function theorem47, equations (5) and (8) imply that there exists one and only one \({\varvec{x}}={\varvec{x}}({\varvec{h}})\) such that$$\begin{aligned} {\varvec{f}}({\varvec{x}}({\varvec{h}}),{\varvec{h}}) = 0 \end{aligned}$$
(9)
for \({\varvec{h}}\) in a neighborhood of \({\varvec{h}}^*\), meaning that \({\varvec{x}}({\varvec{h}})\) is the equilibrium point of the CRN defined by the set of parameters \({\varvec{h}}\). Further,$$\begin{aligned} \nabla _{{\varvec{h}}} {\varvec{x}}({\varvec{h}}^*)= – \big [ \nabla _{{\varvec{x}}} {\varvec{f}}({\varvec{x}}^*, {\varvec{h}}^*) \big ]^{-1} \, \nabla _{{\varvec{h}}} {\varvec{f}}({\varvec{x}}^*, {\varvec{h}}^*) \, , \end{aligned}$$
(10)
where, according to the definition of \({\varvec{f}}\) in (7),$$\begin{aligned} \nabla _{{\varvec{x}}} {\varvec{f}}({\varvec{x}}^*, {\varvec{h}}^*) = \left[ \begin{array}{c} {\varvec{S}}_2 \nabla _{\varvec{x}}{\varvec{v}}({\varvec{x}}^*, {\varvec{k}}^*) \\ {\varvec{N}}\end{array} \right] \end{aligned}$$
(11)
and$$\begin{aligned} \nabla _{{\varvec{h}}} {\varvec{f}}({\varvec{x}}^*, {\varvec{h}}^*) = \left[ \begin{array}{cc} {\varvec{S}}_2 \nabla _{\varvec{k}}{\varvec{v}}({\varvec{x}}^*, {\varvec{k}}^*) &{} \textbf{0}\\ \textbf{0} &{} – \textbf{I}_p\end{array} \right] \, . \end{aligned}$$
(12)
Step 2. Computation of the SSIFollowing the work by Liu and colleagues23, to quantify the overall effect that a relative change in the parameters \({\varvec{h}}\) induces on the protein concentrations at equilibrium, we define the positive semi-definite quadratic form$$\begin{aligned} Q(\hat{\Delta } {\varvec{h}}) = \hat{\Delta } {\varvec{x}}^T \hat{\Delta } {\varvec{x}}\end{aligned}$$
(13)
where$$\begin{aligned} \hat{\Delta } {\varvec{h}}= \left[ \frac{h_1 – h_1^*}{h_1^*}, \cdots , \frac{h_{r+p} – h_{r+p}^*}{h_{r+p}^*} \right] ^T \end{aligned}$$
(14)
and$$\begin{aligned} \hat{\Delta } {\varvec{x}}= \left[ \frac{x_1({\varvec{h}}) – x_1^*}{x_1^*}, \cdots , \frac{x_{n}({\varvec{h}}) – x_{n}^*}{x_{n}^*} \right] ^T \end{aligned}$$
(15)
\({\varvec{x}}({\varvec{h}}) = (x_1({\varvec{h}}), \dots , x_n({\varvec{h}}))^T\) being the equilibrium reached by the network with parameters \({\varvec{h}}\). We observe that, in order to compute \(\hat{\Delta } {\varvec{x}}\), we shall assume that \(x_{i}^* \ne 0\) for all \(i=1, \dots , n\), which is true, e.g., for the original CR-CRN. However, as we shall see in the next sections, if some of the species have a null equilibrium concentration our approach can be applied by simply removing the components of \({\varvec{x}}\) and the rows of the sensitivity matrix corresponding to such species. By applying some straightforward algebraic computation to the first order Taylor’s polynomial of the implicit function \({\varvec{x}}= {\varvec{x}}({\varvec{h}})\), it follows that in a neighbour of \({\varvec{h}}^*\) it holds$$\begin{aligned} \hat{\Delta } {\varvec{x}}= \varvec{\mathcal {S}} \, \hat{\Delta } {\varvec{h}}\ \ , \end{aligned}$$
(16)
where we have introduced the relative local sensitivity matrix29$$\begin{aligned} \varvec{\mathcal {S}} = ({\varvec{D}}_{{\varvec{x}}^*})^{-1} \, \nabla _{{\varvec{h}}} {\varvec{x}}({\varvec{h}}^*) \, {\varvec{D}}_{{\varvec{h}}^*} \ , \end{aligned}$$
(17)
being \({\varvec{D}}_{{\varvec{x}}^*} = \textrm{diag}[ x_1^*, …, \, x_n^*]\) and \({\varvec{D}}_{{\varvec{h}}^*} = \textrm{diag} [ h_1^*, …, \, h_{r+p}^*]\). Replacing equation (16) into (13) leads to$$\begin{aligned} Q(\hat{\Delta } {\varvec{h}}) = \hat{\Delta } {\varvec{h}}^T \varvec{\mathcal {S}}^T \varvec{\mathcal {S}} \hat{\Delta } {\varvec{h}}= \hat{\Delta } {\varvec{h}}^T\, {\varvec{U}}\, {\varvec{\Lambda }}\, {\varvec{U}}^T \, \hat{\Delta } {\varvec{h}}\ , \end{aligned}$$
(18)
where \({\varvec{\Lambda }}= \textrm{diag} [\lambda _1,…,\,\lambda _{r+p}]\), \(\lambda _1 \ge \lambda _2 \ge … \ge \lambda _{r+p} \ge 0\), are the eigenvalues of the matrix \(\varvec{\mathcal {S}}^T \varvec{\mathcal {S}}\) and \({\varvec{U}}= [{\varvec{u}}_1, …,\,{\varvec{u}}_{r+p}]\) collects a set of independent normalized eigenvectors.As a special case, let us assume that only the j-th parameter is perturbed and thus \(\hat{\Delta } {\varvec{h}}\) has \(\hat{\Delta } h_j\) as a unique non-vanishing component. Then, equation (18) reads as$$\begin{aligned} Q(\hat{\Delta } h_j) = \hat{\Delta } h_j^2\sum _{a=1}^{r+p} \lambda _a \, u_{ja}^2 . \end{aligned}$$
(19)
This equation inspires the definition of the j-th statistical sensitivity index (SSI)$$\begin{aligned} e_j^{{\varvec{h}}} = \frac{ \textstyle \sum \limits _{a=1}^{r+p} \lambda _a \, u_{ja}^2 }{\textstyle \sum \limits _{a=1}^{r+p} \lambda _a} \quad j = 1, \dots r+p \ . \end{aligned}$$
(20)
It can be shown that \(e_j^{{\varvec{h}}} \in [0, 1]\) and \(\sum _{j=1}^{r+p} e_j^{{\varvec{h}}} =1\). Further, getting back to the original parameters \({\varvec{k}}\) and \({\varvec{c}}\), the first r SSIs in (20), also denoted as \(e_j^{{\varvec{k}}}\), \(j = 1, \dots , r\), provide a synthetic and compact estimate of the weight of perturbation in reaction j on the whole CRN, while the latter p SSIs in (20), also denoted as \(e_j^{{\varvec{c}}}\), \(j=1, \dots ,p\), refers to the overall sensitivity of the network to the j-th component of vector \({\varvec{c}}\).We observe that the definition of the SSIs in Eq. (20) is similar to that of Liu and colleagues23 , the only difference being that we used the squared version, \(u_{ja}^2\), of the eigenvector components instead of simply using \(u_{ja}\). The motivation of this choice is twofold. On the one hand, it ensures that all the terms of the sums in (20) are positive thus avoiding cancellation effects. On the other hand, it accounts for the fact that the sign of the eigenvectors is ambiguous. In fact, for any eigenvector \({\varvec{u}}_j\), \(j=1, \dots , r+p\), also \(-{\varvec{u}}_j\) is an eigenvector of \(\varvec{\mathcal {S}}^T \varvec{\mathcal {S}}\) and different Matlab routines may actually return eigenvectors with different sign (see Supplementary Note S2 for further details) .Sensitivity analysis of subnetworksSince the SSIs in (20) are defined based on a PCA of the whole relative local sensitivity matrix \(\varvec{\mathcal {S}}\), their values allow discriminating among all the parameters whose variation mainly affects the overall equilibrium point reached by the network. However, it is possible to focus only on a subset of chemical species \(\mathcal {A}=\{A_{i_1},…, A_{i_{\tilde{n}}}\}\), \(\tilde{n}<n\), and on a subset of parameters \(\widetilde{{\varvec{h}}} = \left( h_{j_1}, \dots , h_{j_\tau } \right)\), \(\tau < r+p\), by applying the procedure just described on the submatrix \(\widetilde{\varvec{\mathcal {S}}}\) of \(\varvec{\mathcal {S}}\), formed by the \(\tilde{n}\) rows referred to \(\mathcal {A}\) and the \(\tau\) columns referred to \(\widetilde{{\varvec{h}}}\). \(\widetilde{\varvec{\mathcal {S}}}\) represents the relative local sensitivity matrix of the equilibrium concentrations \((x_{i_1}, \dots , x_{i_{\tilde{n}}})\) of the protein in \(\mathcal {A}\) with respect to the parameters in \(\widetilde{{\varvec{h}}}\). Hence the SSIs computed by using in (20) the eigenvalues and eigenvectors of \(\widetilde{\varvec{\mathcal {S}}}^T\widetilde{\varvec{\mathcal {S}}}\) quantify the impact of each parameters only on the considered proteins.As an example, let use assume \(\widetilde{{\varvec{h}}}={\varvec{h}}\), \(\tilde{n} = 1\), and \(\mathcal {A}=\{A_{i}\}\). Then \(\widetilde{\varvec{\mathcal {S}}}\) is the i-th row of \(\varvec{\mathcal {S}}\), denoted as \(\varvec{\mathcal {S}}_{i:}\), and the matrix \(\widetilde{\varvec{\mathcal {S}}}^T\widetilde{\varvec{\mathcal {S}}}\) has a unique non-vanishing eigenvalue equal to \(||\varvec{\mathcal {S}}_{i:}||^2\) with \(\frac{\varvec{\mathcal {S}}_{i:}}{||\varvec{\mathcal {S}}_{i:}||}\) as a normalized eigenvector. From the definition of the SSIs in (20) it follows that$$\begin{aligned} e_j^{{\varvec{h}}}=\frac{\mathcal {S}_{ij}^2}{||\varvec{\mathcal {S}}_{i:}||} \, , \end{aligned}$$
(21)
i.e., the SSIs coincide with the squared normalized component of the i-th row of \(\varvec{\mathcal {S}}\). This confirm our interpretation of the SSIs. Indeed, the i-th row of the rescaled sensitivity matrix \(\varvec{\mathcal {S}}\) describes the relative sensitivity of the equilibrium concentrations of the species \(\mathcal {A}\) with respect to the parameters of the network29.Sensitivity analysis of mutated CRNsAs described in the Background section, a mutation resulting in the LoF of an elemental species can be implemented by setting to zero the constant of the conservation law involving such a species4,5. Coherently, in this scenario the algorithm sketched in Fig. 1 can be applied by setting to zero the component of \({\varvec{h}}^*\) corresponding to the modified conservation law and then computing \({\varvec{x}}^*\) as the equilibrium point of the CRN with the novel set of parameters. Furthermore, the species involved in the modified conservation law are removed from the set \(\mathcal {A}\) because their equilibrium concentration will clearly be zero.A mutation resulting in the GoF of a species implies that a set of reactions involving it no longer occurs in the network5. Hence, a reduced network is defined by removing the columns of \({\varvec{S}}_2\), the rows of \({\varvec{v}}\), and the rate constants within \({\varvec{h}}^*\) associated to these reactions. The algorithm sketched in Fig. 1 can be applied to the novel network after computing the equilibrium \({\varvec{x}}^*\) and by defining \(\mathcal {A}\) as the set of species with an equilibrium concentration greater than \(10^{-15}\) so as to neglect those species whose concentration is null or close to zero.