Docstoc

web_supp_jan19

Document Sample
web_supp_jan19 Powered By Docstoc
					      Supplementary Material for “Dramatic reduction of dimensionality in large
               biochemical networks due to strong pair correlations”




Table of Contents:                                                                       Page number

1. Supplementary Figures S1-S4 ………………………                                                 1-3
2. Generation of Random Linear Networks ……………..                                           4
3. Nonlinear Feedback Reactions in a
   Reaction Network with Linear Architecture …………                                        5
4. Biological Networks (Fig. S6) ……………………….                                              7-9
5. Supplementary Figures S7-S8 ……………………......                                            9-10
6. Details of the Simulation Method ……………………                                             11-12
7. Supplementary Figures S9 …………………………...                                                13
8. Failure of PCA in reducing dimensionality
   in a simple network …………………………………...                                                  14-15
9. Analytical Calculations for
   the Linear Network ……………………………………                                                     17-19
10. Supplement on Gramians …………………………….                                                  20-21
11. Studying the role of modifications likely to occur in exp-
    eriments .........................................................................   22-26
12. Numbers of species and timepoints used for the biological
    networks …………………………………………........                                                    27
13. Effect of unit variance scaling (Fig. S12) ...................                       28
14. References ...................................................................       29
                         N=64 Linear Cascade: Top Two Components
                       100                                              100

                                                                        99
                       80                                               98

                                                                        97
         Trial Index




                       60
                                                                        96

                                                                        94
                       40
                                                                        93

                       20                                               92

                                                                        91

                        0                                               90
                         0     0.2     0.4          0.6       0.8   1
                                             Time



Fig. S1. Time dependence of the percentage explained by the top two eigenvalues of the
covariance matrix calculated for the linear network for N=64 species. Results for100 trials (each
associated with a set of rate constants and initial conditions randomly drawn from a uniform
distribution of a range of 0.1 to 10) are shown and demonstrate that more than 90% of the total
variance is captured by the top two eigenvalues most of the time. All the trials were started with
a non-zero concentration of the last species ( c64 ! 0 ) and zero concentrations for all other
species. The calculation of the covariance matrix included all the 64 species concentrations.




                                                          1
      (A)                                              (B)
Fig. S2. Minimum percentage explained by the top and the top two eigenvalues of the
covariance matrix for the linear reaction network. (A) Distributions of the minimum of the
percentage explained by the top eigenvalue in the entire time course for any trial for N=8 to
N=512. In each trial a set of rate constants and initial concentrations were randomly drawn from
a uniform distribution between 0.1 and 10. All the 10,000 trials were started with a non-zero
concentration of the last species and zero concentrations for all other species. The calculation of
the covariance matrix included all the species concentrations. (B) Distributions of the minimum
percentage explained by the top two eigenvalues for the same cases as in (A).




(A)                             (B)                                   (C)
Fig. S3 Effect of initial conditions on the distributions of the minimum percentage
explained by the eigenvalues in the linear network. We compare the distributions of the
minimum percentage explained by the largest eigenvalue of the covariance matrix for two
different types of initial conditions for the linear network with N=64. The two different types
correspond to (i) only the concentration of the last species is non-zero ( c64 ! 0 ) at t=0, whereas,
the rest of the species start with zero concentrations, (ii) all the species concentrations have non-
zero values at t=0. For both the cases, the rate constants and the initial concentrations were
varied in the same manner (as in Fig. S2). (A) Shows distributions of the largest eigenvalue. (B)
Distributions for the top two eigenvalues. (C) Distributions for the top three eigenvalues. The top



                                                  2
eigenvalue seem to capture the variances more efficiently when all the species were given non-
zero initial concentrations . The shapes of the distributions start to converge when three
eigenvalues are considered.




    (A)                                                  (B)
Fig. S4. Comparisons between results of calculations including all the species vs. including
only a set of independent species in the covariance matrix for the linear network. (A)
Distributions of the minimum percentage explained by the largest eigenvalue when the
covariance matrix was calculated using all the species (black line) and only the independent
species concentrations (red line) for N=64. The number of the independent species was
calculated by accounting for the conservation laws (e.g., total concentrations of the all the
                                                   N
species remains unchanged in the linear network,   ! c (t) = const ) that produce linear
                                                   i=1
                                                         i


dependence between different species concentrations. The parameters used for the simulation are
the same as described in Fig. S2. (B) Distributions of the minimum percentage explained by the
top three eigenvalues for the same cases as in (A).






                                               3
Generation of Random Linear Networks
1. For a system of N number of species, we choose a subset of n1 ! N species, where n1 is drawn
from a uniform random number distribution ranged from 2 to N. The species in the sub-network,
indexed as, c1,!, cn1 , are connected such that ci reacts with ci+1 with first order forward and
backward reactions for i < n1 . The step produces a connected reaction network with linear
architecture. Then in addition to these reactions, we choose all the pairs that are not directly
connected to each other and introduce a pair of forward and backward first order reaction, each
with probability ½. The second step produces branched structure in the sub-network.


2. The next step involves three sub-steps. A. If, N ! n1 > 2 , then, we choose another subnet-work
of n2 ! N " n1 following the same method as described in step 1. The reactions in the sub-
network are also generated in the same manner. We continue this process as long as, a single
species or a pair of species is left in the system. We follow the next two steps when that situation
is reached.
B. If, N ! n1 = 2 , then we connect the remaining pair of species by first order forward and
backward reactions.
C. If N ! n1 = 1 , then we connect the last species with a randomly chosen species in the last sub-
network by first order forward and backward reactions.


