PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Theor Biol. Author manuscript; available in PMC 2013 June 7.
Published in final edited form as:
J Theor Biol. 2012 June 7; 302: 29–38.
doi:  10.1016/j.jtbi.2012.02.020
PMCID: PMC3613433
NIHMSID: NIHMS365650

Quantification of Degeneracy in Biological Systems for Characterization of Functional Interactions Between Modules

Abstract

There is an evolutionary advantage in having multiple components with overlapping functionality (i.e degeneracy) in organisms. While theoretical considerations of degeneracy have been well established in neural networks using information theory, the same concepts have not been developed for differential systems, which form the basis of many biochemical reaction network descriptions in systems biology. Here we establish mathematical definitions of degeneracy, complexity and robustness that allow for the quantification of these properties in a system. By exciting a dynamical system with noise, the mutual information associated with a selected observable output and the interacting subspaces of input components can be used to define both complexity and degeneracy. The calculation of degeneracy in a biological network is a useful metric for evaluating features such as the sensitivity of a biological network to environmental evolutionary pressure. Using a two-receptor signal transduction network, we find that redundant components will not yield high degeneracy whereas compensatory mechanisms established by pathway crosstalk will. This form of analysis permits interrogation of large-scale differential systems for non-identical, functionally equivalent features that have evolved to maintain homeostasis during disruption of individual components.

Keywords: degeneracy, complexity, robustness, systems biology, biological networks

1. Introduction

In 1999, Hartwell et al. introduced the concept of modular biology, where a functional module is created by interacting molecules collectively performing a discrete function. Such modules provide an advantage for both evolvability (sensitivity to environmental changes) and robustness (insensitivity to perturbations) in a system [1]. Reconciling the divergent requirements of evolvability and robustness requires knowledge regarding the connections between these functional units and an appreciation for how each module integrates information arriving from multiple inputs. In complex biological systems such as neural networks, there has been recent emphasis on features of structural complexity such as degeneracy. As first introduced in [2], structural complexity can be understood in terms of the interplay between specialization of functions in individual modules (functional segregation) and the ability of the modules to interact and perform functions coherently (functional integration). A highly complex system maintains segregation of function while still allowing for functional integration. Degeneracy measures how well functionally independent modules can interact to produce redundant outputs. Modules that are structural duplicates form a completely redundant system and always produce the same output; degenerate systems arise from structurally distinct modules with different outputs interacting to produce the same output under certain conditions.

Systemic features like degeneracy, complexity, robustness are related to one another. It has already been observed via numerical simulations for neural networks that high degeneracy not only yields high robustness, but also it is accompanied by an increase in structural complexity [3]. Thus, it is believed that degeneracy is a necessary feature in the evolution of complex biological systems, partly because genetically dissimilar organisms will drive toward convergence of function while maintaining non-redundant components [4]. It is also believed that the robustness and adaptability that ensue from degeneracy are key features of complex biological systems at multiple scales [5] and organisms have evolved to contain many non-identical structures to produce similar functions. In this manner, degeneracy helps to fulfill the necessary properties of biologically functional modules. While the concept of degeneracy was introduced for neural networks, there is strong evidence that many dynamical biological networks including cellular metabolic and signaling networks exhibit the property of degeneracy; for example, protein kinase isoforms regulated by separate genes phosphorylate the same substrate protein [4]. Furthermore, because regulatory features of protein or metabolic networks often rely on dynamical systems analysis, such theory must be compatible with kinetic description of biochemical reactions rather than neural network [3] or logic-based description [6]. Although some features like regulation and robustness of biochemical networks of signal transduction have been studied quantitatively [7, 6], the features of interest here, such as degeneracy and complexity, have not been formalized mathematically in terms of ordinary differential equation description that is so commonly used for describing protein networks. By developing a method for calculating degeneracy in general biological networks modeled by differential equations, network structures can be explored that fulfill the contrary requirements of evolvability and robustness first posited by Hartwell et al. Understanding these features becomes important as increasingly complex biological systems are being designed ab initio in synthetic biology. In this article, we first mathematically formalize the relationships between degeneracy, complexity and robustness in dynamical systems. We next establish that high degeneracy always yields high complexity. Finally, we illustrate two distinct approaches for numerically calculating degeneracy in dynamical systems, a predator-prey model and a kinetic signaling model. In many instances, the metric of degeneracy is a useful indicator of how easily the system adapts under evolutionary pressure.

2. Mathematical approach

2.1. Random perturbations of ODE system

Our strategy is to inject a fixed amount of stochastic perturbation into a differential system. With such small random perturbations, the corresponding variable sets of modules of the network become stochastic processes. If two modules have strong functional connectivity, these two stochastic processes should have high statistical correlation. Conversely, two functionally independent components must be statistically independent. This statistical connectivity can be measured by the mutual information of the components. It is already known that degeneracy measures the ability of structurally different components to perform the same function, while complexity measures the degree of functional integration and segregation between different components. Using these ideas, we can quantify degeneracy and complexity using linear combinations of mutual information.

To avoid mathematical challenges from injecting small stochastic perturbations into an ODE system, we will assume throughout this work that the differential system is dissipative. In other words, asymptotically the differential system will gravitate to its global attractor such that the generated stochastic process has a stable invariant measure. This invariant measure allows us to obtain the asymptotic correlation between modules. The random perturbation and its invariant measure is described below.

x=f(x),x[set membership]Rn
(1)

i.e., we consider the Ito stochastic differential equation (SDE)

dX=f(x)dt+εσ(x)dWt,
(2)

where Wt is the Wiener process, σ is a non-singular, n × n matrix-valued function and ε is a small parameter. The time evolution of the probability density function associated with the SDE (2) satisfies the so-called Fokker-Planck equation

ρt=12ε2i,j=1n(Aijρ)ij[nabla](fρ),
(3)

where A(x) = σ(x) σT(x) is an n × n symmetric nonnegative definite matrix.

Of particular importance among the solutions of the Fokker-Planck equation are the steady states, which satisfy the stationary Fokker-Planck equation

