Proteins with greater influence on network dynamics\\evolve more slowly but are not more essential
{ %Ryan N. Gutenkunst$^\ast$$^\dagger$}
%\affiliation{$^\ast$Theoretical Biology and Biophysics, Los Alamos National Laboratory; $\dagger$Center for Nonlinear Studies, Los Alamos National Laboratory %} %
A fundamental question for evolutionary biology is why rates of evolution vary dramatically between proteins. Perhaps surprisingly, it is controversial how much a protein's functional importance affects its rate of evolution. In most studies, functional importance has been measured on the coarse scale of protein knock-outs, while evolutionary rate has been measured on the fine scale of amino acid substitutions. Here we introduce dynamical influence, a finer measure of protein functional importance. To measure dynamical influence, we first use detailed biochemical models of particular reaction networks to measure the influence of each reaction rate constant on network dynamics. We then define the dynamical influence of a protein to be the average influence of the rate constants for all reactions it is involved in. Using models of a dozen biochemical systems and sequence data from vertebrates, we show that dynamical influence and evolutionary rate are negatively correlated; proteins with greater dynamical influence evolve more slowly. We also show that proteins with greater dynamical influence are not more likely to be essential. This suggests that there are many cellular reactions whose presence is essential for life, but whose quantitative rate is relatively unimportant to fitness. We also provide evidence that the effect of dynamical influence on evolutionary rate is independent of protein expression level, expression specificity, gene compactness, and reaction degree. Dynamical influence offers a finer view of functional importance, and our results suggest that focusing on essentiality may have previously led to an underestimation of the role function plays in protein evolution. %They also suggest that mathematical models of biochemical networks can successfully capture not only how a system works, but also how it evolves.
Gutenkunst.tex
%\documentclass{mbe}
\documentclass[floatfix,twocolumn,amsmath]{revtex4}
\pdfoutput=1
%\documentclass[twocolumn,letter,10pt]{article}
%\usepackage{fullpage}
%\usepackage{amsmath}
%\usepackage[numbers,square,sort&compress]{natbib}
\usepackage{graphicx}
%\bibpunct{[}{]}{,}{n}{}{;}

%\bibliographystyle{unsrtnat}

\bibliographystyle{plos2009}

%\usepackage{url}
\usepackage{color}
\usepackage{subfigure}
%\usepackage{calc}

% Reviewers: Wagner, Drummond, Jianzhi Zhang, Dennis Vitkup

%\renewcommand{\cite}{\cite}

\begin{document}

% Max 120 characters for Letter
\title{Proteins with greater influence on network dynamics\\evolve more slowly but are not more essential}

%\author{{
%Ryan N. Gutenkunst$^\ast$$^\dagger$}\\
%\affiliation{$^\ast$Theoretical Biology and Biophysics, Los Alamos National Laboratory; $\dagger$Center for Nonlinear Studies, Los Alamos National Laboratory
%}
%}
%\author{Ryan N. Gutenkunst\\
%{\small Theoretical Biology and Biophysics, Los Alamos National Laboratory}\\ {\small Center for Nonlinear Studies, Los Alamos National Laboratory}\\
%\url{ryang@lanl.gov}}

\author{Ryan N. Gutenkunst}
\email{ryang@lanl.gov}
\affiliation{Theoretical Biology and Biophysics, Los Alamos National Laboratory}
\affiliation{Center for Nonlinear Studies, Los Alamos National Laboratory}

\date{\today}

%\authorlist{Gutenkunst}
%% Max 50 characters and spaces
%\shorttitle{Dynamical influence and evolution}