After we generate the network architecture that is consists of a collection of disconnected sub-
networks with branched and/or linearly connected reactions we assign rate constants and initial
conditions to all the reactions from the range [0.1-10] for the rates and [0.1-10] for the initial
concentrations.


At each trial, a new network architecture and a new set of reaction rates and initial conditions are
generated.






                                                 4
Nonlinear feedback reactions in a reaction network with linear architecture




The species in the network react with each other as shown above in the schematic figure. The
wavy lines indicate the feedback reactions. In the above network, species i, interacts with species
i and i+1 as following the mass-action kinetics below,
    dci
        = ki!1ci!1 + k!i ci+1 ! (k!(i!1) + ki )ci ! ki!1ci ci!1 + kif ci ci+1
                                                      f

    dt

The last two terms in the above ODE denote the feedback reactions.




Fig. S5. Minimum percentage explained by the top four eigenvalues in the linear reaction
networks with non-linear feedback reactions. Distributions of the minimum percentage
explained shown for N=8 to N=512 for 10,000 trials. In each trial a set of rate constants and
initial concentrations were randomly drawn from a uniform distribution in the range of 0.1 to
10). Each trial was started with a non-zero concentration for the last species and zero
concentrations for all other species.






                                                                     5
    (A) Ras Activation Network




(B) EGFR Signaling Network (Kholodenko et al. (1999))






                                          6
(C) EGFR Signaling Network (Blinov et al. (2006))






                                           7
(D) NF-kB Activation Network


Fig. S6. Biological networks. (A) Ras activation network. Ras activation (transition from
RasGDP to RasGTP) is carried out by two enzymes (Guanine nucleotide Exchange Factors or
GEFs) SOS and Rasgrp1. SOS is recruited to the plasma membrane by the T cell receptor (TCR)
signalosome (not shown in the figure) which is formed when TCRs interact with strong affinity
antigens. In addition to a catalytic site (marked by a yellow circle) that mediates Ras activation,
SOS also has an allosteric site which can bind to both RasGDP and RasGTP (1). When the
allosteric site is occupied by RasGTP the rate of Ras activation via SOS becomes almost 75
times faster. This is the positive feedback. Rasgrp1 is recruited to the plasma membrane by
Diacylglycerol (DAG) which is produced by signaling events downstream to TCR triggering.
Unlike SOS, Ras activation by Rasgrp1 does not have a positive feedback. Ras deactivation is
mediated by the enzyme, RasGAP. (B) EFGR signaling network (coarse grained). EGFR
signaling network where multiple receptor activation states are coarse-grained into a single state.
The details regarding the signaling pathway shown here can be found at
http://models.cellml.org/exposure/468382d300fa0027425f72c195afd11a
and Ref. (2). (C)
EFGR signaling network (including detailed receptor activation states). Details of the
signaling network can be found in Ref. (3). The figure above was taken from Ref. (3) with the
permission of the authors. (D) NF-kB activation network. In a resting state NF-κB is bound to
the three isoforms of IκB (α, β, ε) in cytoplasm. Following stimulation, IKK-complexes get




                                                 8
activated. IKK phosphorylates and degrades IκB isoforms, releasing NF-κB from its bound state.
Free NF-κB translocates to the nucleus(4). While the transcription factors for the IκBβ and IκBε
are produced at a steady rate, the transcription factor for the IκBα gets positively regulated by
nuclear NF-κB. The nuclear IκB isoforms then bind with their nuclear NF-κB partners and
translocate them back to cytoplasm, where in presence of IKK, the same cycle repeats itself,
giving rise to oscillation in NF-κB’s nuclear localization. The BioNetGen files describing the
reactions and the parameters used can be found at the link
http://planetx.nationwidechildrens.org/~jayajit/pairwise_correlations/

.





                    Ras−Sos Network: Top Component                                       Ras−Sos Network: Top Two Components
              100                                                                  100                                             100
                                                           100
                                                                                                                                   99
                                                           95
                                                                                   80                                              98
              80                                           89
                                                                                                                                   97
                                                           84
                                                                     Trial Index




                                                                                   60
Trial Index




              60                                                                                                                   96
                                                           78
                                                                                                                                   94
                                                           73                      40
              40                                                                                                                   93
                                                           68
                                                                                   20                                              92
              20                                           62
                                                                                                                                   91
                                                           57
                                                                                    0                                              90
               0                                           52                        0        0.2    0.4          0.6   0.8    1
                0     0.2     0.4          0.6   0.8   1                                                   Time
                                    Time




Fig. S7 (A) Temporal evolution of the percentage explained by the largest eigenvalue for the
Ras activation network until the steady state is reached for 100 different trials (denoted as trial
index in the y axis). In each trail, the kinetic rates and initial concentrations were varied around
the base parameters (web supplement) by up to 50% from a uniform random distribution. (B)
Time evolution of the percentage explained by the top two eigenvalues for the same trails as in
(A).






                                                                 9
                              Minimum Percent Explained by Top Two Components
                                                Ras-Sos Network
                        100
                                  Stable
                                  Bistable
                        80
    Percent of Trials




                        60


                        40


                        20


                         0
                         20       30    40    50     60       70        80   90   100
                                             Percent Explained




(A)                                                                                          (B)




                                                                  (C)