{12ε2i,j=1n(Aijρ)ij[nabla](fρε)=0,ρ(x)>0,Rnρε(x)dx=1.
(4)

A smooth solution ρε of the stationary Fokker-Planck equation (4) is known to uniquely exist [8, 9] if

  • f is differentiable and σ is twice differentiable on Rn; and
  • there exists a Lyapunov function V (x) > 0, V → +∞ such that 12ε2i,j=1nAij(x)[partial differential]ij2V(x)+f(x)·V(x)<γ for some constant γ > 0 and all |x| sufficiently large.

We remark that while the ODE (1) may have many complicated invariant measures without even having density functions, the steady-state of the SDE (2) is nevertheless unique and smooth. We also note that when ε is fixed, we denote the invariant solution as ρ instead of ρε.

2.2. Definitions of Degeneracy, Complexity and Robustness

2.2.1. Degeneracy and Complexity

Inspired by, but differing from [3], our definition of degeneracy is divided into several steps. In the first step, we define projected density, entropy and mutual information associated with any subspace. Then we fix a subspace as the “output” set and define its associated degeneracy by considering the complementary subspace as the “input” set. Lastly, we define the degeneracy of the entire system by varying the output sets and taking the maximum among all degeneracies of these sets. efinitions of projected density, entropy and mutual information are provided in Appendix A. In a biological network, mutual information between two components I1 and I2, MI(I1; I2), measures the functional connectivity between the components. Using mutual information, degeneracy can be defined as follows.

Let An external file that holds a picture, illustration, etc.
Object name is nihms365650ig1.jpg be a fixed subspace of Rn, viewed as an output set. We denote I as the complementary subspace to An external file that holds a picture, illustration, etc.
Object name is nihms365650ig1.jpg, viewed as the input set. In other words, the set An external file that holds a picture, illustration, etc.
Object name is nihms365650ig1.jpg is a fixed set of “observables” when the system (2) is excited by noise. To measure the impact of noise on all possible components of the input set, we consider any subspace Ik of I and denote its complementary set in I by Ikc. The interacting information among Ik, Ikc and An external file that holds a picture, illustration, etc.
Object name is nihms365650ig1.jpg is defined by

D(k)=MI(I;Ik;O)=MI(Ik;O)+MI(Ikc;O)MI(I;O).
(5)

The interacting information measures how much Ik and Ikc are structurally different but perform the same function as signified by the output set An external file that holds a picture, illustration, etc.
Object name is nihms365650ig1.jpg.

We note that unlike the mutual information between two subspaces, the interacting information among three subspaces can take negative values ([10]).

Similar to the case of neural networks, we define the degeneracy associated with An external file that holds a picture, illustration, etc.
Object name is nihms365650ig1.jpg by averaging all the interacting information among all possible subspaces of I, i.e.,

D(O)=left angle bracketMI(I;Ik,O)right angle bracket=Ik12Cknmax{MI(I;Ik;O),0}
(6)

Similar to degeneracy, complexity C( An external file that holds a picture, illustration, etc.
Object name is nihms365650ig1.jpg) associated with An external file that holds a picture, illustration, etc.
Object name is nihms365650ig1.jpg could be obtained by averaging all the mutual information between Ik and Ikc, i.e.,

C(O)=left angle bracketMI(Ik;Ikc)right angle bracket=Ik12CknMI(Ik;Ikc).
(7)

This value measures how much codependency in a network appears among different modules rather than different elements (units that constitute a module).

Now, for a fixed diffusion matrix σ and ε > 0, we define the degeneracy Dε,σ and structural complexity Cε,σ of the system (1) as

Dε,σ=MaxOD(O),Cε,σ=MaxOC(O),

We call a differential system (1) degenerate (resp. complex) with respect to a diffusion matrix σ if there exists ε0, such that Dε,σ>0 (resp. Cε,σ>0) for all 0 < ε < ε0.

We would like to make the following remarks

  • In many applications, one can often choose σ(x) as the identity matrix, so that the noise perturbation becomes purely white. But a variable diffusion matrix σ(x), associated with a colored noise perturbation, should play an important role in detecting the key output set mainly responsible for the degeneracy.
  • For a particular biological system, one often has a natural choice of “observable” variables to be used as the output set An external file that holds a picture, illustration, etc.
Object name is nihms365650ig1.jpg. If one can select a special subspace Ik0 of the complementary subspace I so that the interacting information MI(I;Ik0; An external file that holds a picture, illustration, etc.
Object name is nihms365650ig1.jpg) among the three is positive with respect to a fixed diffusion matrix, then it follows from the definition that the whole system has a certain level of degeneracy. Since the interacting information could be negative, we take the average of max{MI(I; Ik; An external file that holds a picture, illustration, etc.
Object name is nihms365650ig1.jpg), 0} to measure the degeneracy of the system.

2.2.2. Robustness – System robustness and functional robustness

Our notion of robustness will be defined in a way that reflects the strength of attraction of the global attractor of system (1). Recall that the system (1) was assumed to be dissipative so that a global attractor already exists. We denote the global attractor by An external file that holds a picture, illustration, etc.
Object name is nihms365650ig4.jpg.

To define the robustness, we require in this paper that An external file that holds a picture, illustration, etc.
Object name is nihms365650ig4.jpg is a strong attractor in the following sense. The attractor An external file that holds a picture, illustration, etc.
Object name is nihms365650ig4.jpg is said to be a strong attractor with nonnegative index α if there exists a compact neighborhood N with C1 smooth boundary and a Lyapunov function V (x) such that

[nabla]V(x)f(x)αdist(x,A),forallx[set membership]N.

For a strong attractor An external file that holds a picture, illustration, etc.
Object name is nihms365650ig4.jpg, the system robustness of An external file that holds a picture, illustration, etc.
Object name is nihms365650ig4.jpg is the following quantity

R=inf{1α:αisanindexofA}.

The system is said to be robust if An external file that holds a picture, illustration, etc.
Object name is nihms365650ig4.jpg is a strong attractor and R is finite.

If the performance function p(x) of the system is given with the following assumptions: (1) p(x) = 1 ∀x [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms365650ig4.jpg; (2) 0 < p(x) < 1 if x [negated set membership] An external file that holds a picture, illustration, etc.
Object name is nihms365650ig4.jpg, then, following Kitano [7], one can define functional robustness Rf (ε) as

Rf(ε)=ρε(x)p(x)dx.

Using this notation, the system has the best performance when the perturbation vanishes and robustness is interpreted as the ability to preserve phenotype rather than maintaining a fixed steady state. We remark that if a system is robust and some property of the performance function is also known, then some estimate on the functional robustness can be made.

2.3. Connection between degeneracy, complexity and robustness

Degeneracy, complexity and robustness are not isolated concepts. More and more examples suggest some internal connections among them (see [4, 5]). In fact, these relationships could not only be observed in biological experiments, but also be verified by our mathematical quantification of degeneracy, complexity and robustness. From the definition, degeneracy always implies complexity. Further, with some additional conditions, a robust system must have certain level of degeneracy (and hence, also complexity).

2.3.1. Degeneracy and Complexity

It has been observed in neural networks that a degenerate system must have a complex structure [3]. The same property holds for a differential system, which is commonly used to describe biochemical networks. A simple calculation shows that

MI(I;Ik;O)min{MI(I;Ik),MI(Ikc;O),MI(I;O)}.
(8)

If we compare equations (6) and (7), by taking the average among all possible subsets Ik, we obtain

C(O)D(O).

because MI(I;Ikc;O)MI(I;Ikc). In other words, with respect to a fixed diffusion matrix, degeneracy implies complexity.

This explains the observation in [4] that biological systems selected for high degeneracy are accompanied by high complexity.

Remark

We can prove equation (8) in the following way:

MI(X;Y;Z)=H(X)+H(Y)+H(Z)H(X,Y)H(Y,Z)H(X,Z)+H(X,Y,Z)=H(X)+H(Y)H(X,Y)(H(X,Z)+H(Y,Z)H(Z)H(X,Y,Z))=MI(X;Y)MI(X;Y[mid ]Z)

Since the l mutual information is nonnegative: MI(X; Y |Z) ≥ 0, we have MI(X; Y; Z) ≤ MI(X; Y ). (The non-negativity of conditional mutual information is a direct corollary of Kullback’s inequality, or see [11])

Similarly we can prove MI(X; Y; Z) ≤ MI(X; Z) and MI(X; Y; Z) ≤ MI(Y; Z), from which (8) follows.

2.3.2. Degeneracy and Robustness

We would like to examine the connections between degeneracy and robustness for an ODE system (1). Robustness alone does not necessarily imply degeneracy of the system; this is because one can certainly have a system with zero complexity (e.g., a system with many symmetric components) which is however robust. By (8), such a system must be non-degenerate. Therefore, for a robust system to be degenerate, the system must be complex and such structural complexity often gives rise to some kind of embedding complexity of the global attractor into the phase space. Roughly speaking, the components of a complex system interact strongly with each other and as a result, the global attractor is twisted in the phase space such that it does not lie in any hyperplane. To characterize the twist property of the global attractor, it is natural to consider its projections on certain hyperplanes and measure the dimensions of the corresponding projections. We note that the attractor as well as its projections may only be fractal sets, hence they should be measured with respect to the Minkowski dimension, also called box counting dimension [12].

For a subspace An external file that holds a picture, illustration, etc.
Object name is nihms365650ig5.jpg of Rn, we denote by dV the co-dimension of An external file that holds a picture, illustration, etc.
Object name is nihms365650ig4.jpg in An external file that holds a picture, illustration, etc.
Object name is nihms365650ig5.jpg, i.e., the dimension of An external file that holds a picture, illustration, etc.
Object name is nihms365650ig5.jpg subtracts the Minkowski dimension of the projection of An external file that holds a picture, illustration, etc.
Object name is nihms365650ig4.jpg to An external file that holds a picture, illustration, etc.
Object name is nihms365650ig5.jpg.

The twisted attractor is defined as follows. The global attractor An external file that holds a picture, illustration, etc.
Object name is nihms365650ig4.jpg is said to be twisted if there is a linear decomposition Rn = An external file that holds a picture, illustration, etc.
Object name is nihms365650ig6.jpg [plus sign in circle] An external file that holds a picture, illustration, etc.
Object name is nihms365650ig7.jpg [plus sign in circle] An external file that holds a picture, illustration, etc.
Object name is nihms365650ig1.jpg such that

dI+dJ+dO+dRn<dI[plus sign in circle]J+dI[plus sign in circle]O+dJ[plus sign in circle]O.

We have the following theorem

The invariant probability density function ρε is said to be regular for An external file that holds a picture, illustration, etc.
Object name is nihms365650ig4.jpg if there exists some function C(K) > 0 that is independent with respect to ε, such that

min(ρε(x))Cmax(ρε(x));xwithdist(x,A)Kε

for all 0 < ε < ε0 and K > 0.

Theorem 1

If the system (1) is robust with a twisted global attractor, and if the ε-invariant density function ρε is regular for An external file that holds a picture, illustration, etc.
Object name is nihms365650ig4.jpg, then there exists an ε0 > 0, such that An external file that holds a picture, illustration, etc.
Object name is nihms365650ig2.jpg > 0 for all 0 < ε < ε0

For a proof of the theorem see Appendix B.

2.3.3. Degeneracy at equilibrium

Degenerate behavior could occur not only at the twisted attractor, but also at certain equilibria, or what a biologist may regard as homeostasis. Here, we introduce another theorem on the connection between robustness and degeneracy. If an ODE system has a unique equilibrium point and in the neighborhood of this equilibrium point the reactions to random perturbations have certain level of diversity, then we claim that it is a degenerate system. More precisely, if different directions demonstrate different sensitivities under random perturbation, then it is a degenerate system. Since it is known that a large number of chemical reaction networks have unique stable equilibrium points, the degeneracy near equilibrium may be more applicable for biological reaction networks.

Assume that system (1)

x=f(x)

has a unique stable fixed point, say x0. Let B denote the Jacobian matrix of f(x) at x0. Since we have assumed the robustness already, it is obvious that the eigenvalues of B only have negative real parts. After some calculation, one can find the solution to the stationary Fokker-Planck equation (4):

ρ=1KezTS1z/2ε+o([mid ]z[mid ]2),
(9)

where z = xx0. The symmetric positive definite matrix S solves the Lyapunov equation uniquely:

SBT+BS+A=0

where A = σ(x0) σT(x0).

With the stationary solution ρ, we can find the marginals on target subspaces. It is known that the marginal of a normal distribution is also normal, whose covariance matrix is the corresponding sub-matrix of S. More precisely, if X = span{xa1, · · ·, xak} is a subspace, then the sub-matrix S(a1, · · ·, ak; a1, · · ·, ak) is the covariance matrix of the projection of ρε on subspace X. For simplicity, we denote S(a1, · · ·, ak; a1, · · ·, ak) as S(X).

Then we can compute the degeneracy with split X = I1 [plus sign in circle] I2 [plus sign in circle]O. Since equation (9) approximates a multivariate normal distribution, calculation of degeneracy yields the following theorem: if

Γ:=log[mid ]S(I1)[mid ][mid ]S(I2)[mid ][mid ]S(O)[mid ][mid ]S(X)[mid ][mid ]S(I1,I2)[mid ][mid ]S(I1,O)[mid ][mid ]S(I2,O)[mid ]>0.
(10)

then the system is degenerate.

In fact, it can be shown that as ε → 0, the degeneracy of ρε with respect to decomposition I1, I2, O converges to Γ.

Two approaches can be taken for calculating D for a coupled differential system. One relies on Monte Carlo simulations (Appendix C and Algorithm 1 in Fig. 1), while the other is based on stochastic analysis by the Freidlin-Wentzell quasi-potential method [17] (Appendix C and Algorithm 2 in Fig. 1). We demonstrate the utility of each with biological examples below.

Figure 1
Algorithms for calculating degeneracy when fixed points are unknown (Algorithm 1) or known (Algorithm 2).

3. Illustrations

3.1. Implications of degeneracy in a signal transduction pathway

Consider a simple example consisting of three modules A, B and C shown in Fig. 2. A and B serve as inputs while C is the output. If module A has a functional relationship with the output module C, the mutual information between the two is high. Similarly, if modules B and C share high mutual information, they are functionally related as well. However, both modules A and B being functionally related to the output is not enough for degeneracy. We also require A and B to be structurally different. This can be checked by treating A and B as a single unit, measuring its mutual information with the output and comparing it with the mutual information A and B share individually with C. The value MI(A; C) + MI(B; C) − MI({A, B}; C), thus measures the degeneracy, or how much more correlation the inputs A and B share with the output C than expected. Defining degeneracy enables us to explore the applications of degeneracy quantitatively.

Figure 2
A toy example of a modular biological network

Using a simplified model of crosstalk in protein signal transduction, we illustrate the calculation of degeneracy using Algorithm 2 (Fig. 1 and Appendix C). We also demonstrate how certain biological features of the signaling network affect the numerical value of degeneracy.

For illustrative purposes we have chosen the JAK-STAT signaling pathway since this system presents features that are useful for illustrating the concept of degeneracy. The JAK-STAT pathway is a two-step intracellular signaling pathway in which a member of the JAK family of kinases, typically bound to a transmembrane receptor, is activated by phosphorylation following ligation of the receptor with extracellular cytokine. The activated JAK molecule phosphorylates STAT which can then dimerize and act as a transcription factor. The signaling pathway is regulated by several mediators including phosphatases that dephosphorylate JAK and STAT molecules, thereby inhibiting their catalytic activity [19]. It has been shown previously that cytokine receptor activation can be accompanied by production of reactive oxygen species (ROS) in response to multiple kinds of cytokines [20]. The generated ROS reacts with some phosphatases to oxidize them reversibly, resulting in temporary inactivation of the phosphatases. Phosphatase inactivity results in amplification of STAT phosphorylation. Sharma et al. demonstrated that different cytokines, signaling through their respective receptors, can crosstalk in an ROS-mediated manner to amplify signals coming through other cytokines [20]. This is a result of oxidative inactivation of phosphatases regulating the different JAK-STAT pathways.

Based on this information we have constructed a simplified model of crosstalk between IL-4 and Epo signaling. IL-4 signals through the IL-4 receptor and activates the JAK3/STAT6 pathway [21]. Epo signals through the Epo receptor and activates the JAK2/STAT5 signaling pathway [22]. Multiple phosphatases can regulate these pathways. To illustrate the crosstalk between the pathways we have chosen one phosphatase, PTP1B, which is important in both signaling pathways. In the IL-4 pathway it directly dephosphorylates STAT6 whereas the substrate of PTP1B in the Epo pathway is JAK2 [23, 24]. PTP1B is also susceptible to ROS-mediated oxidative inactivation [20]. This information was compiled to get the IL-4/Epo crosstalk model shown in Fig. 3A. For the sake of parsimony we have treated phosphorylated STAT as the output and ignored phospho-STAT dimerization. The details of model implementation are provided in Appendix E and supplementary information.

Figure 3
Illustration using IL-4R and EpoR crosstalk model. A. The core JAK-STAT modules of IL-4R and EpoR pathways with crosstalk. Both modules are regulated by PTP1B, both generate ROS which oxidatively inactivates PTP1B. B. Crosstalk enhances degeneracy. The ...

The model we constructed was used to empirically study relationships between the signaling pathway and the computed degeneracy. We chose the receptors (IL-4R and EpoR) as a pair of inputs to the system and activated STAT molecules (STAT5* and STAT6* in Fig. 3A) as the output. Using the method outlined in Algorithm 2 (Fig. 1), the model in Fig. 3A was found to be degenerate with a value of D equal to 0.4267. Since our theoretical results relate increased complexity with increased degeneracy, we sought to verify if this was reflected in our model of JAK-STAT crosstalk. The cross-talk between the two linear JAK-STAT pathways is the source of increased complexity of the system. To reduce the complexity of the system, we abrogated all cross-talk by switching off ROS production and regulation by the common phosphatase PTP1B to get the independent signaling systems shown in Fig. 3B. The calculated value of degeneracy decreased by more than 99% for this system as compared to the pathway in Fig. 3A and the value of D was calculated to be 0.0016. This demonstrates that cross-talk between signaling pathways results in increased complexity which could result in increased degeneracy.

Redundancy in signaling systems can also lead to complexity in the pathway, in the sense that there can be significant amount of cross-talk between parallel pathways. However, a redundant system is by definition not degenerate because the redundant modules perform identical functions under any given condition. To test how a redundant system compares with a degenerate system, we modified the pathway in Fig. 3A to that shown in Fig. 3C by inserting some hypothetical connections. This was done to ensure that the two modules were structurally identical and affected the output (STAT5* and STAT6*) identically. The rate parameters were also identical for the two modules resulting in a completely redundant system where EpoR and IL-4R affect STAT5 and STAT6 identically. The redundant system was found to still have a positive D but the magnitude was reduced by more than 68% as compared to the value calculated for the system in Fig. 3A. This agrees with the understanding that redundancy does not lead to degeneracy and our calculation of D successfully reflects this.

3.2. Degeneracy in a Lotka-Volterra system

We provide a three dimensional example to demonstrate how to verify degeneracy using the Monte Carlo method (see Appendix C). Consider the following competitive Lotka-Volterra system

x.1=x1(3x1x2x3)x.2=x2(4x1x22x3)x.3=x3(7.2212.61x11.611x23x3).

This system represents a simple three-species competitive population model. The system has a limit cycle as described previously (Fig. 4) [18]. Using the theory of quasi-potential functions, one can rewrite the vector field as −[nabla]Ψ(x1, x2, x3) + l(x1, x2, x3), where Ψ is called a quasi-potential function and l is a small perturbation in a definite sense with [nabla]Ψ · l = 0. It is well known that for such a system admitting a limit cycle, Ψ is a Lyapunov function which is as regular as the vector field. It then follows from definition that the system is robust. Furthermore, the condition in Theorem 1 is also satisfied due to the regularity of the quasi-potential function.

Figure 4
Limit cycle of the Lotka-Volterra system showing a twisted attractor.

Numerical simulations show that the limit cycle is not parallel to any coordinate axis. In fact, it follows that dx = dy = dz = 0, dxyz = 2, dxy = dxz = dyz = 1. Hence the attractor is also twisted. Now applying the theorem on twisted attractors, we conclude that the system is degenerate. Further details regarding this illustration are provided in Appendix D.

3.3. Degeneracy enhances evolvability

It has been argued that not only is degeneracy an outcome, but also an important driver of evolution [4]. We performed computational simulations of adaptive evolution to study the interplay between degeneracy and evolvability (see Appendix F for details). We use the term evolvability to generally represent the ease with which a biological network can adapt to an environmental change through the process of evolution to increase its fitness. We studied a network modeled using ODEs, consisting of a fixed number of nodes and evolving by mutating the strengths of connections between nodes. Using empirical studies with stochastic simulations, we observed that evolution was often, but not always, accompanied by increase in degeneracy. However, systems with higher initial degeneracy exhibited greater evolvability. Since adaptive evolution is random and does not follow any design, it is to be expected that not all mutations that increase fitness will also lead to increased degeneracy. On the other hand, a system with high initial degeneracy consists of multiple interwoven modules performing the same function in different ways. As a consequence, mutations in each of these modules allow the system to evolve a different phenotypic response resulting in an enhanced ability to explore the phenotype space. We note that our observation that evolution does not always lead to increased degeneracy may not be true in the longer term. Since systems with greater degeneracy have an increased ability to adapt to environmental changes, mutations leading to increased fitness in a particular environment but reduced degeneracy, lose out on their ability to respond to further environmental changes. In some sense, degeneracy can be thought of as a measure of the evolvability of a system. This means that the ability to quantify degeneracy gives us a sense of the evolvability of networks, at least on a relative scale, simply by analyzing its structure.

4. Discussion

Degeneracy can be generally understood as the ability of structurally distinct components of a system to behave similarly under certain conditions, while the behavior may be different under other conditions. Increasingly large numbers of instances of degeneracy are being found in biological systems at all scales ranging from molecular to animal population levels [4]. Particularly, in the context of cellular signaling networks, there are multiple examples of degenerate behavior. Different members of the interleukin (IL) family can activate the same transcription factor. For instance, IL-2, IL-7 and IL-21 can all activate STAT1 [25]. Growth factors can bind to multiple types of receptors in the EGF receptor family [26]. MAPK signaling induced by growth factors and stress exhibits promiscuous interaction between MEKK and MAPK proteins where multiple types of MEKKs can activate the same MAPK and a single MEKK can activate multiple MAPKs [27]. Recently, experimental studies have indicated that a significant role for genetic buffering by non-homologous genes (i.e. functional redundancy or degeneracy) [28] exists and may confer a selective advantage over paralogs for regulation. The widespread appearance of degeneracy in biological systems across scales suggests that it is a property favored by adaptive evolution. Because of the nature of evolutionary systems, desired changes are not intelligently engineered into them. Instead, they evolve by incorporating random changes some of which improve the fitness for survival. Also, since no part of an adaptive biological system is outside the scope of evolution, multiple components of the system can evolve to achieve the same kind of adaptation. It is, therefore, reasonable to expect an evolutionary system to increase in complexity over time which could result in enhanced degeneracy. Thus, the theory presented here provides a method not only for looking at structural characteristics of biological networks, but also for exploring the links between degeneracy and evolvability in biological systems modeled using ODEs.

Complexity arising as a consequence of the evolutionary process means that biological systems rarely, if at all, operate in isolation. Chen et al. showed that the behavior of a signaling pathway in isolation is different from the behavior it exhibits when put in the context of a more complex intracellular environment [26]. Systems biologists are aware that cellular signaling pathways which are classically seen as isolated, and often linear, chains of biochemical modifications rarely operate in this simple fashion. The connections between signaling pathways give rise to networks with much greater complexity. These resulting systems can exhibit degenerate behavior in that one signaling pathway, a module by itself, interacts with another structurally different module with both of them regulating the same output, resulting in similar or different outcomes depending on conditions. This is important from the point of view of applications such as drug targeting. For instance, despite major efforts, very few drugs specifically targeting the PI3K signaling pathway, which exhibits strong crosstalk with a number of other pathways, make it to the clinical trial stage [29]. The complexity arising from cross-talk is thought to be one reason for the failure of specific inhibitors to work successfully in cells. Determining what points should be targeted in a complex signaling network is a critical question for drug design. It is therefore desirable that the extent of compensation between connected pathways be defined quantitatively. A quantitative measure of degeneracy in the network can be exploited to identify candidate points in the network most suitable for drug intervention. Quantification of degeneracy can also have applications in synthetic biology for designing system modules that are structurally distinct but can be made to perform similar functions when needed.

Like degeneracy, complexity of a system can also be quantified using mutual information between modules. A complex system is one with high overall mutual information between different modules (and not between simpler elements making the modules). For example, in the context of crosstalk between the IL-4 and Epo pathways, the IL-4 and Epo pathways separately can be thought of as modules while the individual molecules constitute the basic elements of the system. This means that a complex system maintains a modular structure and also has co-dependence between the existing modules. Maintenance of functional modules means the system is functionally segregated. At the same time, co-dependence between the modules means there is integration of function between the modules. A complex system, therefore, preserves both these properties which is captured by the measure of complexity we have presented.

Using a model of cross-talk in interleukin signaling we have demonstrated the biological significance of this numerical measure of degeneracy. The IL-4R and EpoR signaling pathways, by virtue of phosphatase and ROS mediated crosstalk, give rise to a degenerate system when the receptors are treated as input and STAT activation as output. By computationally manipulating the signaling pathway we have empirically shown the relationships between degeneracy and some network features. As first demonstrated by Tononi et al. for neural networks, we found that independent signaling modules exhibit very low degeneracy [3]. Reduced connectivity also means that the modules have a reduced ability to influence each other via crosstalk, which our definition of degeneracy is able to reflect. Simply increasing the complexity of the network may not be sufficient to guarantee a degenerate system. For instance, fully redundant signaling modules with a high degree of crosstalk result in a structurally complex system. In our computational analysis, when the modules were made fully redundant by making them identical in the structure of the network and in the strengths of the internal connections, the calculated degeneracy dropped despite an increase in network edges. This agrees with the notion that redundancy and degeneracy are functionally distinct, and demonstrates that our definition of degeneracy is able to distinguish between a truly degenerate system against one with high complexity but low degeneracy.

We also present a definition of robustness in the context of differential equation models of biological systems. The stability of a differential system can be measured by its robustness under random perturbation. A robust system strongly resists change under fixed random perturbation. Moreover, as suggested by [7], if we know the performance function of a system, we don’t have to require that the system offer this resistance everywhere – the system only needs to be stable at places where the performance function decreases dramatically. Biologically, this means that a robust system is not necessarily one that is able to maintain a fixed steady state; instead it is a system that is able to maintain its phenotype in the face of perturbations [7]. We have provided a definition of functional robustness that takes this into account. Further, we have shown that robustness and degeneracy are connected– a robust system has positive degeneracy when it satisfies certain conditions (Theorem 1).

While these illustrations with simplistic biological models provide some insight into the significance of our theoretical framework for defining degeneracy, several aspects remain to be explored. For instance, does the calculated degeneracy provide an estimate of the ability of cross-talking pathways to compensate for each other under perturbation? System dynamics are of great importance in understanding cellular signaling networks. Our method for calculating degeneracy takes into account only the fixed points of the differential system. In thinking about the meaning of calculated degeneracy in the context of cell signaling, it is important to keep system dynamics in mind. The outcome of a signaling event is not always dictated by the steady state value, instead instantaneous rates of changes or integrated values of signals may be of relevance in a given system. For this reason it is important to explore the relationships between system dynamics and degeneracy. Given the “no free lunch” concept in control systems in which operating performance of one control function comes at the cost of fragility elsewhere [30, 31], the consequences of degenerate network properties over redundant components can be explored further. These concepts may be exploited in the design of synthetic biological circuits to ensure a desired functional outcome under a variety of biological contexts. Although several issues remain to be addressed, the methods presented in this paper are significant in providing a theoretical framework to the concept of degeneracy and functional robustness for the class of systems represented by differential equations.

Highlights

  1. Methods to quantify degeneracy and robustness in differential systems.
  2. Cross talk between signaling pathways enhances degeneracy.
  3. Redundant signaling systems show low degeneracy.
  4. Robustness implies degeneracy under certain conditions.

Supplementary Material

01

Acknowledgments

WH was supported by NSFC, Fok Ying Tung Education Foundation, FANEDD (Grant 200520) and the Fundamental Research Funds for the Central Universities. MLK was supported by a NIH New Innovator Award DP2OD006483. YY was supported by NSF grants DMS0708331, DMS1109201 and a scholarship from Jilin University.

Appendix A. Definition of Projected Density, Entropy and Mutual Information

Let An external file that holds a picture, illustration, etc.
Object name is nihms365650ig5.jpg be the variable set of equation (1). Biologically An external file that holds a picture, illustration, etc.
Object name is nihms365650ig5.jpg means the set of elements or species of the network.

Let ρ be a smooth solution of (4) for fixed ε and σ. For any subspace I of Rn coordinated by u [set membership] I, we define the marginal distribution with respect to I by

ρI(u)=Jρ(u,v)dv,

where J is the complementary subspace of I coordinated by v [set membership] J. The coordinates of I are a subset of the variable set An external file that holds a picture, illustration, etc.
Object name is nihms365650ig5.jpg, so biologically I represents a subset of the whole network. For instance, in R3 = {(x1, x2, x3)} if I = {(0, u, 0)} and J = {(v1, 0, v2)}, then u = x2 and

ρI(u)=Jρ(x1,x2,x3)dx1dx3.

The projected entropy associated with the projected density above is defined by

H(ρI)=IρI(u)logρI(u)du.

For any two subspaces I1 and I2, the direct sum I = I1 [plus sign in circle] I2 is also a sub-space. We then define their joint entropy H(I1, I2) simply by the projected entropy H(I1 [plus sign in circle] I2) associated with the direct sum, i.e.,

H(I1,I2)=H(I1[plus sign in circle]I2)=I1[plus sign in circle]I2ρI1,I2(u,v)logρI1,I2(u,v)dudv,

where

ρI1,I2(u,v)=Jρ(u,v,w)dw

with J being the complementary subspace of I1 [plus sign in circle] I2. The mutual information among subspaces I1, I2 is defined by

M(I1;I2)=H(I1)+H(I2)H(I1,I2).

It is easy to see that

MI(I1;I2)=I1[plus sign in circle]I2ρI1,I2(u,v)logρI1,I2(u,v)ρI1(u)ρI2(v)dudv
(A-1)

Statistically, the mutual information (A-1) measures the correlation between marginal distributions with respect to subspaces I1 and I2.

Appendix B. Proof of Theorem 1

Proof

This Theorem is a corollary of the Entropy-Dimension identity proved in [13]. Under the given conditions, we have

limε0H(ρε(x))logε=Nd
(A-2)

where H(ρ) means the entropy of ρ. Then using the definitions of degeneracy and a twisted attractor, we can prove the positivity of the degeneracy Dε,σ.

Remark

We note that ρε is always regular if there exists a quasi-potential function W(x) of An external file that holds a picture, illustration, etc.
Object name is nihms365650ig4.jpg, such that for every 0 < ε < ε*, we have

ρε(x)=1KeW(x)/ε2+o(ε)

where

K=RNeW(x)/ε2dx

and o(ε) means high order terms of ε.

From [14, 15], we can find the desired function W(x) whenever the Freidlin-Wentzell quasi-potential function W(x) has second order derivatives. From [15, 16], we know that the Freidlin-Wentzell quasi-potential function W(x) has high regularity in the neighborhood of stable nodes and limit cycles. Thus the example discussed below in 3.2 satisfies the conditions of Theorem 1. For more detailed introduction of Freidlin-Wentzell quasi-potential function, see [17].

Appendix C. Calculating degeneracy

The Monte Carlo method

According to the definitions previously provided in equations (5) and (7), the degeneracy and complexity can be computed for a general ODE system if we can calculate the mutual information between two components. We used the following way to calculate the mutual information (Algorithm 1 in Fig. 1). First, a rough bound of the attractor was determined numerically. This was done using a simple Monte-Carlo simulation with some statistics. Using the Monte-Carlo simulation we obtained a sample set of solutions of equation (2) by randomly choosing a set of points S in the space and letting it evolve with equation (2) until some large enough time T. Assuming {a1, · · ·, aN } is a sample of variable x1, and μ and σ are the mean and standard deviation of the sample, [μ − 3σ, μ + 3σ] was chosen as a rough bound of the attractor. The rough bounds of the other variables were determined similarly.

Then we generated a numerical grid in the rough bound of the attractor and ran another Monte-Carlo simulation to create another large sample. This sample was required to be large enough such that an approximate probability distribution could be computed numerically using the sample set. With the approximate probability distribution functions available, numerical integration over target variables was used to calculate entropy, mutual information and interacting information. Note that degeneracy is the interacting information and complexity is the mutual information.

The Lyapunov method

For most ODE systems generated from chemical reaction networks, there exists a unique stable equilibrium. The invariant measure for such a system can be obtained using some linear algebra calculations (Fig. 1). We first calculated the steady-state solution x0. The Jacobian matrix B was then obtained either numerically or analytically. The invariant measure was approximated with a multivariate normal distribution with covariant matrix S, where S solves the Lyapunov equation

SBT+BS+A=0

S could be solved analytically or numerically. Several software are available for solving the Lyapunov equation numerically. The degeneracy was then obtained using equation (10).

Appendix D. Details of implementation of the twisted attractor illustration

We have applied Monte-Carlo simulation to compute the degeneracy for a fixed value of ε. The approach, as represented in Algorithm 1 (Fig. 1) is to first use the Monte-Carlo method to obtain the invariant measure, then to project the measure onto relevant subspaces to compute their entropies. In the Monte-Carlo simulation, we set noise matrix δ as identity. The degeneracy is a linear combination of the entropies. For example, when ε = 0.001, the degeneracy is computed as D = 2.3144. Note that the calculated D is the degeneracy with respect to subspaces coordinated by x1, x2, x3, which is less than the degeneracy of this dynamical system when we take the maximum over all possible subspace splits.

This number will increase when ε decreases. The accuracy of the Monte-Carlo method is of order N−1/2, where N is the total number of grid points used, the accuracy of integration is compatible to the grid size. In our simulation, we have taken N = 4000000 and a grid size of 0.0005. This gives an accuracy of around 10−3 for our computation.

Appendix E. ODE model of IL4-R/EpoR crosstalk

The system shown in Supplementary Figure 1 was modeled using coupled ordinary differential equations. All reactions were modeled using mass action kinetics. Furthermore, the pathway was modeled as an open system where new receptors and STAT proteins were synthesized at a constant rate. Activated receptors were lost due to receptor internalization and degradation [32]. Activated STAT degraded by the proteasome was modeled as slow first order decay [33]. The details of receptor-mediated ROS production were collapsed into a single reaction whereby active receptor produced ROS which could oxidize PTP1B or get degraded by cellular ROS scavengers (not modeled explicitly). Oxidized PTP1B could be reduced back to its active form.

The parameters of the model were estimated by fitting only the Epo signaling module (along with ROS production and PTP1B oxidation) to STAT5 phosphorylation data previously published for the Epo signaling pathway [34]. The same parameter estimates were used for the IL4 module since we expect similar qualitative behavior in both modules. The species used in the model and their initial values are listed in Supplementary Table 1. The reaction rate parameters used and the differential equation system are tabulated in Supplementary Tables 2 and 3, respectively. The qualitative fits are shown in Supplementary Figure 2.

The model without crosstalk (Fig. 3B) was obtained by setting the rate constants k8, k11, k13 and k14 to 0 (see Supplementary Figure 1).

The hypothetical redundant model (Fig. 3C) was obtained by making the following modifications to the model in Supplementary Figure 1: i) JAK3* catalyzed phosphorylation of STAT5 was added (k = 0.8); ii) JAK2* catalyzed phosphorylation of STAT6 was added (k = 0.8); iii) PTP1B catalyzed dephosphorylation of STAT6 was turned off; iv) PTP1B catalyzed dephosphorylation of JAK3* was added (k = 1.2).