% 350 words for Research Article
% 150 words for Letter
%\mbeabstract{
\begin{abstract}
A fundamental question for evolutionary biology is why rates of evolution vary dramatically between proteins.
Perhaps surprisingly, it is controversial how much a protein's functional importance affects its rate of evolution.
In most studies, functional importance has been measured on the coarse scale of protein knock-outs, while evolutionary rate has been measured on the fine scale of amino acid substitutions.
Here we introduce dynamical influence, a finer measure of protein functional importance.
To measure dynamical influence, we first use detailed biochemical models of particular reaction networks to measure the influence of each reaction rate constant on network dynamics.
We then define the dynamical influence of a protein to be the average influence of the rate constants for all reactions it is involved in.

Using models of a dozen biochemical systems and sequence data from vertebrates, we show that dynamical influence and evolutionary rate are negatively correlated; proteins with greater dynamical influence evolve more slowly.
We also show that proteins with greater dynamical influence are not more likely to be essential.
This suggests that there are many cellular reactions whose presence is essential for life, but whose quantitative rate is relatively unimportant to fitness.
We also provide evidence that the effect of dynamical influence on evolutionary rate is independent of protein expression level, expression specificity, gene compactness, and reaction degree.
Dynamical influence offers a finer view of functional importance, and our results suggest that focusing on essentiality may have previously led to an underestimation of the role function plays in protein evolution.
%They also suggest that mathematical models of biochemical networks can successfully capture not only how a system works, but also how it evolves.
\end{abstract}

\maketitle

% 3 to 6
%\keywords{dynamical influence, evolutionary rate, essentiality, expression, correlation.}

%\email{ryang@lanl.gov}


%\mbestyle

\section*{Introduction}

Every protein evolves at a characteristic rate~\cite{bib:Zuckerkandl1965}, and these rates vary dramatically from protein to protein.
Understanding this variation is a fundamental challenge for evolutionary biology.
A predictive theory for evolutionary rates might also guide drug design toward targets that evolve slowly, lessening the problem of evolved drug resistance in pathogens.

The recent revolution in genomics has enabled many factors to be examined for correlation with protein evolutionary rate; among them are expression level~\cite{bib:Drummond2005}, gene compactness~\cite{bib:Liao2006}, knock-out essentiality~\cite{bib:Wang2009}, and interaction degree~\cite{bib:Fraser2002}.
Perhaps surprisingly, it is controversial how much a protein's functional importance impacts its rate of evolution~\cite{bib:Park2009}.
Theoretical considerations suggest that functional importance should play a substantial role~\cite{bib:Wilson1977}.
Studies in bacteria~\cite{bib:Jordan2002}, yeast~\cite{bib:Wang2009}, and mammals~\cite{bib:Liao2006} have shown that the essentiality of a gene, i.e. the fitness effect of gene deletion, has only a weak (but statistically significant) correlation with its evolutionary rate.

Protein evolutionary rate is most commonly by the ratio dN/dS of the rate of nonsynonymous (amino acid changing, dN) to synonymous (amino acid preserving, dS) genetic mutation.
Protein importance is most commonly assessed by the fitness effect of knocking that protein out.
Individual substitutions may, however, produce only small changes in a protein's structure and, consequently, small changes in its functional efficacy.
It is thus not clear that the fitness effect of a protein knock-out will correlate well with the fitness effect of small changes in protein efficacy, and this may explain the weak correlation seen between protein importance measured by knock-out and evolutionary rate.

Here we take advantage of another recent revolution in biology, the rise of mathematical modeling, to take an alternative and much finer view of a protein's functional importance.
We assess how small changes in the efficacy of a protein's reactions impact the dynamics of the network in which the protein is embedded.
To do so, we use computational models of specific molecular networks. 
Each of these models condenses the results of many biochemical studies, and we use them to probe the fine scale influence of variation in reaction rate constants on the network dynamics~\cite{bib:Gutenkunst2007}.
We construct a measure of a protein's influence on network dynamics from the influence of each constant for each reaction in which it is involved, and we correlate this with evolutionary rate.
The models we consider were all constructed in mice, rats, or humans, but we expect conservation among vertebrates of network structure~\cite{bib:Fokkens2009} and dynamics, allowing us to assess evolutionary rate using more species.

We find that proteins with greater dynamical influence evolve more slowly.
Moreover, we find that a protein's dynamical influence is uncorrelated with its essentiality as measured by knock-out phenotype.
We do find that in our data that evolutionary rate is also negatively correlated with essentiality, but this correlation is typically weaker than with dynamical influence.
We also examine the correlation of dynamical influence with previously established predictors of protein evolutionary rate.
In particular, we find that although dynamical influence and expression level are positively correlated, they appear to exert independent effects on evolutionary rate.

The advent of systems biology has led to much finer knowledge of the role played by specific proteins in their network contexts.
Here we have leveraged this knowledge to measure protein functional importance by the influence of that protein on the dynamics of its network.
We find that the dynamical influence of a protein is negatively correlated with its evolutionary rate and uncorrelated with its knock-out essentiality.
These observations offer a finer perspective on the role of functional importance in the protein evolution.

\section*{Methods and Materials}

\subsection*{Data and statistics}

%\begin{table}
%\centering
%\input{mapping}
%\caption{Model citations.\label{tbl:citations}}
%\end{table}

Dynamical biochemical models encoded in the Systems Biology Markup Language~\cite{bib:Hucka2003} and formulated as systems of ordinary differential equations were obtained from the 14th release of the BioModels database~\cite{bib:LeNovere2006}.
We considered all deterministic models which involved 8 or more distinct proteins, as annotated by their UniProt~\cite{bib:Bairoch2005} identifiers, yielding 10 models.
One additional model~\cite{bib:Borisov2009} was obtained from the Molecular Systems Biology structured data archive, and another~\cite{bib:Albeck2008} was converted to SBML from published Matlab scripts.
%The systems considered and their corresponding publications are listed in table~\ref{tbl:citations}.

Protein sequence alignments were obtained from the Homologene database~\cite{bib:Wheeler2007}, after mapping UniProt identifiers to UniGene identifiers. Corresponding DNA sequences were obtained from GenBank.

\begin{figure}
\centering
\includegraphics[width=2.5in]{alltree}
\caption{Species tree used for measuring evolutionary rate. This tree is that obtained from concatenating all sequences studied, and the scale bar is in synonymous substitutions.\label{fig:tree}}
\end{figure}

PAML~\cite{bib:Yang2007paml} was used to obtain a maximum-likelihood estimate of dN/dS for each protein over the tree shown in Figure~\ref{fig:tree}.
Only proteins for which 3 or more sequences were available were considered.

Following~\cite{bib:Liao2006}, proteins were identified as essential or non-essential based on targeted deletion experiments compiled in the Mouse Genome Database (\url{http://www.informatics.jax.org}).
A protein was deemed essential if it resulted in a phenotype of ``lethality-prenatal/perinatal'' or ``lethality-postnatal'', including more specific classifications.
We expect this essential vs non-essential characterization to be well-conserved among vertebrates.

Mouse mRNA average expression and tissue specificity data were obtained from~\cite{bib:Liao2006}.

Correlations are assessed using Pearson's linear correlation coefficient. Reported p-values are from one-sided permutation tests.
They correspond to the probability by chance of the correlation having the expected sign and being of the observed magnitude or larger.

\subsection*{Dynamical influence}

\begin{figure}
\begin{center}
\subfigure[Example system]{\input{rxn.pdf_t}}
\subfigure[Dynamical trajectories]{\includegraphics[width=1.5in]{trajectories}}\label{fig:trajectories}
\end{center}
\caption{(a) Reactions indicated with bold arrows in this example network are considered when calculating the dynamical influence of protein A. Note also that the external protein ligand L would not be considered in our analysis. (b) Dynamical trajectories are shown for the active forms of proteins A and C. The bold lines are trajectories given the basal rate constants, while the dark and light shaded regions indicate the change in dynamics upon increasing $k_1^+$ or $k_4^-$, respectively. In this case $k_1^+$ has a larger effect than $k_4^-$, so the dynamical influence $\kappa_1^+$ of $k_1^+$ is greater than the influence $\kappa_4^-$ of $k_4^-$.}\label{fig:illus}
\end{figure}

To quantify the change in network dynamics as the reaction rate constants $\mathbf{k}$ vary from their published values $\mathbf{k}^*$, we introduce the sensitivity measure $\kappa_i$:%~\cite{bib:Gutenkunst2007}:
\begin{equation}
\kappa_i^2 %= &\, \frac{d^2 \chi^2}{d^2 \log k_i}\bigg|_{\mathbf{k} = \mathbf{k^*}}\\
 = \sum_{\substack{\text{species }y \\ \text{conditions } c}} \frac{1}{T_c \sigma_y^2} \int_0^{T_c} \bigg(\frac{d y_c(t, \mathbf{k})}{d \log k_i}\bigg)^2 \bigg|_{\mathbf{k} = \mathbf{k^*}} dt\label{eqn:kappa}.
\end{equation}
%\chi^2(\mathbf{k}) = \sum_{\substack{\text{species }y \\ \text{conditions } c}} \frac{1}{2 T_c} \int_0^{T_c} \left(\frac{y_c(t,\mathbf{k}) - y_c(t, \mathbf{k^*})}{\sigma_y}\right)^2 dt.\label{eqn:chi2}
%\end{equation}
This measure sums over all molecular species $y$ in the model, and we sum over all conditions $c$ considered in the original model publication.
Most commonly, the conditions considered are different levels of stimulation by an external ligand.
The normalization $\sigma_y$ is taken to be the maximum of $y$ over all conditions considered. 
See Figure~\ref{fig:illus}(b) for an illustration of how varying rate constants changes dynamics.
%We defined the dynamical influence of rate constant $k_i$ by
%This measure %XXX describe it
The $\kappa_i$ were calculated using SloppyCell (\url{http://sloppycell.sourceforge.net}), and the scripts used are provided in supplementary data.

We define the  dynamical influence of a protein to be the geometric mean of the influence of all rate constants for all reactions it is involved in:
\begin{equation}
D = \left<\kappa_i\right>_\mathrm{geom}\label{eqn:I}.
\end{equation}
%XXX: Arithmetic mean => Results not as convincing.
The $\kappa_i$ are only considered for reactions in which the protein participates by itself, not in a complex.
See Figure~\ref{fig:illus}(a) for an example.
Several models consider perturbations by externally manipulated protein ligands; these ligands
% (EGF, NGF, heregulin, insulin, and IL-6)
were excluded from consideration, as the manipulations considered in the models may not be physiological.
Additionally, the TNF receptor was excluded from our analysis of the apoptosis model~\cite{bib:Albeck2008}, because it was modeled non-mechanistically.

\section*{Results}

\subsection*{Proteins with greater dynamical influence evolve more slowly}

To measure the influence of a particular protein on network dynamics, we first consider the influence of each rate constant of each reaction on network dynamics.
The influence $\kappa_i$ of rate constant $i$ is defined as the summed and integrated change in \emph{all} molecular species trajectories as that rate constant is changed (Equation~\ref{eqn:kappa}).
See Figure~\ref{fig:illus}(b) for an illustration of how changes in rate constants are reflected in changes in dynamics.
We sum influence over all molecular species because it is not clear which aspects of each network's dynamics are relevant to an animal's fitness.
Results are, however, qualitatively similar if we consider changes in only a select subset of ``output'' molecular species (see online supporting data).
Given the influence of each rate constant, the influence $D$ of a protein is then defined as the geometrical mean of the influence of all rate constants for all reactions it is involved in (Figure~\ref{fig:illus}(a)).

We measure the evolutionary rate of each protein via the ratio of the rate of nonsynonymous (amino acid changing, dN) to synonymous (amino acid preserving, dS) genetic mutation, as inferred over the tree seen in Figure~\ref{fig:tree}.
(Results are qualitatively similar if we consider only the tree of mammals, rather than vertebrates (see supporting data).)

\begin{figure}
\centering
\includegraphics[width=\columnwidth]{correlation}
\caption{Relationship between protein evolutionary rate and dynamical influence. Each point represents a particular protein in each system. Proteins deemed essential in mouse knock-out studies are plotted in black, non-essential in grey. There was no essentiality data available for proteins plotted in white.\label{fig:correlation}}
\end{figure}

Figure~\ref{fig:correlation} plots, for each system, protein evolutionary rate dN/dS versus protein influence $D$ on network dynamics.
We consider each system separately, because it is not clear that changes in the dynamics of, for example, arachidonic acid metabolism and ErbB signaling would yield comparable fitness effects.
Table~\ref{tbl:main} summarizes these results.

\begin{table*}
\centering
\input{table1}
\caption{Correlations between evolutionary rate dN/dS, dynamical influence $D$, and essentiality $E$. p-values are from one-sided permutation tests, and N is the number of proteins in each correlation.\label{tbl:main}}
\end{table*}

We expect that proteins with greater influence on network dynamics are functionally more important and thus should evolve more slowly.
For 11 of the 12 systems studied, there is indeed a negative correlation between protein evolutionary rate and dynamical influence (Table~\ref{tbl:main}).
The single exception is the $\beta$-adrenergic signaling model, which may be unreliable because it is a well-mixed approximation to a spatial model that was actually fit to data.
(Technical limitations presently preclude us from considering the partial-differential-equation spatial model in this analysis.)
Individual p-values are not dramatic; each plot includes few points because building such detailed dynamic models is currently low-throughput, and the models are thus relatively small.
If all the systems were completely independent, the odds of 11 of 12 showing a negative correlation by chance are $12/2^{12}$, which is less than 0.3\%.
The first 7 systems in Table~\ref{tbl:main} do share some proteins between them, because the EGF signaling pathway is the most often-modeled system in higher eukaryotes.
If we eliminate all proteins shared between two or more systems, then 10/11 show a negative correlation, expected by chance at the 0.5\% level.
(The smaller NGF/EGF signaling model is eliminated in this comparison because all save one of the proteins in it are shared.)
See supporting data for details.

The ratio dN/dS is affected by both changes in the substitution rate dN and in the neutral mutation rate dS.
In the supporting data, we show that the correlation between dN/dS and $D$ we observe here is not due to correlation with dS.
%Table~\ref{tbl:main} shows that this correlation is not due to variation in the mutation rate, as indicated by dS between human and mouse.

\subsection*{Proteins with greater dynamical influence are not more essential}

The most common measure of a protein's functional importance is essentiality, whether or not an organism can survive without it.
As in~\cite{bib:Liao2006}, we assessed essentiality by whether targeted knock-out of a protein yielded lethality in mice.
Table~\ref{tbl:main} shows that 10 of 12 systems show the expected negative correlation between essentiality $E$ and evolutionary rate, and that the correlations are typically weaker than that between dynamical influence and evolutionary rate.

Perhaps surprisingly, we find little correlation between dynamical influence and essentiality.
The correlation is only positive for only 7 of 12 models (p-value 0.39).
This suggests that the correlation we observe between $D$ and dN/dS is independent of the correlation between $E$ and dN/dS.
As seen in Table~\ref{tbl:main}, the partial correlations between evolutionary rate and dynamical influence typically remain negative when controlling for essentiality.
Of the 11 systems with a negative correlation between dN/dS and $D$, three do change sign when controlling for essentiality.
It is difficult to assess the significance of these sign flips, however, because in two of these three cases only one protein was essential or nonessential.

\subsection*{Other potential correlates}

\begin{table*}
\centering
\footnotesize
\input{expression}
\caption{Correlations between evolutionary rate dN/dS, dynamical influence D, and expression-related variables. p-values are from one-sided permutation tests, and N is the number of proteins in each correlation\label{tbl:other}}
\end{table*}

The most well-established correlate with protein evolutionary rate is level of expression. Table~\ref{tbl:other} reports correlations in our data between evolutionary rate and expression-related variables.

To measure expression level, we used the average expression level across 61 mouse tissues, as compiled by~\cite{bib:Liao2006}.
As found in other studies~\cite{bib:Drummond2005}, evolutionary rate is consistently negatively correlated with expression level in our data (Table~\ref{tbl:other}).
Moreover, dynamical influence and expression level are positively correlated in 8 of 12 cases.
Of the 11 systems for which evolutionary rate and dynamical influence are negatively correlated, in only 2 cases does controlling for expression level via partial correlation change the sign.
This suggests that dynamical influence and expression level, although correlated with each other, exert mostly independent effects on evolutionary rate.

In vertebrates, expression level may vary from tissue to tissue, so we also we consider data on tissue specificity $\tau$ compiled by~\cite{bib:Liao2006}.
Tissue specificity is consistently positively correlated with evolutionary rate (Table~\ref{tbl:other}); proteins expressed in few tissues evolve faster.
Dynamical influence and specificity are negatively correlated in 9 of 12 cases, indicating that influential proteins tend to be expressed in many tissues.
Controlling for tissue specificity, dynamical influence and evolutionary rate remain negatively correlated for 10 of the 11 systems for which they began negatively correlated, suggesting that these factors independently influence the rate of protein evolution.

We have also examined correlations with gene compactness~\cite{bib:Liao2006} and the number of reactions a protein participates in (see supporting data).
We find no consistent correlation between a protein's dynamical influence and it's genetic compactness.
We do find that there is a negative correlation between a protein's dynamical influence and the number of reactions it participates in.
In both cases, however, it appears that dynamical influence exerts an independent effect on evolutionary rate.

\section*{Discussion}

We have introduced a fine scale measure of protein functional importance, based on the influence each protein has on the dynamics of the biochemical network in which it is embedded.
We have shown that in the systems we study proteins with greater dynamical influence evolve more slowly, although they are no more essential as determined by mouse knock-outs.
Moreover, the correlation of dynamical influence with evolutionary rate appears to be independent of expression level and specificity.
It also appears to be independent of gene compactness and number of reactions a protein is involved in.

Of the 12 systems tested, 11 show a negative correlation between protein dynamical influence and evolutionary rate.
The one exception, a model of $\beta$-adrenergic signaling, is a well-mixed approximation to a spatial model for which we do not currently have the techniques to analyze parameter sensitivity, and it was the spatial model that was actually fit to biochemical data.
The approximate model may thus not properly reflect the actual network dynamics, explaining why this system is an exception.
This emphasizes the importance of model quality in assessing dynamical influence.
We measured dynamical influence using computational models because it is at present technically very challenging to experimentally introduce small controlled perturbations to protein activity.
All models, however, are approximations to reality, which introduces noise into our measure of dynamical influence.
In particular, the models currently available of necessity omit many potential network components and interactions.
These omissions may explain the many cases in which a protein has low measured dynamical influence but evolves slowly. 
(Points clustered near 0,0 in Figure~\ref{fig:correlation}).
The rate constant values used in the current models are also likely inaccurate.
It appears, however, that these models nevertheless capture the correct sensitivity to changes in rate constants.
As more complete and refined models are developed, we may find that the observed correlation between dynamical influence and evolutionary rate grows stronger.

Many proteins are modular, with different domains performing different reactions.
Here we measure sequence divergence of whole proteins and define dynamical influence using all reactions in which a protein is involved.
Considering the evolution of each protein domain separately and computing its influence using only the reactions it is actively involved in may make our analysis more powerful.
Similarly, we eliminated protein complexes from analysis, but future studies may identify which components of a complex participate in which reactions.
Both these refinements are difficult at present, because models are not typically annotated with the necessary information.
As domain-oriented rules-based models~\cite{bib:Hlavacek2006} grow in prominence, this limitation may disappear.

We measure protein dynamical influence by first assessing how changing each reaction rate constant affects network dynamics, and we measure evolutionary rate by sequence divergence.
However, changes in a protein sequence do not necessarily correspond to changes in rate constants~\cite{bib:Kiel2008}.
Ideally, we would directly measure the evolutionary divergence in rate constants, but at present these quantities are difficult to measure in even a single organism.
Computational methods that predict the impact of particular sequence changes on protein structure (e.g.~\cite{bib:Sunyaev2001}) may provide a more relevant measure of evolutionary change.
However, given that that we see a correlation between dynamical influence and protein sequence evolution, we expect that the correlation between dynamical influence and more direct measures of rate constant evolution would be even stronger.

A further underlying assumption in our analysis is that the structure of the network remains constant over the evolutionary distances considered.
A recent study of network structure evolution in eukaryotes indeed found that it is typically well-conserved among vertebrates.
The stability of the networks we consider is further supported by the fact that the correlations we observe are stronger when the tree of vertebrates is used rather than the tree of mammals.
The reduction in noise resulting from using more sequence in assessing dN/dS thus appears to outweigh noise from possible changes in network structure or function.

Here we have examined evolutionary rate using biochemical pathway models.
Other authors have also taken a pathway perspective evolutionary rate.
Ramsay et al.\ showed that proteins acting earlier in particular plant metabolic synthesis pathways evolved more slowly than those downstream~\cite{bib:Ramsay2009}.
Similarly, Greenberg et al.\ showed that proteins in the metabolic ``core'' or that share metabolites with other enzymes evolve more slowly~\cite{bib:Greenberg2008}.
Other authors have used flux-balance analysis to consider the steady-state properties of the yeast metabolic network:
Vitkup et al.\ showed that proteins mediating greater metabolic flux evolved more slowly~\cite{bib:Vitkup2006}, while Wang et al. showed that proteins with larger knock-out effect on the steady-steady metabolic network evolve more slowly~\cite{bib:Wang2009}.
Here we consider not only structure or steady-state behavior, but also dynamics.
We also consider a different class of networks, as most of our examples are signaling networks, the exceptions being models of the cell cycle and arachidonic acid metabolism.

It is perhaps surprising that in these systems dynamical influence is independent of knock-out essentiality.
This may reflect the presence of many reactions which must take place for survival, but whose precise rate matters little, perhaps due to other forms of regulation.
Proteins involved in such reactions would have dramatic effects if knocked out (and thus be deemed essential), but the network would be relatively insensitive to changes in their sequence (and thus they would have low dynamical influence).
%Also note that dynamical influence may differentiate between proteins which are essential.
%A particular example from this study is the~\citet{bib:Borisov2009} system (Figure~\ref{fig:correlation}).

A great deal of effort has been devoted to teasing apart which correlates of evolutionary rate exert independent effects.
This has proven very difficult, and there is no well-established statistical test that may be used.
Drummond et al.\ have shown that partial correlation may overestimate independence between noisy factors, and they argue in favor of principle components regression~\cite{bib:Drummond2006}.
Principle components regression, however, has similar difficulties when variables differ in their noisiness~\cite{bib:Plotkin2007}.
Here we have used partial correlation analysis, as it is simple and can be applied to the small number of data points we have for each system.
Our results are thus suggestive of independent effects between dynamical influence and other factors, but they are not definitive.
It does seem, however, that essentiality and dynamical influence are themselves almost uncorrelated, implying that their effects on evolutionary rate are independent.

We have introduced a measure of protein functional importance, dynamical influence, that assesses the effects of small changes in protein reaction efficacy on network dynamics.
Dynamical influence is negatively correlated with evolutionary rate as measured by dN/dS, but it is independent of essentiality, a commonly used but much coarser measure of importance.
By taking a finer view of protein importance, our results shed light on the controversy as to how much functional importance matters for evolution, suggesting that more important proteins do indeed evolve more slowly.

\begin{acknowledgments}
I thank Sarah Stockwell, Tony Greenberg, Tanja Kortemme and particularly Simon Gravel for helpful comments and suggestions. I also thank Jianzhi Zhang and Ben-Yang Liao for graciously sharing their compiled expression data.
\end{acknowledgments}

%which sums over all species in the network and integrates over time how the dynamics how the dynamics of the network have changed when parameters change to a value $\mathbf{k}$ from their basal (published) values $\mathbf{k}^*$~\cite{bib:Gutenkunst2007}.

%To measure change in network dynamics, we take an agnostic approach that sums over all dynamical species $y$ in the model, to define a change function $\chi^2$ which measures how the dynamics of the network have changed when parameters change to a value $\mathbf{k}$ from their basal (published) values $\mathbf{k}^*$~\cite{bib:Gutenkunst2007}.

%We view the dynamics of the network through the lens of the biochemical rate constants that quantify the reations, for example bindnig constants between proteins, or enzymatic rates. Mutations that change protein amino acid sequence will change these reaction rate constants. Thus we begin by examining the sensitivity of the network dynamics to variations in rate constants.

%\bibliographystyle{MBErefs}
\bibliography{../ProteinEvol}

\end{document}

Nrxn.tex
\begin{tabular}{cccc}
system & $r_{dN/dS,Nrxn}$ (p-val, N) & $r_{D,Nrxn}$ (p-val, N) & $r_{dN/dS,D|Nrxn}$\\
\hline
Insulin/EGF signaling \cite{bib:Borisov2009} & +0.12 (0.65, 13) & +0.32 (0.86, 13) & -0.42\\
Heregulin/EGF signaling \cite{bib:Birtwistle2007} & +0.03 (0.55, 14) & -0.60 (0.01, 14) & -0.31\\
NGF/EGF signaling (smaller) \cite{bib:Brown2004} & -0.16 (0.35, 8) & -0.40 (0.17, 8) & -0.79\\
NGF/EGF signaling (larger) \cite{bib:Sasagawa2005} & +0.05 (0.57, 17) & +0.23 (0.82, 17) & -0.45\\
EGF receptor endocytosis \cite{bib:Ung2008} & -0.40 (0.02, 25) & -0.27 (0.09, 25) & -0.24\\
Wnt/ERK crosstalk \cite{bib:Kim2007} & -0.43 (0.11, 10) & -0.47 (0.09, 10) & -0.66\\
IL6 signaling \cite{bib:Singh2006} & -0.42 (0.03, 10) & +0.37 (0.82, 10) & -0.14\\
G1 cell cycle progression \cite{bib:Haberichter2007} & -0.32 (0.20, 10) & -0.32 (0.18, 10) & -0.45\\
Rho-kinase activation \cite{bib:Maeda2006} & -0.02 (0.52, 16) & -0.26 (0.10, 16) & -0.22\\
Arachidonic acid metabolism \cite{bib:Yang2007} & +0.03 (0.53, 10) & -0.26 (0.25, 10) & -0.54\\
TNF-driven apoptosis \cite{bib:Albeck2008} & --- & --- & ---\\
$\beta$-adrenergic signaling \cite{bib:Neves2008} & +0.15 (0.68, 10) & -0.54 (0.05, 10) & +0.57\\
\hline
expected sign & $-$ & $+$ & $-$ \\
 total of expected sign & 6 & 3 & 10 \\
\end{tabular}
compactness.tex
\begin{tabular}{cccc}
system & $r_{dN/dS,C}$ (p-val, N) & $r_{D,C}$ (p-val, N) & $r_{dN/dS,D|C}$\\
\hline
Insulin/EGF signaling \cite{bib:Borisov2009} & +0.11 (0.62, 11) & +0.50 (0.94, 11) & -0.49\\
Heregulin/EGF signaling \cite{bib:Birtwistle2007} & +0.20 (0.75, 14) & +0.01 (0.53, 14) & -0.27\\
NGF/EGF signaling (smaller) \cite{bib:Brown2004} & +0.29 (0.74, 7) & -0.39 (0.20, 7) & -0.59\\
NGF/EGF signaling (larger) \cite{bib:Sasagawa2005} & +0.00 (0.52, 16) & -0.11 (0.38, 16) & -0.41\\
EGF receptor endocytosis \cite{bib:Ung2008} & -0.01 (0.50, 25) & +0.39 (0.96, 25) & -0.11\\
Wnt/ERK crosstalk \cite{bib:Kim2007} & +0.73 (0.99, 10) & -0.31 (0.19, 10) & -0.15\\
IL6 signaling \cite{bib:Singh2006} & +0.07 (0.69, 10) & +0.63 (0.96, 10) & -0.41\\
G1 cell cycle progression \cite{bib:Haberichter2007} & +0.41 (0.88, 10) & +0.25 (0.76, 10) & -0.46\\
Rho-kinase activation \cite{bib:Maeda2006} & +0.74 (1.00, 16) & -0.34 (0.06, 16) & +0.07\\
Arachidonic acid metabolism \cite{bib:Yang2007} & -0.55 (0.06, 10) & +0.28 (0.81, 10) & -0.47\\
TNF-driven apoptosis \cite{bib:Albeck2008} & -0.44 (0.07, 12) & +0.04 (0.58, 12) & -0.14\\
$\beta$-adrenergic signaling \cite{bib:Neves2008} & +0.22 (0.75, 10) & -0.41 (0.14, 10) & +0.54\\
\hline
expected sign & $+$ & $-$ & $-$ \\
 total of expected sign & 9 & 5 & 10 \\
\end{tabular}
expression.tex
\begin{tabular}{ccccccc}
& \multicolumn{3}{c}{average expression $X$} & \multicolumn{3}{c}{expression specificity $\tau$}\\
system & $r_{dN/dS,X}$ (p-val, N)& $r_{D,X}$ (p-val, N)& $r_{dN/dS,D|X}$ & $r_{dN/dS,\tau}$ (p-val, N) & $r_{D,\tau}$ (p-val, N) & $r_{dN/dS,D|\tau}$\\
\hline
Insulin/EGF signaling \cite{bib:Borisov2009} & -0.78 (0.01, 11) & +0.77 (1.00, 11) & +0.42 & +0.21 (0.71, 11) & -0.08 (0.40, 11) & -0.43\\
Heregulin/EGF signaling \cite{bib:Birtwistle2007} & -0.45 (0.08, 10) & +0.38 (0.85, 10) & -0.07 & +0.14 (0.66, 10) & -0.05 (0.47, 10) & -0.23\\
NGF/EGF signaling (smaller) \cite{bib:Brown2004} & -0.39 (0.22, 5) & +0.69 (0.88, 5) & -0.88 & +0.84 (0.98, 5) & -0.66 (0.11, 5) & -0.73\\
NGF/EGF signaling (larger) \cite{bib:Sasagawa2005} & -0.53 (0.03, 12) & +0.50 (0.94, 12) & -0.46 & +0.44 (0.92, 12) & -0.28 (0.19, 12) & -0.56\\
EGF receptor endocytosis \cite{bib:Ung2008} & -0.43 (0.01, 23) & +0.19 (0.82, 23) & -0.14 & +0.43 (0.98, 23) & -0.00 (0.50, 23) & -0.23\\
Wnt/ERK crosstalk \cite{bib:Kim2007} & -0.87 (0.00, 7) & -0.38 (0.19, 7) & +0.00 & +0.17 (0.67, 7) & +0.83 (0.98, 7) & +0.34\\
IL6 signaling \cite{bib:Singh2006} & -0.58 (0.02, 9) & -0.12 (0.39, 9) & -0.47 & +0.16 (0.72, 9) & -0.10 (0.45, 9) & -0.30\\
G1 cell cycle progression \cite{bib:Haberichter2007} & -0.07 (0.43, 10) & +0.46 (0.90, 10) & -0.31 & +0.55 (0.95, 10) & -0.43 (0.11, 10) & -0.10\\
Rho-kinase activation \cite{bib:Maeda2006} & -0.09 (0.39, 13) & +0.09 (0.66, 13) & -0.15 & +0.45 (0.93, 13) & +0.01 (0.56, 13) & -0.18\\
Arachidonic acid metabolism \cite{bib:Yang2007} & -0.54 (0.12, 7) & +0.47 (0.77, 7) & -0.23 & +0.00 (0.51, 7) & -0.06 (0.58, 7) & -0.43\\
TNF-driven apoptosis \cite{bib:Albeck2008} & -0.50 (0.08, 9) & -0.04 (0.43, 9) & -0.39 & +0.49 (0.90, 9) & -0.28 (0.22, 9) & -0.22\\
$\beta$-adrenergic signaling \cite{bib:Neves2008} & -0.66 (0.03, 8) & -0.44 (0.13, 8) & +0.42 & +0.09 (0.59, 8) & +0.17 (0.66, 8) & +0.57\\
\hline
expected sign & $-$ & $+$ & $-$ & $+$ & $-$ & $-$ \\
 total of expected sign & 12 & 8 & 9 & 12 & 9 & 10 \\
\end{tabular}
supp_dN_dS.tex
\begin{tabular}{ccc}
system & $r_{dN,D}$ (p-val, N) & $r_{dS,D}$ (p-val, N)\\
\hline
Insulin/EGF signaling \cite{bib:Borisov2009} & -0.29 (0.17, 13) & -0.40 (0.09, 13)\\
Heregulin/EGF signaling \cite{bib:Birtwistle2007} & -0.25 (0.20, 14) & -0.30 (0.15, 14)\\
NGF/EGF signaling (smaller) \cite{bib:Brown2004} & -0.65 (0.03, 8) & -0.46 (0.11, 8)\\
NGF/EGF signaling (larger) \cite{bib:Sasagawa2005} & -0.41 (0.03, 17) & -0.01 (0.46, 17)\\
EGF receptor endocytosis \cite{bib:Ung2008} & -0.07 (0.39, 25) & -0.16 (0.22, 25)\\
Wnt/ERK crosstalk \cite{bib:Kim2007} & -0.34 (0.16, 10) & +0.29 (0.79, 10)\\
IL6 signaling \cite{bib:Singh2006} & -0.04 (0.50, 9) & +0.04 (0.55, 9)\\
G1 cell cycle progression \cite{bib:Haberichter2007} & -0.41 (0.13, 9) & -0.22 (0.30, 9)\\
Rho-kinase activation \cite{bib:Maeda2006} & -0.32 (0.09, 15) & +0.36 (0.93, 15)\\
Arachidonic acid metabolism \cite{bib:Yang2007} & -0.17 (0.32, 10) & +0.56 (0.92, 10)\\
TNF-driven apoptosis \cite{bib:Albeck2008} & +0.21 (0.74, 11) & -0.33 (0.15, 11)\\
$\beta$-adrenergic signaling \cite{bib:Neves2008} & +0.35 (0.80, 8) & +0.07 (0.56, 8)\\
\hline
expected sign & $-$ & $-$ \\
 total of expected sign & 10 & 7 \\
\end{tabular}
supp_important.tex
\begin{tabular}{ccc}
system & $r_{dN/dS,D}$ (p-val, N) & ``output'' variables\\
\hline
Insulin/EGF signaling \cite{bib:Borisov2009} & -0.24 (0.23, 12) & active ERK\\
Heregulin/EGF signaling \cite{bib:Birtwistle2007} & -0.33 (0.16, 12) & active ERK and active Akt\\
NGF/EGF signaling (smaller) \cite{bib:Brown2004} & -0.73 (0.03, 8) & active ERK\\
NGF/EGF signaling (larger) \cite{bib:Sasagawa2005} & -0.40 (0.05, 16) & active ERK\\
EGF receptor endocytosis \cite{bib:Ung2008} & -0.19 (0.21, 24) & active ERK\\
Wnt/ERK crosstalk \cite{bib:Kim2007} & -0.02 (0.54, 10) & $\beta$-catenin/TCF complex\\
IL6 signaling \cite{bib:Singh2006} & -0.05 (0.69, 9) & nuclear STAT3\\
G1 cell cycle progression \cite{bib:Haberichter2007} & -0.57 (0.05, 9) & hypo- and hyper-phosphorylated Rb\\
Rho-kinase activation \cite{bib:Maeda2006} & -0.19 (0.28, 15) & phosphorylated myosin light chain\\
Arachidonic acid metabolism \cite{bib:Yang2007} & -0.57 (0.05, 10) & PGE2 and LTB4\\
TNF-driven apoptosis \cite{bib:Albeck2008} & -0.22 (0.34, 12) & the apoptosome\\
$\beta$-adrenergic signaling \cite{bib:Neves2008} & -0.18 (0.32, 9) & active MAP kinase\\
\hline
expected sign & $-$ \\
 total of expected sign & 12 \\
\end{tabular}
supp_mammals.tex
\begin{tabular}{cc}
system & $r_{dN/dS,D}$ (p-val, N)\\
\hline
Insulin/EGF signaling \cite{bib:Borisov2009} & -0.28 (0.17, 13)\\
Heregulin/EGF signaling \cite{bib:Birtwistle2007} & -0.23 (0.21, 14)\\
NGF/EGF signaling (smaller) \cite{bib:Brown2004} & -0.58 (0.07, 8)\\
NGF/EGF signaling (larger) \cite{bib:Sasagawa2005} & -0.27 (0.15, 17)\\
EGF receptor endocytosis \cite{bib:Ung2008} & -0.03 (0.45, 25)\\
Wnt/ERK crosstalk \cite{bib:Kim2007} & -0.37 (0.14, 10)\\
IL6 signaling \cite{bib:Singh2006} & -0.30 (0.21, 10)\\
G1 cell cycle progression \cite{bib:Haberichter2007} & -0.31 (0.19, 10)\\
Rho-kinase activation \cite{bib:Maeda2006} & -0.10 (0.40, 16)\\
Arachidonic acid metabolism \cite{bib:Yang2007} & -0.49 (0.08, 10)\\
TNF-driven apoptosis \cite{bib:Albeck2008} & -0.11 (0.41, 12)\\
$\beta$-adrenergic signaling \cite{bib:Neves2008} & +0.70 (0.98, 9)\\
\hline
expected sign & $-$ \\
 total of expected sign & 11 \\
\end{tabular}
supp_unique.tex
\begin{tabular}{cc}
system & $r_{dN/dS,D}$ (p-val, N)\\
\hline
Insulin/EGF signaling \cite{bib:Borisov2009} & -0.61 (0.08, 6)\\
Heregulin/EGF signaling \cite{bib:Birtwistle2007} & -0.63 (0.06, 6)\\
NGF/EGF signaling (smaller) \cite{bib:Brown2004} & --- (---, 1)\\
NGF/EGF signaling (larger) \cite{bib:Sasagawa2005} & -0.39 (0.21, 7)\\
EGF receptor endocytosis \cite{bib:Ung2008} & -0.01 (0.49, 10)\\
Wnt/ERK crosstalk \cite{bib:Kim2007} & -0.42 (0.16, 7)\\
IL6 signaling \cite{bib:Singh2006} & -0.28 (0.24, 10)\\
G1 cell cycle progression \cite{bib:Haberichter2007} & -0.31 (0.19, 10)\\
Rho-kinase activation \cite{bib:Maeda2006} & -0.26 (0.20, 14)\\
Arachidonic acid metabolism \cite{bib:Yang2007} & -0.53 (0.06, 10)\\
TNF-driven apoptosis \cite{bib:Albeck2008} & -0.14 (0.36, 12)\\
$\beta$-adrenergic signaling \cite{bib:Neves2008} & +0.26 (0.71, 7)\\
\hline
expected sign & $-$ \\
 total of expected sign & 10 \\
\end{tabular}
table1.tex
\begin{tabular}{ccccc}
system & $r_{dN/dS,D}$ (p-val, N) & $r_{dN/dS,E}$ (p-val, N) & $r_{D,E}$ (p-val, N) & $r_{dN/dS,D|E}$ \\
\hline
Insulin/EGF signaling \cite{bib:Borisov2009} & -0.35 (0.12, 13) & -0.11 (0.37, 12) & -0.14 (0.36, 12) & -0.37\\
Heregulin/EGF signaling \cite{bib:Birtwistle2007} & -0.27 (0.18, 14) & +0.20 (0.74, 13) & -0.43 (0.08, 13) & -0.24\\
NGF/EGF signaling (smaller) \cite{bib:Brown2004} & -0.65 (0.05, 8) & -0.22 (0.43, 7) & +0.05 (0.43, 7) & -0.47\\
NGF/EGF signaling (larger) \cite{bib:Sasagawa2005} & -0.43 (0.04, 17) & -0.16 (0.30, 15) & +0.17 (0.70, 15) & -0.64\\
EGF receptor endocytosis \cite{bib:Ung2008} & -0.10 (0.32, 25) & -0.44 (0.04, 17) & -0.36 (0.08, 17) & -0.44\\
Wnt/ERK crosstalk \cite{bib:Kim2007} & -0.32 (0.18, 10) & -0.98 (0.14, 7) & +0.46 (1.00, 7) & +0.65\\
IL6 signaling \cite{bib:Singh2006} & -0.28 (0.24, 10) & -0.08 (0.34, 9) & +0.18 (0.56, 9) & -0.08\\
G1 cell cycle progression \cite{bib:Haberichter2007} & -0.31 (0.20, 10) & -0.26 (0.31, 7) & -0.28 (0.25, 7) & -0.56\\
Rho-kinase activation \cite{bib:Maeda2006} & -0.21 (0.24, 16) & -0.04 (0.45, 12) & +0.26 (0.75, 12) & -0.23\\
Arachidonic acid metabolism \cite{bib:Yang2007} & -0.53 (0.06, 10) & -0.36 (0.29, 7) & +0.98 (1.00, 7) & +0.67\\
TNF-driven apoptosis \cite{bib:Albeck2008} & -0.14 (0.36, 12) & +0.09 (0.60, 9) & +0.11 (0.61, 9) & +0.08\\
$\beta$-adrenergic signaling \cite{bib:Neves2008} & +0.39 (0.87, 10) & -0.40 (0.24, 7) & -0.21 (0.34, 7) & +0.31\\
\hline
expected sign & $-$ & $-$ & $+$ & $-$ \\
 total of expected sign & 11 & 10 & 7 & 8 \\
\end{tabular}
Gutenkunst.bbl
\begin{thebibliography}{10}
\providecommand{\url}[1]{\texttt{#1}}
\providecommand{\urlprefix}{URL }
\expandafter\ifx\csname urlstyle\endcsname\relax
  \providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}\else
  \providecommand{\doi}{doi:\discretionary{}{}{}\begingroup
  \urlstyle{rm}\Url}\fi
\providecommand{\bibAnnoteFile}[1]{%
  \IfFileExists{#1}{\begin{quotation}\noindent\textsc{Key:} #1\\
  \textsc{Annotation:}\ \input{#1}\end{quotation}}{}}
\providecommand{\bibAnnote}[2]{%
  \begin{quotation}\noindent\textsc{Key:} #1\\
  \textsc{Annotation:}\ #2\end{quotation}}
\providecommand{\eprint}[2][]{\url{#2}}

\bibitem{bib:Zuckerkandl1965}
Zuckerkandl E, Pauling L (1965) Evolutionary divergencs and convergence in
  proteins.
\newblock In: Bryson V, Vogel H, editors, Evolving Genes and Proteins, New
  York: Academic Press. pp. 97--166.
\bibAnnoteFile{bib:Zuckerkandl1965}

\bibitem{bib:Drummond2005}
Drummond DA, Bloom JD, Adami C, Wilke CO, Arnold FH (2005) Why highly expressed
  proteins evolve slowly.
\newblock P Natl Acad Sci USA 102: 14338-14343.
\bibAnnoteFile{bib:Drummond2005}

\bibitem{bib:Liao2006}
Liao BY, Scott NM, Zhang J (2006) {Impacts of gene essentiality, expression
  pattern, and gene compactness on the evolutionary rate of mammalian
  proteins}.
\newblock Mol Biol Evol 23: 2072--2080.
\bibAnnoteFile{bib:Liao2006}

\bibitem{bib:Wang2009}
Wang Z, Zhang J (2009) Why is the correlation between gene importance and gene
  evolutionary rate so weak?
\newblock PLoS Genet 5: e1000329.
\bibAnnoteFile{bib:Wang2009}

\bibitem{bib:Fraser2002}
Fraser HB, Hirsh AE, Steinmetz LM, Scharfe C, Feldman MW (2002) Evolutionary
  rate in the protein interaction network.
\newblock Science 296: 750--752.
\bibAnnoteFile{bib:Fraser2002}

\bibitem{bib:Park2009}
Park D, Choi SS (2009) Why proteins evolve at different rates: The functional
  hypothesis versus the mistranslation-induced protein misfolding hypothesis.
\newblock FEBS Lett 583: 1053--1059.
\bibAnnoteFile{bib:Park2009}

\bibitem{bib:Wilson1977}
Wilson AC, Carlson SS, White TJ (1977) Biochemical evolution.
\newblock Annu Rev Biochem 46: 573-639.
\bibAnnoteFile{bib:Wilson1977}

\bibitem{bib:Jordan2002}
Jordan IK, Rogozin IB, Wolf YI, Koonin EV (2002) Essential genes are more
  evolutionarily conserved than are nonessential genes in bacteria.
\newblock Genome Res 12: 962--968.
\bibAnnoteFile{bib:Jordan2002}

\bibitem{bib:Gutenkunst2007}
Gutenkunst RN, Waterfall JJ, Casey FP, Brown KS, Myers CR, et~al. (2007)
  {Universally sloppy parameter sensitivities in systems biology models}.
\newblock PLoS Comput Biol 3: 1871--1878.
\bibAnnoteFile{bib:Gutenkunst2007}

\bibitem{bib:Fokkens2009}
Fokkens L, Snel B (2009) Cohesive versus flexible evolution of functional
  modules in eukaryotes.
\newblock PLoS Comput Biol 5: e1000276.
\bibAnnoteFile{bib:Fokkens2009}

\bibitem{bib:Hucka2003}
Hucka M, Finney A, Sauro HM, Bolouri H, Doyle JC, et~al. (2003) The systems
  biology markup language ({SBML}): a medium for representation and exchange of
  biochemical network models.
\newblock Bioinformatics 19: 524-531.
\bibAnnoteFile{bib:Hucka2003}

\bibitem{bib:LeNovere2006}
Le~Novere N, Bornstein B, Broicher A, Courtot M, Donizelli M, et~al. (2006)
  {BioModels} database: a free, centralized database of curated, published,
  quantitative kinetic models of biochemical and cellular systems.
\newblock Nucl Acids Res 34: D689-691.
\bibAnnoteFile{bib:LeNovere2006}

\bibitem{bib:Bairoch2005}
Bairoch A, Apweiler R, Wu CH, Barker WC, Boeckmann B, et~al. (2005) The
  universal protein resource ({UniProt}).
\newblock Nucl Acids Res 33: D154-159.
\bibAnnoteFile{bib:Bairoch2005}

\bibitem{bib:Borisov2009}
Borisov N, Aksamitiene E, Kiyatkin A, Legewie S, Berkhout J, et~al. (2009)
  Systems-level interactions between insulin-{EGF} networks amplify mitogenic
  signaling.
\newblock Mol Syst Biol 5: 256.
\bibAnnoteFile{bib:Borisov2009}

\bibitem{bib:Albeck2008}
Albeck JG, Burke JM, Spencer SL, Lauffenburger DA, Sorger PK (2008) Modeling a
  snap-action, variable-delay switch controlling extrinsic cell death.
\newblock PLoS Biol 6: e299.
\bibAnnoteFile{bib:Albeck2008}

\bibitem{bib:Wheeler2007}
Wheeler DL, Barrett T, Benson DA, Bryant SH, Canese K, et~al. (2007) Database
  resources of the {N}ational {C}enter for {B}iotechnology {I}nformation.
\newblock Nucl Acids Res 35: D5-12.
\bibAnnoteFile{bib:Wheeler2007}

\bibitem{bib:Yang2007paml}
Yang Z (2007) {PAML} 4: {P}hylogenetic {A}nalysis by {M}aximum {L}ikelihood.
\newblock Mol Biol Evol 24: 1586-1591.
\bibAnnoteFile{bib:Yang2007paml}

\bibitem{bib:Birtwistle2007}
Birtwistle MR, Hatakeyama M, Yumoto N, Ogunnaike BA, Hoek JB, et~al. (2007)
  Ligand-dependent responses of the erbb signaling network: experimental and
  modeling analyses.
\newblock Mol Syst Biol 3: 144.
\bibAnnoteFile{bib:Birtwistle2007}

\bibitem{bib:Brown2004}
Brown KS, Hill CC, Calero GA, Myers CR, Lee KH, et~al. (2004) The statistical
  mechanics of complex signaling networks: nerve growth factor signaling.
\newblock Phys Biol 1: 184--195.
\bibAnnoteFile{bib:Brown2004}

\bibitem{bib:Sasagawa2005}
Sasagawa S, Ozaki Y, Fujita K, Kuroda S (2005) {Prediction and validation of
  the distinct dynamics of transient and sustained {ERK} activation}.
\newblock Nat Cell Biol 7: 365--373.
\bibAnnoteFile{bib:Sasagawa2005}

\bibitem{bib:Ung2008}
Ung CY, Li H, Ma XH, Jia J, Li BW, et~al. (2008) {Simulation of the regulation
  of {EGFR} endocytosis and {EGFR}-{ERK} signaling by endophilin-mediated
  {R}ho{A}-{EGFR} crosstalk}.
\newblock FEBS Lett 582: 2283--2290.
\bibAnnoteFile{bib:Ung2008}

\bibitem{bib:Kim2007}
Kim D, Rath O, Kolch W, Cho KH (2007) {A hidden oncogenic positive feedback
  loop caused by crosstalk between {W}nt and {ERK} pathways}.
\newblock Oncogene 26: 4571--4579.
\bibAnnoteFile{bib:Kim2007}

\bibitem{bib:Singh2006}
Singh A, Jayaraman A, Hahn J (2006) {Modeling regulatory mechanisms in {IL}-6
  signal transduction in hepatocytes}.
\newblock Biotechnol Bioeng 95: 850--862.
\bibAnnoteFile{bib:Singh2006}

\bibitem{bib:Haberichter2007}
Haberichter T, Mädge B, Christopher RA, Yoshioka N, Dhiman A, et~al. (2007) {A
  systems biology dynamical model of mammalian {G}1 cell cycle progression}.
\newblock Mol Syst Biol 3: 84.
\bibAnnoteFile{bib:Haberichter2007}

\bibitem{bib:Maeda2006}
Maeda A, Ozaki Y, Sivakumaran S, Akiyama T, Urakubo H, et~al. (2006)
  {Ca$^{2+}$-independent phospholipase {A}2-dependent sustained {R}ho-kinase
  activation exhibits all-or-none response}.
\newblock Genes Cells 11: 1071--1083.
\bibAnnoteFile{bib:Maeda2006}

\bibitem{bib:Yang2007}
Yang K, Ma W, Liang H, Ouyang Q, Tang C, et~al. (2007) {Dynamic simulations on
  the arachidonic acid metabolic network}.
\newblock PLoS Comput Biol 3: e55.
\bibAnnoteFile{bib:Yang2007}

\bibitem{bib:Neves2008}
Neves SR, Tsokas P, Sarkar A, Grace EA, Rangamani P, et~al. (2008) {Cell shape
  and negative links in regulatory motifs together control spatial information
  flow in signaling networks}.
\newblock Cell 133: 666--680.
\bibAnnoteFile{bib:Neves2008}

\bibitem{bib:Hlavacek2006}
Hlavacek WS, Faeder JR, Blinov ML, Posner RG, Hucka M, et~al. (2006) Rules for
  modeling signal-transduction systems.
\newblock Sci STKE 2006: re6.
\bibAnnoteFile{bib:Hlavacek2006}

\bibitem{bib:Kiel2008}
Kiel C, Aydin D, Serrano L (2008) Association rate constants of ras-effector
  interactions are evolutionarily conserved.
\newblock PLoS Comput Biol 4: e1000245.
\bibAnnoteFile{bib:Kiel2008}

\bibitem{bib:Sunyaev2001}
Sunyaev S, Ramensky V, Koch I, Lathe~III W, Kondrashov AS, et~al. (2001)
  Prediction of deleterious human alleles.
\newblock Hum Mol Genet 10: 591--597.
\bibAnnoteFile{bib:Sunyaev2001}

\bibitem{bib:Ramsay2009}
Ramsay H, Rieseberg LH, Ritland K (2009) The correlation of evolutionary rate
  with pathway position in plant terpenoid biosynthesis.
\newblock Mol Biol Evol 26: 1045--1053.
\bibAnnoteFile{bib:Ramsay2009}

\bibitem{bib:Greenberg2008}
Greenberg AJ, Stockwell SR, Clark AG (2008) Evolutionary constraint and
  adaptation in the metabolic network of {D}rosophila.
\newblock Mol Biol Evol 25: 2537-2546.
\bibAnnoteFile{bib:Greenberg2008}

\bibitem{bib:Vitkup2006}
Vitkup D, Kharchenko P, Wagner A (2006) Influence of metabolic network
  structure and function on enzyme evolution.
\newblock Genome Biol 7: R39.
\bibAnnoteFile{bib:Vitkup2006}

\bibitem{bib:Drummond2006}
Drummond DA, Raval A, Wilke CO (2006) A single determinant dominates the rate
  of yeast protein evolution.
\newblock Mol Biol Evol 23: 327-337.
\bibAnnoteFile{bib:Drummond2006}

\bibitem{bib:Plotkin2007}
Plotkin JB, Fraser HB (2007) Assessing the determinants of evolutionary rates
  in the presence of noise.
\newblock Mol Biol Evol 24: 1113-1121.
\bibAnnoteFile{bib:Plotkin2007}

\end{thebibliography}