Fig. S8 Percentage explained by more than one principal components for the nonlinear
networks. (A) Shows the Ras activation network, (B) EGFR signaling and (C) NF-kB activation
networks. The parameters used are the same as in Fig. 2 in the main text.






                                                                                        10
Details of the Simulation Method:
We use the rule based modeling software, BioNetGen(5), to generate time courses for the species
kinetics in a biochemical network. This software produces a set of ODEs corresponding to the
mass-action kinetics describing bio-chemical reactions in the networks. The ODEs are solved
numerically in the software using the CVODE solver. We have used a time step ( dt )
proportional to the inverse of the fastest time scale ( dt = ! fastest / 2 ) associated with the reactions.
For the linear networks we calculate ! fastest from the largest rate constant, ie., ! fastest = 1 / klargest ,
for each trial as all the reactions are first order reactions. For the networks with non-linear
reactions, which consist of second order binding unbinding and first order reactions, we calculate
the largest possible rate ( klargest ) of the reactions for every trial by finding the maximum of the
first order rates and the product of the binding rates with the largest species concentration.
Following the same scheme we calculate the minimum possible rate ( ksmallest ) of a reaction in the
linear ( ksmallest =minimum of all the rates) and the non-linear networks ( ksmallest =minimum of all
the first order rates and the product of binding rates and the largest species concentration). The
BIONETGEN files describing the reactions for all the networks studied can be found at the link
http://planetx.nationwidechildrens.org/~jayajit/pairwise_correlations/. We have checked in
specific cases that the results obtained from BionetGen agree with results from higher accuracy
calculations using Mathematica.


Calculation of the steady state: We calculated the eigenvalues of the correlation matrix, Cij ,
for times until the system reached the steady state. When reaction kinetics is described by ODEs,
the system reaches the steady state in principle at t ! " . However, as the system approaches the
fixed point the relative changes in species concentrations compared to the steady state
concentrations decreases with time. In our simulations, we ended the kinetics when we reached
such a limit. We denote this time as the estimated steady state time (tS ) . For all the networks
(linear and non-linear, except the NF-κB network) tS is calculated from the smallest time scale
(! ksmallest ) in the system. We justified this choice by evolving the system for a longer time to
                                                                  N                         N
check if the relative Euclidean distance, i.e., d(t1, t2 ) =     " (ci (t1 ) ! ci (t2 ))2
                                                                 i=1
                                                                                            " c (t = 0) between
                                                                                            i=1
                                                                                                  i


the species concentrations is smaller than a tolerance of 10-6 when t1 = tS and t2 = 2tS . The table
below (Table I) shows the values of these variables for the different networks we use. For some
of the networks, the system reached close to steady state concentrations even for times t < tS , in
such cases, we calculated the Euclidean distance, d(tend , tS ) and end the simulation when
d(tend , tS ) ! 10 "6 . The table below also provides details regarding this.



NF-κB network: This network can produce oscillations when the kinetic rates and the initial
conditions fall within a particular range. We stop the simulations when 50 oscillations are
detected in a trial. Our simulations indicated it was reasonable to assume that the system is well




                                                       11