Appendix F. Simulation of adaptive evolution

An n-node network with directed edges, representing a hypothetical signaling pathway, was constructed (n values 3 and above were used). The nodes represent molecules and the edges represent reactions. The pathway was modeled using mass action kinetics using either linear or non-linear reaction rates. Two nodes were arbitrarily chosen as input and output nodes and degeneracy was calculated using the Lyapunov method described above (Appendix C). To model evolution, an initial population of genetically similar individuals was created. Genetic similarity here is represented by similar strengths of connections, or reaction rates, between the nodes. Further, the initial population was constructed to meet certain constraints, such as all individuals could be required to have a positive degeneracy. This allowed comparison between different initial configurations. To simulate adaptation to a new environment, the initial population was made to evolve towards a new steady state value. Evolution followed cycles of mutations and selections. Mutations were modeled by making small random changes to the reaction rates. Selection was based on fitness of individuals in the new environment; fitter individuals had a higher chance of passing their traits to the next generation. Population averages of the steady state value of the output node and degeneracy were monitored over evolution.

Footnotes

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

References

1. Hartwell L, Hopfield J, Leibler S, Murray A, et al. From molecular to modular cell biology. Nature. 1999;402(6761):47. [PubMed]
2. Tononi G, Sporns O, Edelman G. A measure for brain complexity: relating functional segregation and integration in the nervous system. Proceedings of the National Academy of Sciences. 1994;91(11):5033. [PubMed]
3. Tononi G, Sporns O, Edelman G. Measures of degeneracy and redundancy in biological networks. Proceedings of the National Academy of Sciences of the United States of America. 1999;96(6):3257. [PubMed]
4. Edelman G, Gally J. Degeneracy and complexity in biological systems. Proceedings of the National Academy of Sciences. 2001;98(24):13763. [PubMed]
5. Stelling J, Sauer U, Szallasi Z, Doyle F, III, Doyle J. Robustness of cellular functions. Cell. 2004;118(6):675–685. [PubMed]
6. Rizk A, Batt G, Fages F, Soliman S. A general computational method for robustness analysis with applications to synthetic gene networks. Bioinformatics. 2009;25(12):i169. [PMC free article] [PubMed]
7. Kitano H. Towards a theory of biological robustness. Molecular systems biology. 3(1) [PMC free article] [PubMed]
8. Bogachev V, Krylov N, Röckner M. Elliptic and parabolic equations for measures. Russian Mathematical Surveys. 2009;64:973.
9. Huang W, Liu Z, Ji M, Yi Y. Fokker-Planck Equations. 2011. Stochastic Stability of Invariant Measures: Part I. Preprint.
10. Sun Han T. Multiple mutual informations and multiple interactions in frequency data. Information and Control. 1980;46(1):26–45.
11. Yeung R. A first course in information theory. Vol. 1. Plenum Pub Corp; 2002.
12. Pesin Y. Dimension theory in dynamical systems: contemporary views and applications. University of Chicago Press; 1997.
13. Li Y, Yi Y. Random perturbation in the vicinity of attractor. preprint.
14. Ludwig D. Persistence of dynamical systems under random perturbations. Siam Review. 1975:605–640.
15. Day M, Darden T. Some regularity results on the ventcel-freidlin quasi-potential function. Applied Mathematics and Optimization. 1985;13(1):259–282.
16. Day M. Regularity of boundary quasi-potentials for planar systems. Applied mathematics & optimization. 1994;30(1):79–101.
17. Freĭdlin M, Wentzell A. Random perturbations of dynamical systems. Vol. 260. Springer Verlag; 1998.
18. Xiao D, Li W. Limit cycles for the competitive three dimensional lotka-volterra system. Journal of Differential Equations. 2000;164(1):1–15.
19. Shuai K, Liu B. Regulation of jak–stat signalling in the immune system. Nature Reviews Immunology. 2003;3(11):900–911. [PubMed]
20. Sharma P, Chakraborty R, Wang L, Min B, Tremblay M, Kawahara T, Lambeth J, Haque S. Redox regulation of interleukin-4 signaling. Immunity. 2008;29(4):551–564. [PMC free article] [PubMed]
21. Kelly-Welch A, Hanson E, Boothby M, Keegan A. Interleukin-4 and interleukin-13 signaling connections maps. Science. 2003;300(5625):1527. [PubMed]
22. Constantinescu S, Ghaffari S, Lodish H. The erythropoietin receptor: structure, activation and intracellular signal transduction. Trends in endocrinology and metabolism. 1999;10(1):18–23. [PubMed]
23. Lu X, Malumbres R, Shields B, Jiang X, Sarosiek K, Natkunam Y, Tiganis T, Lossos I. Ptp1b is a negative regulator of interleukin 4–induced stat6 signaling. Blood. 2008;112(10):4098. [PubMed]
24. Myers M, Andersen J, Cheng A, Tremblay M, Horvath C, Parisien J, Salmeen A, Barford D, Tonks N. Tyk2 and jak2 are substrates of protein-tyrosine phosphatase 1b. Journal of Biological Chemistry. 2001;276(51):47771. [PubMed]
25. Röchman Y, Spolski R, Leonard W. New insights into the regulation of t cells by γc family cytokines. Nature Reviews Immunology. 2009;9(7):480–490. [PMC free article] [PubMed]
26. Chen W, Schoeberl B, Jasper P, Niepel M, Nielsen U, Lauffenburger D, Sorger P. Input–output behavior of erbb signaling pathways as revealed by a mass action model trained against dynamic data. Molecular systems biology. 5(1) [PMC free article] [PubMed]
27. Oda K, Matsuoka Y, Funahashi A, Kitano H. A comprehensive pathway map of epidermal growth factor receptor signaling. Molecular systems biology. 1(1) [PMC free article] [PubMed]
28. van Wageningen S, Kemmeren P, Lijnzaad P, Margaritis T, Benschop J, de Castro I, van Leenen D, Groot Koerkamp M, Ko C, Miles A, et al. Functional overlap and regulatory links shape genetic interactions between signaling pathways. Cell. 2010;143(6):991–1004. [PMC free article] [PubMed]
29. Hennessy B, Smith D, Ram P, Lu Y, Mills G. Exploiting the pi3k/akt pathway for cancer drug discovery. Nature Reviews Drug Discovery. 2005;4(12):988–1004. [PubMed]
30. Lander A. Pattern, growth, and control. Cell. 2011;144(6):955–969. [PMC free article] [PubMed]
31. Doyle J, Csete M. Rules of engagement. Nature. 2007;446(7138):860–860. [PubMed]
32. Becker V, Schilling M, Bachmann J, Baumann U, Raue A, Maiwald T, Timmer J, Klingmüller U. Covering a broad dynamic range: information processing at the erythropoietin receptor. Science. 2010;328(5984):1404. [PubMed]
33. Wang D, Moriggl R, Stravopodis D, Carpino N, Marine J, Teglund S, Feng J, Ihle J. A small amphipathic α-helical region is required for transcriptional activities and proteasome-dependent turnover of the tyrosine-phosphorylated stat5. The EMBO Journal. 2000;19(3):392–399. [PubMed]
34. Swameye I, Müller T, Timmer J, Sandra O, Klingmüller U. Identification of nucleocytoplasmic cycling as a remote sensor in cellular signaling by databased modeling. Proceedings of the National Academy of Sciences of the United States of America. 2003;100(3):1028. [PubMed]