inside the limit cycle after 50 oscillations. For the kinetic trajectories that did not show
oscillations we stopped our kinetics at a time t=200 units which guaranteed that the distributions
for the minimum percentage explained did not change significantly even when the system was
evolve farther in time. We used a less stringent criterion for estimating the steady states for this
network compared to the other networks since this case required a much smaller discretization
time step ( dt = 4 !10 "4 ) which substantially increased the computation time for the range of
parameters studied here.


Table I: Details regarding the solutions of the ODEs and steady state calculations



                          Network            Estimated            Simulation end
                                            steady state            time (tend )
                                              time (tS )
                       Linear cascade        10/ ksmallest          Same as tS
                                          if max{ci } < 1 .
                                              When,
                                            max{ci } ! 1
                                          (10 / ksmallest )cmax
                       Random linear          1/ ksmallest          Same as tS
                          cascade
                         Nonlinear           10/ ksmallest          Same as tS
                         feedback
                                          if max{ci } < 1 .
                                              When,
                                            max{ci } ! 1
                                          (10 / ksmallest )cmax
                         Ras-SOS              1/ ksmallest          Same as tS
                       Kholodenko et          2/ ksmallest        tend < 2 / ksmallest
                            al.
                        Blinov et al.         2/ ksmallest        tend < 2 / ksmallest






                                                   12
Fig. S9. Comparison with the Latin Hypercube Sampling method. The linear network for
N=64 was investigated with values of the rate constants and initial concentrations generated
using the Latin Hypercube(6) sampling method (red histograms) for 10,000 trials. The
parameters were varied in the same manner as described in Fig. S2. The histograms generated by
using the built-in random number generator in Perl are shown (black) for comparison. Both
sampling with a uniform distribution and Latin hypercube sampling generate very similar
distributions.






                                              13
Failure of Principal Component Analysis in reducing dimensionality in a simple linear
network


Consider a set of disconnected first order reactions shown below:



       "          "                "
X1 !k! #, X 2 !k! #,$ $ $ $ Xi !k! #,$ $ $ $ X m !k! # .
     1          2                i                 m
                                                     "                                         (1)

The kinetics of the above reactions is determined by a set of uncoupled linear ordinary rate
equations given by,


    dc1            dc             dc
        = !k1c1 ,!, i = !ki ci ,!, m = !km cm
    dt             dt              dt

If ci 0 (ti ) is the initial concentration of the ith species ( Xi ) that starts decaying at, t = ti , then the
concentration of the species at any time, t(> ti ) is given by, ci (t) = ci 0 (ti )e! ki (t !ti ) .



Case I
First we consider a situation where all the reactions start at the same time, t1 = t 2 = ! = t m = 0 ,
with the same values of initial concentrations, i.e., c10 (0) = ! = cm 0 (0) = c0 . We also assume that
all the species decay at the same rate, therefore, k1 = k2 = ! = km = 1 . With the above
simplifications, all the species follow identical kinetics, where, ci (t) = c0 e!t . Thus, the correlation
matrix, Cij , defined as,

                T
          1
Cij (T ) = ! dt ci (t)ci (t)
          T 0
               2 T
              c0             c2
              T !
          =        dt e"2t = 0 (1 " e"2T )
                 0
                            2T



Note that all the elements of the Cij matrix are equal, and there is only one eigenvalue (
        2
     mc0
!1 =      (1 " e"2T ) ) that is non-zero, and, rest of the eigenvalues vanish, i.e.,
     2T
!2 = !3 = ! = !m = 0 . Therefore, the largest eigenvalue captures 100% of the variance! This
behavior is expected, because, all the species started decaying at the same time with the same



                                                       14
rate (or “velocity”), therefore, the trajectory of the system in the manifold spanned by the species
concentrations is a straight line and an eigenvector along that direction will capture all the
variations in the trajectory at all times. Now, we will construct a case where we will introduce
each of the reactions at different times and show that a particular arrangement of those times will
create a trajectory that will require all the eigenvectors capture equal amount of variation in the
trajectory.


Case II
Here we consider a situation where the reactions in (1) are introduced sequentially at different
times ( {ti } ), such that, t m > t m !1 > t m ! 2 > ! > t 2 > t1 = 0 . As in case I, we assume all the rates are
set to unity and all the species start with the same initial concentration, c0 . The correlation
matrix, Cij , is given by,

                     T                     2 T                            2
                     1                    c0                             c0 (ti +t j ) "2ti
Cij (T , ti , t j ) = ! dt ci (t)c j (t) = ! dt e "(t "ti ) "(t "t j )
                                                           e           =    e         # e " e"2T %
                                                                                      $          &
                     T 0                  T ti                           2T

where, ti > t j and ci (t) = 0 when, t < ti and c j (t) = 0 when t < t j .



It will be instructive to look at the case with two species (m=2) first and then we will generalize
the result for any number of species. In this case,
          1                             1 !(t2 !t1 )                                    1
C11 =       "1 ! e!2(T !t1 ) $ , C12 =
            #                %            "e
                                          #          ! e!2T +(t1 +t2 ) $ = C21 , C22 =
                                                                       %                  "1 ! e!2(T !t2 ) $ , where,
         2T                            2T                                              2T #                %
t 2 > t1 . The corresponding eigenvalues are,

         (C11 + C22 ) ± (C11 " C22 )2 + 4C12
                                          2

!1,2 =
                         2

The % explained by the largest component is given by

    1"    (C11 ! C22 )2 + 4C12     %
                            2

     $1 +                          ' ( 100% .
    2$        C11 + C22            '
     #                             &



When, T ! t1 ,t 2 , the above expression simplifies to,

&1       !2(t 2 !t1 )             )
( "1 + e
   #                  $ + O(e!2T )+ , 100% . Therefore, when t 2 ! t1 ! 1 , the % explained by the top
                      %
'2                                *
component is very close to 50%, and both the eigenvectors capture the total variance almost
equally.



                                                            15
Now we can easily generalize the result for any m > 2 . When, T ! t1 ,",t m ,

                       1 !(ti !t j )
Cij (T ,ti ,t j ) =      e           + O(e!2T ), ti > t j
                      2T

and,
                       1
Cii (T ,ti ,ti ) =       + O(e!2T )
                      2T

Therefore, when the times ({tt }) for introducing the reactions are chosen such that ti ! t j ! 1,
then the off-diagonal elements of the Cij matrix decreases exponentially fast compared to the
diagonal elements. Thus, in this case we get a set of degenerate eigenvalues all equal to 1/2T in
the leading order and the percentage explained by the any eigenvalue is equal to 1 / m (to the
leading order).






                                                            16
    Linear Model
    In this section we provide analytic arguments to support the numerical result pre-
sented in the text that one or two principal components are sufficient to explain almost
all the variance in an overwhelming number of cases for the linear cascade and its variants.

   In general the concentrations of the linear system satisfy the rate equations
                                 dci (t)   ∑
                                         =      Mij cj (t)                                    (1)
                                   dt        j

where the reaction matrix elements Mij are constants that are linear in the the rate con-
stants.

    The concentrations that obey the set of reactions shown in Figure 1 in the text yield
                                                                 −
the matrix elements Mi,i+1 = k−i , Mi,i−1 = ki−1 and Mii = −ki∑ k−(i−1) for 1 ≤ i ≤ n
if we set k−1 = kn = 0. The system obeys the conservation law j cj (t) is independent
of time.

    The eigenvalues and eigenvectors of M determine the time evolution of the system.
One of the eigenvalues say λ0 will be zero reflecting the conservation law and the rest of
them can be proved to be negative. Denote the eigenvalues by −λα (λα ≥ 0) and the
right- and left-eigenvectors by vα,i and uα,j respectively. The solution can be written in
the form
                                          ∑
                                          N −1
                                 ci (t) =      e−λα t Dα vα,i                          (2)
                                                α=0
                                                                                ∑
where the coefficients Dα are determined by the initial conditions: Dα =             i uα,i ci (t =
0) . The first coefficient D0 yields the sum of the concentrations. We integrate the
solutions from t = 0 up to t = τ to obtain the expression for the elements of the covariance
matrix
                ∑∑                     (                                               )
                                         1 − e−(λα +λβ )τ   (1 − e−λα τ )(1 − e−λβ τ )
     Cij (τ ) =        Dα Dβ vα,i vβ,j                    −                                    (3)
                α   β
                                           (λα + λβ )τ            λα τ × λβ τ

Only the terms in which both λα and λβ nonzero contribute in the above sums and so the
covariance matrix is effectively an (N − 1) × (N − 1) matrix.

   We write it in the form
                                       ∑∑
                          Cij (τ ) =             Dα Dβ vα,i vβ,j × Fα,β (τ )
                                       α    β
                                 (                                                   )
                                     1 − e−(λα +λβ )τ   (1 − e−λα τ )(1 − e−λβ τ )
             where Fαβ (τ ) ≡                         −                                  .    (4)
                                       (λα + λβ )τ            λα τ × λβ τ
Since Fαβ (τ ) is independent of initial conditions and depends only on the rate constants
the matrix denoted by F is easier to analyze.




                                                 17
                                                                                         k
   An extreme limit of the cascade is one with degenerate time scales given by Ci            Cj
                                                                                         k
for all pairs i ̸= j. In this case the reactant matrix M has elements Mi,j = k for i ̸= j and
Mii = −(N − 1)k. The eigenvalues are easy to determine: there is one zero eigenvalue
and N − 1 degenerate eigenvalues with value −N k. It is convenient to let the indices run
from 0 to N − 1: λ0 = 0 and λα = −N k for α = 1, · · · (N − 1). It is easy to see that
F0α = 0 for all α. So we consider Fαβ with α and β ranging from 1 to (N − 1). Since all
the λα are identical all the matrix elements of F are the same and we will denote it by f .
It is easy to determine the eigenvalues of F : it has one non-vanishing eigenvalue (N − 1)f
and the rest vanish. Since a single eigenvalue exhausts the trace it is straightforward to
show that PCA works perfectly for all initial conditions. Given that in this extreme case
PCA works exactly, for small deviations one can argue that the performance of PCA will
be very good using matrix perturbation theory.

    The other case we consider is one in which the rate constants are randomly distributed
yielding a uniform distribution for the eigenvalues of M . We have verified numerically
that the eigenvalues of M for the case displayed in Figure 1 of the text yields an ap-
proximately uniform distribution with an excess at λα close to zero for large N . We
provide a variational bound for the percentage explained by the first principal compo-
nent of F . Denote the largest eigenvalue of F by Λmax . We use the variational bound
         ∑
Λmax ≥ α,β vα Fαβ vβ for any normalized vector v.

    For a given value of N we consider a model in which the eigenvalues of λ of the reac-
tion matrix are distributed uniformly between 1/N and 1. Since F is a positive-definite
matrix, we know from the Perron-Frobenius theorem that the eigenvector corresponding
to the largest eigenvlaue has all positive elements. This suggests choosing the trial vector
⃗ = √N −1 (1, 1, 1, ·, · · · 1, 1) that is guaranteed to have a positive overlap with the exact
v        1
                                                                  ∑ ∑
eigenvector corresponding to Λmax . We have Λmax > N 1         −1    α   β Fαβ . We replace the
sum on the eigenvalues by an integral with a uniform density. Similarly we can evaluate
the trace and find a lower bound on the percentage explained by by the first principal
component. The averaging over the eigenvalue distribution can be evaluated numerically
and the results for N = 1000 are displayed in the Figure.
    In the limit N → ∞ the integrals can be evaluated exactly in terms of incomplete
Gamma functions (not displayed) and the results are similar. We emphasize that our
results are for the matrix F and not the covariance matrix itself and therefore, for special
initial conditions this bound may be violated.

   We have studied the large τ behavior as follows. Consider the matrix F for τ >>
1/λmin where λmin is the eigenvalue of M with the least magnitude. Recall that F is an
N − 1 by N − 1 matrix with a zero eigenvalue removed. Since all the exponentials are
negligible, to leading order in τ −1 it has the form of a Cauchy matrix with
                                              1    1
                                     Fα,β =             .
                                              τ λα + λβ




                                              18
                                  0.9


                                  0.8




                      explained
                                  0.7


                                  0.6


                                  0.5
                                        1   5     10            50   100   500 1000
                                                            Τ


Figure S10: Bound on the percentage explained by the first principal component as a
function of τ . We have chosen τ ≥ 1 since the largest eigenvalue is O(1).

We use the variational bound as before and obtain
                                        1
                                            ∑ ∑       1
                             Λmax      N −1  α   β λα +λβ
                                   ≥        ∑ 1
                             T rF             α 2λα

Now we write the numerator as a sum of rows and use the fact that the harmonic mean
is less than the arithmetic mean. For row α we have
                           1            1 ∑
                          ∑         ≤          (λα + λβ ) = µ + λα
                      1
                     N −1
                                1
                           β λα +λβ   N −1 β
                ∑
where T r M =     α   λα is the trace of the reaction matrix and µ ≡                   1
                                                                                      N −1
                                                                                           Tr M   we obtain
                                                   ∑    1
                                            Λmax    α λα + µ
                                                 ≥ ∑ 1
                                            T rF     α 2λα

So far this is exact; if we assume that the eigenvalues λα of the reaction matrix are
distributed between 1/N and 1 uniformly for large N we can simplify this further. We
have µ = 1/2 and the numerator can be seen to be approximately N log 3 for large N .
The denominator involves a harmonic sum and is logarithmic: is (1/2) N log N . We have

                                                Λmax   2 log 3
                                                     ≥         .
                                                T rF   log N

Thus the fraction explained by the first component is bounded at least by 1/ log N . At
N = 100, N = 500 and N = 1000 this yields a bound of 0.477, 0.3535 and 0.31. These
bounds are consistent with our numerical results. We have been unable to establish better
analytic bounds.




                                                       19
    Supplement on Gramians

    Consider a set of k vectors {⃗j } (j = 1, · · · , k) in an N -dimensional real vector space. The
                                 v
(symmetric) Gram matrix or Gramian is defined by Gij = ⃗i · ⃗j . Note that the k vectors
                                                                     v v
are linearly independent if there exist k constants ck (not all of which are zero) such that
∑k                                                        ∑n
          v
   j=1 cj ⃗ j = 0 . This corresponds to the equation         j=1 Gij cj = 0 and linear independence
occurs if and only if det G ̸= 0. Thus the Gram determinant arises naturally in determining
linear independence. The Gram determinant is non-negative and it can be shown to be the
squared (hyper)-volume of the k-dimensional parallelepiped (or parallelotope) formed by the
vectors. This is easy to demonstrate in three dimensions.

                                                                                  c
    As in the differential geometry of three dimensions we consider a curve ⃗ in N -dimensions
                                                                    ′ , ⃗′′ , ⃗′′′ , · · · , ⃗ (k) } where
and define its properties in terms of the N -dimensional vectors {⃗ c c
                                                                   c                         c
the prime denotes differentiation with respect to the arc length s the arc length traced by the
curve from a chosen initial point. The curve is planar (embedded in N -dimensions) if only
the first two of the vectors are linearly independent and so on. Thus the local dimensionality
of the curve can be obtained by computing the above vectors successively and calculating the
corresponding Gram determinant, We have done this for k = 2 and k = 3 for the linear and
nonlinear networks.

    For completeness we provide some of the equations used in the program below. Let
      ...
c, c, c denote the first, second, and third derivatives of the concentration with respect to
˙ ¨                                                                              √
                                                                                    ˙ ˙
time respectively. The arc length is defined by (ds)2 = d⃗ · d⃗. Therefore, s = ⃗ · ⃗ .
                                                           c c               ˙     c c
                                                                    √
    We repeatedly use the fact f   ′ (s) ≡ df = f˙ where s =  ˙        ˙ ˙
                                                                      ⃗ · ⃗ to obtain the the
                                                                      c c
                                            ds    s˙
derivatives with respect to the arc length in terms of the time derivatives. We have
                                            ˙
                                            c
                                            ⃗
                                 ⃗ ′ (s) =
                                 c
                                            ˙
                                            s
                                            ⃗¨
                                             c     ˙
                                                   c¨
                                                   ⃗s
                                ⃗ ′′ (s) = 2 − 3
                                 c
                                            ˙
                                            s... s ˙
                                            s⃗ − ⃗ s
                                            ˙c     c˙ ...     ¨
                                                              s
                                ⃗ ′′′ (s) =
                                c                         − 3 2 ⃗ ′′ (s)
                                                                c
                                                s4
                                                ˙            s˙
                                                  c
These can be used to calculate the derivatives of ⃗ with respect to the arc length s knowing
the time derivatives. Now the equation of motion yields ⃗:c˙
                               ∑                ∑
                          ˙
                          ci =     Mij cj (t) +      Bijk cj (t)ck (t)
                                      j                 j,k

for system with a quadratic nonlinearity as in the case of NFκB network, for example. For
                                                                                     c
the linear network B vanishes. We can obtain the higher-order time derivatives of ⃗ by




                                                   20
differentiation iteratively. The time derivatives of s are displayed below:
                                         √
                                    ˙
                                   s =       ˙ ˙
                                            ⃗·⃗
                                            c c
                                          ¨ ˙
                                         ⃗·⃗
                                          c c
                                   ¨
                                   s =
                                           ˙
                                           s       ...
                                  ...     ¨·⃗ +⃗ · ⃗
                                          c ¨ ˙
                                         ⃗ c c c          ¨
                                                          s
                                   s =                 − s. ¨
                                               ˙
                                               s          ˙
                                                          s
   Knowing the vectors ⃗′ , ⃗′′ , ⃗′′′ the Gram determinant can be evaluated using the dot
                         c c c
products numerically as a function of time.




                                              21
Studying the role of modifications likely to occur in experiments.

The EGFR model of Kholodenko et al (Ref. 27 in the main text) was used to address the role of
the effects that are likely to occur in experimental datasets where it is not possible to assay all
molecular species and measure protein activations at very small time intervals in modifying the
ability of the principal components to capture variance in the data sets. In particular, we wanted
to investigate why the variations captured (71%) by the top three components for the
experiments in Gaudet et al. (Ref. 8 in the main text) was significantly less compared to the
variation explained by the top three components for the in-silico networks we studied. We
compared our results with the experiments in Gaudet et al. (Ref. 8 in the main text) where the
EFGR signaling network is one of the cross-reacting signaling networks. All the simulations
were performed using 1200 different parameter sets.

(i) PCA was performed on a sub set of species instead of all the species that are generated by this
signaling network. This is more like to occur in experiments where is it difficult to assay every
molecular species that is generated due to the biochemical reactions. Here, we considered the
species EGF, Grb2, Grb2_Sos, Shc, Sos, and R as described in the BioNetGen file
(http://bionetgen.org/index.php/Egfr_path) for the above network that contains 18 species.

(ii) We performed PCA on data with the concentrations chosen in (i) scaled by their
concentrations at t=0. This probes the effect of using fold changes instead of absolute
concentrations of proteins to measure activation as done in the analysis of experiments by
Gaudet et al. (Ref. 8 in the main text).

(iii) Next we probed the effect of choosing fewer data points in time separated by larger time
intervals than that used for the results in Fig. 2D. There were more than 1000 time points in each
kinetic trajectory used for Fig. 2D. We used 13 time points, 0 t S /1440, 5 t S /1440, 15 t S /1440, 30
t S /1440, 60 t S /1440, 90 t S /1440, 120 t S /1440, 240 t S /1440, 480 t S /1440, 720 t S /1440, 960 t S /1440,
1200/ t S 1440, 1440 t S /1440, where the time t S denotes the time where the system reaches the
steady state. The specific time points were selected to better match the time points in the data
given in Gaudet et al. (Ref. 8 in the main text).

(iv) We added random noise in the species concentrations. This was done to investigate the effect
of having random variability in the data arising from experimental measurements in Gaudet et al.
(Ref. 8 in the main text). In our simulations, for each protein concentration obtained from the
solution of the ODEs, a random number is drawn from either [.001,1] or [1,1000] with equal
probability, and multiplied any concentration.
Conclusion: Our results show that introducing the effects described above can in some cases
reduce the percentage explained by the top eigenvectors. On the other hand, we find that, overall
the above modifications increased the percentage explained by the top eigenvector, e.g. the peak
of the minimum percentage explained distributions increased by about 30%. It is possible that
for a specific set of parameters the minimum percentage explained decreased after the
modifications. However, even in that case, the minimum value of the minimum percentage
explained distributions is around 40% for all the modifications. Therefore, the effects considered
above are unlikely to underlie the rather low value observed in the experimental dataset given in




                                                       22
Gaudet et al. (Ref. 8 in the main text) where the top three principal components capture only
71% of the total variance.



Fig. S11: Results showing effects of different modifications in the EGFR network
(Kholodenko et al. )
                                                                                    Minimum Percent Explained by First Component
                                                                                              Kholodenko et. al., 1999
                                                   140
                                                                                                                                                       Experiment 0


                                                   120



                                                   100



                                                    80
                                  Frequency




                                                    60



                                                    40



                                                    20



                                                     0
                                                         40             50                   60                       70                  80                    90                  100
                                                                                     Minimum Percent Explained by Top Component




A. Result of PCA when all the species in the network described by Kholodenko et al were
considered. The time points used are the same as in Fig. 2D. 1200 different parameter sets were
used. This graph represents the ‘control’ in our numerical experiments.
                             Minimum Percent Explained by Top Component
                                       Kholodenko et. al., 1999                                                                    Change in Minimum Percent Explained by Top Component Vs Control
            800                                                                                                                                        Kholodenko et. al., 1999
                                                                              Experiment 1                          120
                                                                                                                                                                                          Experiment 1
            700

                                                                                                                    100
            600


            500                                                                                                      80
Frequency




                                                                                                        Frequency




            400
                                                                                                                     60

            300

                                                                                                                     40
            200


            100                                                                                                      20


              0
                  60   65   70                75           80      85        90        95         100                 0
                                                                                                                       -30   -20           -10         0          10           20       30           40   50
                                                         Percent
                                                                                                                                                           Change in Percent




B. (Left) Results from the modification described in (i). (Right) Distribution of the percent
changes for any parameter set as (i) was carried out. A positive change implies increase in the
percent explained when (i) was carried out and vice versa.




                                                                                                                    23
                                                          Minimum Percent Explained by First Component
                                                                    Kholodenko et. al., 1999
                              400
                                                                                                            Experiment 2

                              350


                              300


                              250




                  Frequency
                              200


                              150


                              100


                               50


                                0
                                    40        50                60             70             80                  90            100
                                                           Minimum Percent Explained by Top Component


C. Results from the simulations when both the effects described in (i) and (ii) were included.




                                                        Minimum Percent Explained by First Component
                                                                  Kholodenko et. al., 1999
                              700
                                                                                                       Experiment 3


                              600



                              500



                              400
                 Frequency




                              300



                              200



                              100



                                0
                                    80   82        84         86        88         90       92         94        96        98
                                                        Minimum Percent Explained by Top Component

D. Results from simulations when both the modifications described in (i) and (iii) are included.




                                                                              24
                                              Minimum Percent Explained by First Component
                                                        Kholodenko et. al., 1999
                         500
                                                                                             Experiment 4
                         450

                         400

                         350

                         300
             Frequency
                         250

                         200

                         150

                         100

                          50

                           0
                               50   55   60      65      70       75        80     85          90       95   100
                                              Minimum Percent Explained by Top Component




E. Results from simulations when the modifications described in (i), (ii), and (iii) are included.

                                              Minimum Percent Explained by First Component
                                                        Kholodenko et. al., 1999
                         400
                                                                                             Experiment 5

                         350


                         300


                         250
             Frequency




                         200


                         150


                         100


                          50


                           0
                               50   55   60      65      70       75        80     85          90       95   100
                                              Minimum Percent Explained by Top Component



F. Results from simulations when the modifications described in (i), (ii), and (iv) are included.




                                                                       25
                                                  Minimum Percent Explained by First Component
                                                            Kholodenko et. al., 1999
                             120
                                                                                                 Experiment 6


                             100



                              80
                 Frequency




                              60



                              40



                              20



                               0
                                   55   60   65        70        75         80       85          90        95   100
                                                  Minimum Percent Explained by Top Component




G. Results from simulations when all the four modifications described in (i), (ii), (iv), and (iv)
are included.




                                                                   26
Table II: Number of species and the time points used in the data matrix for the biological
networks studied in the paper.



                Network                                        Number of species
         Ras Activation Network                                      14
  EGFR Signaling Network: Smaller version                            18
       (Kholodenko et al., Ref. 27)
  EGFR Signaling Network: Larger version                               356
          (Blinov et al., Ref. 26)
    NF-κB Activation Network (Ref. 29)                                 26


The smallest time intervals for calculating the species concentrations and the time the system
takes to reach the steady state concentrations were chosen using the largest and smallest time
scales estimated from the rate constants and initial concentrations in a reaction network.
Therefore, we had a large number (>103) of time points in our simulations. See the section on,
“Details of the simulation method” and Table 1 in the supplementary material for more details.




                                               27
             Fig. S12: Effect of unit variance scaling on percentage explained by top principal
             components in the Ras activation network

                                 Minimum Percent Explained by Top Component                                                                                         Minimum Percent Explained by Top Two Components
                                              Ras-Sos Network                                                                                                                       Ras-Sos Network
            35                                                                                                                           25
                                                                                 Unit-Variance Scaling                                                                                                     Unit-Variance Scaling

            30
                      (A)                                                                                                                20
                                                                                                                                                    (B)
            25



            20                                                                                                                           15
Frequency




                                                                                                                           Frequency
            15
                                                                                                                                         10

            10


                                                                                                                                          5
             5



             0
                 35    40   45     50           55               60    65         70       75         80         85                       0
                                                            Percent                                                                           70               75           80            85             90              95        100
                                                                                                                                                                                        Percent
                                                                                 Minimum Percent Explained by Top Three Components
                                                                                                  Ras-Sos Network
                                                       80
                                                                                                                                       Unit-Variance Scaling

                                                       70
                                                                 (C)
                                                       60


                                                       50
                                           Frequency




                                                       40


                                                       30


                                                       20


                                                       10


                                                        0
                                                            88              90                   92                   94                           96               98
                                                                                                           Percent




             We considered the kinetics of the molecular species involved in the Ras activation
             network (Fig. S6A). Each species concentration at any time point was scaled by the
             variance for that particular species calculated by using concentrations of that species at
             all the time points (from t=0 to the end time of the simulation) as independent variables.
             Then we performed principal component analysis on the scaled kinetic trajectories at all
             the time points as described in the main text. We show results for the distributions of the
             minimum percentage explained by the top three components calculated over 130 trials.
             The minimum percentage explained (shown in (A)) by the top component was peaked at
             a smaller value compared to case when the variables were not scaled (see Fig. 2C in the
             main text). However, as we included the top two (shown in (B)) and three (shown in (C))
             principal components, the principal components in the scaled dataset captured significant
             amount (peak at 98% for three components) of variations in the original data set.




                                                                                                                      28
1.
   Das
J,
et
al.
(2009)
Digital
signaling
and
hysteresis
characterize
ras
activation
in

      lymphoid
cells.
Cell
136(2):337‐351.

2.
   Kholodenko
BN,
Demin
OV,
Moehren
G,
&
Hoek
JB
(1999)
Quantification
of
short

      term
signaling
by
the
epidermal
growth
factor
receptor.
J
Biol
Chem
274(42):30169‐
      30181.

3.
   Blinov
ML,
Faeder
JR,
Goldstein
B,
&
Hlavacek
WS
(2006)
A
network
model
of
early

      events
in
epidermal
growth
factor
receptor
signaling
that
accounts
for

      combinatorial
complexity.
Biosystems
83(2‐3):136‐151.

4.
   Hoffmann
A,
Levchenko
A,
Scott
ML,
&
Baltimore
D
(2002)
The
IkappaB‐NF‐kappaB

      signaling
module:
temporal
control
and
selective
gene
activation.
Science

      298(5596):1241‐1245.

5.
   Blinov
ML,
Faeder
JR,
Goldstein
B,
&
Hlavacek
WS
(2004)
BioNetGen:
software
for

      rule‐based
modeling
of
signal
transduction
based
on
the
interactions
of
molecular

      domains.
Bioinformatics
20(17):3289‐3291.

6.
   Mckay
MD,
Beckman
RJ,
&
Conover
WJ
(1979)
Comparison
of
3
Methods
for

      Selecting
Values
of
Input
Variables
in
the
Analysis
of
Output
from
a
Computer
Code.

      Technometrics
21(2):239‐245.










                                            29

				
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
views:0
posted:2/17/2012
language:
pages:30