16752

J. Phys. Chem. B 2008, 112, 16752–16758

Predicting Perfect Adaptation Motifs in Reaction Kinetic Networks
Tormod Drengstig
Department of Electrical Engineering and Computer Science, UniVersity of StaVanger, StaVanger, Norway

Hiroki R. Ueda
Center for DeVelopmental Biology, RIKEN, Kobe, Japan

Peter Ruoff*
Centre for Organelle Research, Faculty of Science and Technology, UniVersity of StaVanger, N-4036
StaVanger, Norway
ReceiVed: July 31, 2008; ReVised Manuscript ReceiVed: October 02, 2008

Adaptation and compensation mechanisms are important to keep organisms fit in a changing environment.
“Perfect adaptation” describes an organism’s response to an external stepwise perturbation by resetting someof its variables precisely to their original preperturbation values. Examples of perfect adaptation are found inbacterial chemotaxis, photoreceptor responses, or MAP kinase activities. Two concepts have evolved for
how perfect adaptation may be understood. In one approach, so-called “robust perfect adaptation”, the adaptation
is a network property (due to integral feedback control), which is independent of rate constant values. In theother approach, which we have termed “nonrobust perfect adaptation”, a fine-tuning of rate constant values
is needed to show perfect adaptation. Although integral feedback describes robust perfect adaptation in generalterms, it does not directly show where in a network perfect adaptation may be observed. Using control theoretic
methods, we are able to predict robust perfect adaptation sites within reaction kinetic networks and show thata prerequisite for robust perfect adaptation is that the network is open and irreversible. We applied the method
on various reaction schemes and found that new (robust) perfect adaptation motifs emerge when consideringsuggested models of bacterial and eukaryotic chemotaxis.
Introduction
In order to cope with a changing environment, organisms
have the ability to adapt1,2 to persistent external perturbations.
Several modes of adaptation exist,3,4 ranging from noadaptation at all to partial adaptation, perfect adaptation, or
overadaptation (Figure 1a). The term perfect adaptation
describes a system’s response when certain variables returnto precisely their original preperturbation values during a
stepwise perturbation.
Perfect adaptation has been found in a variety of instances,for example, in bacterial5-11 and eukaryotic12 chemotaxis, incertain photoreceptor responses,13,14 and in MAP-kinase
regulation.15-17 Two concepts have evolved for how perfectadaptation may be understood. In one approach, a fine-tuning/balancing of the rate constants is needed to show perfectadaptation,8,18-21 while in the other approach, the adaptation process is a network property which needs nofine-tuning.12,16,22-24 Yi et al.23 coined the term for thisbehavior as “robust perfect adaptation” and indicated that
this form of adaptation can be formulated in terms of integral
feedback,16,23 a concept employed in control theory.25-29
Although integral feedback can describe robust perfectadaptation in general terms, it does not directly predict where
in a network it will occur.
We became interested in the question of how to predictthe location of robust perfect adaptation sites in reaction* To whom correspondence should be addressed. Tel.: (47) 5183-1887.
Fax: (47) 5183-1750.

kinetic networks. By transforming an earlier global approachfor the condition of temperature compensation30,31 to a single
reaction site within a network receiving a signal from a
receptor, it became evident that a necessary (but not
sufficient) condition for robust perfect adaptation is related
to the presence of a zero control (sensitivity) coefficient32-35
of the actual steady-state concentration with respect to therate constant at this site.36
In this paper, we combine kinetic and control engineering
methods to predict where in a kinetic network perfect adaptation
motifs can occur and how the adaptation kinetics, either robust
or nonrobust, can be characterized by their transfer functions
and frequency responses.
Computational Methods
MATLAB (www.mathworks.com) was used to determine
transfer functions as well as step and frequency responses onlinearized models. For comparison with the actual nonlinear model,
rate equations were solved numerically by using the FORTRAN
subroutine LSODE (Livermore Solver of Ordinary Differential
Equations)37 and SIMULINK (www.mathworks.com).
Results
Transfer Function and Adaptation Kinetics. We considera reaction kinetic network with M chemical components (Ij) andN elementary reactions steps, where each step i is associatedwith a rate constant ki. The network can be stimulated bychanging one of the rate constants (ki) by means of a step

10.1021/jp806818c CCC: $40.75  2008 American Chemical Society
Published on Web 12/02/2008

Predicting Perfect Adaptation Motifs

J. Phys. Chem. B, Vol. 112, No. 51, 2008 16753
state concentration Ij,ss and rate constant ki are connected to eachother. Such a description is given by the transfer function matrix
H(s).25-27,29 In order to find H(s), eq 1 is first linearized aroundthe steady-state (ss) values Iss ) [I1,ss,..., IM,ss] and the preperturbation values of k ) [k1,..., kN], giving the following linear
model∆İ(t) ) A · ∆I(t) + B · ∆k(t)(2)

where ∆k(t) ) [∆k1(t),..., ∆kN(t)]T is a vector of small deviations
around the preperturbation values of k, and ∆I(t) ) [∆I1(t),...,∆IM(t)]T is a vector of small deviations around the steady-state
concentrations Iss. The matrix A is the system matrix (M × M)
defined by Aij ) (∂fi/∂Ij)|ss, and B is the input matrix (M × N)
defined by Bij ) (∂fi/∂kj)|ss. Applying the Laplace transform (fromthe real-valued t axis to the complex-valued s plane) to thelinearized model in eq 2 gives the M × N transfer function
matrix H(s) as25

H(s) )

Figure 1. (a) Step responses, (b) zeros, and (c) frequency responses
for four different transfer functions showing (1) no adaptation, (2) partial
adaptation, (3) perfect adaptation, and (4) overadaptation. The fourtransfer functions are described by Hi(s) ) ni(s)/d(s), i ) 1,..., 4, withthe same denominator d(s) ) (0.2s + 1)(0.35s + 1)(0.45s + 1), that
is, a third-order overdamped system. In order to visualize the impact
of different zeros, the numerators are n1(s) ) 1.5, n2(s) ) (s + 1),
n3(s) ) s, and n4(s) ) s - 1.

function. This stimulation is considered to occur due to a signal
from a receptor acting specifically and proportionally (at least
at low stimulation amplitudes) on ki due to an external stepwise
perturbation of the receptor. In this respect, the rate constants
are considered to be dependent upon time t. The kinetics of the
network are described by the rate equations for each chemicalcomponent Ij

dIj(t)) fj(k1(t), ..., kN(t), I1(t), ..., IM(t))dt(1)

In order to find the steady-state concentrations Ij,ss, the derivatives in eq 1 are set to zero, and the resulting equations are
solved with respect to the steady-state concentrations.
Perfect adaptation sites in a network are identified by applying
a stepwise perturbation on each rate constant ki and finding thecorresponding perfectly adapted responses for some of theintermediates Ij. A necessary (but not sufficient) condition forrobust perfect adaptation is that the control coefficient32,38 C Ij,ss)ki∂ ln Ij,ss/∂ ln ki becomes zero. This condition indicates that thesteady-state concentration Ij,ss is independent of (changes in)
ki. C Ij,ss ) 0 is only a necessary condition because there may bekisteady-state
concentrations which are completely independentof certain rate constants. For example, the change of a rate
constant in an irreversible chain of consecutive reactions will
not affect any steady-state level located upstream in the chain.
We call such a relationship between rate constants ki and steadystate levels Ij,ss disconnected. Therefore, an additional functional
relationship is needed, which describes if and how the steady-∆I(s)) (sI - A)-1B
∆k(s)

(3)

H(s) describes the connections between a small change in all
possible inputs, here the array of Laplace transformed rate
constants ∆k(s) ) [∆k1(s),..., ∆kN(s)]T, and the resulting smallchanges in all possible outputs, that is, the array of all Laplacetransformed concentration changes ∆I(s) ) [∆I1(s),..., ∆IM(s)]T.
I is the identity matrix (M × M). The matrix (sI - A)-1 iscalculated by means of the cofactors of (sI - A) divided by thedeterminant of (sI - A). The condition det(sI - A) ) 0
determines the poles of the system, which are identical to theeigenvalues of (sI - A). In general, the elements Hji(s) of thetransfer function matrix H(s) are written asnHji(s) )∆Ij(s))∆ki(s)K·∏ (s - zr)r)1m∏ (s - pq)

)

n(s)
d(s)

(4)

q)1

where zr are the transfer function’s zeros, pq are the poles, andK is the gain. As indicated in eq 4, Hji(s) is described as theratio between two polynomials, the numerator polynomial n(s)and the denominator polynomial d(s).
Robust and Nonrobust Perfect Adaptation. When for a given
Hji(s) n(s) has a zero in origo of the s plane regardless of thevalues in any of the rate constants, the system shows robust
perfect adaptation in Ij(t) with respect to a stepwise increase in
the rate constant ki (Figure 1a). The reason for this is that thesolution s ) 0 for n(s) ) 0, as shown in case 3 in Figure 1b, israte-constant-independent. Hence, a variation in ki will not affectthe steady-state concentration of Ij(t). In situations when thesolution to n(s) ) 0 has a zero in origo of the s plane but thissolution is dependent on the values of certain rate constants,
the system shows nonrobust (rate-constant-dependent) perfectadaptation in Ij(t) with respect to a stepwise change in ki.
Partial Adaptation and OWeradaptation. Partial adaptation(Figure 1a) is observed when the solution to n(s) ) 0 is real
and lies in the left half of the s plane (Figure 1b). Overadaptation
is observed when the solution to n(s) ) 0 is real and lies in theright half of the s plane (Figure 1b). In both cases, the solution
to n(s) ) 0 will be rate-constant-dependent.

16754 J. Phys. Chem. B, Vol. 112, No. 51, 2008

Drengstig et al.
be open, that is, capable of maintaining a steady state. Closed
systems, which are confined by the Principle of Detailed
Balance,44 approach equilibrium and will generally not exhibit,
as illustrated below, robust perfect adaptation. However, certain
closed systems which contain positive and/or negative feedback
loops can transiently exhibit oscillations and steady states55 andmay therefore, during such a transient, show robust perfectadaptation. The second condition to be met in order to get robust
perfect adaptation is that the network contains one or more
irreversible reaction steps. The placement of these irreversible
reactions within the network is crucial for getting robust perfectadaptation (see below).
To illustrate this, we consider the following reversible twocomponent system of first-order reactions, where the wavy line
indicates a (receptor-mediated) stepwise stimulation in k2

The corresponding transfer function HB,k2(s) is
HB,k2(s))
Figure 2. (a) Step responses, (b) zeros, and (c) frequency responses
for three different transfer functions, all showing perfect adaptation,that is, they all have at least a zero in origo. All three transfer functions
are given as Hi(s) ) ni(s)/d(s), i ) 1,..., 3, with the same denominatoras that in Figure 1. The basic transfer function has only a zero in origo,
that is, n2(s) ) s. In order to visualize the impact of an additional zeroin the left and right half-planes, respectively, the numerators of thetwo other transfer functions are n1(s) ) (s + 1)s and n3(s) ) (1 - s)s.

No Adaptation. If n(s) is a constant and does not have any
zeros, no adaptation can be observed (Figures 1a,b).
Figure 1 summarizes the various adaptation types with therespective zeros for n(s) and the corresponding frequency
responses shown as Bode plots (Figure 1c, only magnitudes areshown). In the case Hji(s) ) 0, Ij,ss is (what we called)“disconnected” from the rate constant ki, and consequently, C Ij,sski) 0.
Influence of the Transfer Function’s Poles and Zeros on
Adaptation Kinetics. The solution to d(s) ) 0 identifies the polesof the transfer function. If the poles are negative and real, thesystem is overdamped, and the step response has no overshooting. If the poles are complex conjugated with a negative realpart, the system is underdamped, and the step response shows
overshooting. However, the total response of a system will be
determined by both the poles and the zeros of the transfer
function. Figure 2 shows how the step response of a systemwith a zero in origo (midfigure) is influenced by an additional
zero in the left half-plane (left figure) and by an additional zeroin the right half-plane (right figure), leading to different perfectadaptation responses.
In the following, we will illustrate the above conditions fordifferent “perfect adaptation motifs” and show how changes inthese motifs can alter the adaptation response.
Conditions for Robust Perfect Adaptation. To get robust
perfect adaptation in reaction kinetic networks, two conditions
are necessary. The first condition is that the network needs to

)

∆B(s)
∆k2(s)Ass(s + k-1)

s + (k3 + k-1 + k-2 + k2)s + k3(k2 + k-1) + k-1k-2
2(5)
Since the solution to n(s) ) 0 for HB,k2(s) is given by s ) -k-1,
a step in k2 will show partial adaptation in B(t) as schematicallyindicated in case 2 of Figure 1b. When k-1 ) 0, HB,k2(s) has azero in origo, and a stepwise change in k2 will now lead torobust perfect adaptation in B(t). A more detailed descriptionof the steady states of this system together with a derivation ofeq 5 is given in the Supporting Information.
Thus, a simple motif that can show robust perfect adaptationis scheme M1 where the irreversibility in the first step is essential

for getting robust perfect adaptation in B(t) through a stepwise
stimulation of k2 (or k-2). To get a rapid overview where, inthe network, robust perfect adaptation occurs, circular tags areplaced at the rate constants (inputs) and at the associatedperfectly adapting intermediates (outputs), as indicated inscheme M1.
In the following, we analyze the adaptation response for avariety of irreversible (nonequilibrium) steady-state networks,
where some of them are related to adaptation dynamics observedin biological systems.
Network Size and Frequency Response. We consider thefollowing chain of consecutive reactions. Circular tags indicate

the rate constants (inputs) and the associated intermediates(outputs) which show robust perfect adaptation. For example,

Predicting Perfect Adaptation Motifs

J. Phys. Chem. B, Vol. 112, No. 51, 2008 16755

robust perfect adaptation is observed in every Bi when k2
undergoes a stepwise change. The transfer function from the
input ∆k2(s) to the various ∆Bi’s is given by∆Bl(s)
HBl,k2(s) ))∆k2(s){

Ass · sl+2

∏ (s + kj)

if l ) 1j)2

l+1

Ass · s

∏ ki
i)3

l+2

∏ (s + kj)(6)
if 2 e l e nj)2In both cases of eq 6, the only solution to the numerator
polynomial n(s) ) 0 is s ) 0, that is, a zero in origo independentof the rate constants.
Besides studying the step response of a system, we cananalyze how a system responds to external periodic stimulations.25,39
This is generally done with a sinusoidal perturbation with avarying frequency ω, which is applied to the linearized model.
The amplitudes and phases are determined as a function of ωand presented in form of a Bode plot.25,27
Figure 3 shows the step responses and Bode plots for M2
with n ) 6. Although each of the Bi’s shows perfect adaptation(Figure 3a) with similar adaptation kinetics, the frequency
responses for the Bi’s differ when a sinusoidal perturbation isapplied to k2. With increasing n in scheme M2 the negative
slope at high frequencies ω (Figure 3b) will decrease by anadditional 20 dB per decade change in ω for each Bi that is
added to the network. This means that for a given experimental
result where the step response for an intermediate looks likeone of the responses in Figure 3a, the size (or the order) of the
network cannot be determined unless a frequency response as
shown in Figure 3b,c is obtained.
Application of ConsecutiWe Reaction Models. Lumping allintermediates Bi of scheme M2 into a single variable B gives
motif M1*, which has also been studied by Hao et al.16 The

step response of B(t) is shown in Figure 4a, and this response
is strikingly similar to the photoadaptation in Limulus as

Figure 4. (a, upper panel) Robust perfect adaptation in B(t) of scheme
M1* with rate equations dA/dt ) k1 - k2 · A and dB/dt ) k2 · A - k3 · B.
The time and concentrations are given in arbitrary units (a.u.). Rate constantvalues are k1 ) 1.0, k2 ) 0.55, and k3 ) 0.5. At time t ) 30, (indicated by1), k2 is increased from 0.55 to 1.0 and then reset to 0.55 at t ) 60 (at 2).
Hatched areas indicate the amount of released/absorbed B during adaptation.
The amount of released/absorbed B during adaptation is given by (f 1)k1/fk2k3, where f describes the stimulation strength, that is, the ratio
perbetween the perturbed k2 (kper
2 ) and its preperturbed value (k2) with k2 )f · k2 (see Supporting Information). This indicates that for high stimulations(f . 1), the amount of released/absorbed B becomes independent of thestimulation strength f. (a, lower panel) Same system as that in the upper
panel, but at t ) 30 (at 1), k2 is decreased from 0.55 to 0.316. At 2 (t )60), k2 is reset to 0.55. (b) Photoadaptation with respect to a stepwise
perturbation of a single photoreceptor in Limulus. Shown are changes in
the optic nerve discharge frequency for a 2 s incremental flash (upper panel)or a 2 s decremental flash (lower panel). Open and closed circles represent
experimental determinations of the nerve’s discharge frequencies. The
bottom lines in the upper and lower panels indicate the onset andtermination of the stimulus. Redrawn after Figure 16 from Ratliff et al.13(c) Partial adaptation in B and B′ of two coupled M1* motifs (schemeM1**) showing close resemblance to experimental observations (Figure4d). The following rate equations and rate constants have been used: dA/dt ) k1/(KI1 + A′) - k2 · A; dB/dt ) k2 · A - k3 · B/(KI2 + A′); dA′/dt )k4/(KI1 + A) - k5 · A′; and dB′/dt ) k5 · A′ - k6 · B′/(KI2 + A), with k1 ) k4
) 1.0, k2 ) k5 ) 0.55, k3 ) k6 ) 0.5, KI1 ) 0.5, and KI2 ) 0.35. At t )50, k2 is increased from 0.55 to 1.0 and reset to 0.55 at t ) 100. (d) Partial
adaptation in two adjacent photoreceptors of Limulus (redrawn after
Figure 17 from Ratliff et al.13). One receptor was constantly illuminated(black dots), while the other receptor (open circles) received an increased
illumination between 0 and 2 s and was reset to its preperturbation value
after 2 s.

determined by measuring the activity (discharge frequency) ofa single optic nerve13 shown in Figure 4b.
Remarkably, when two robust perfect adaptation motifs M1*
are coupled together by mutual inhibition as shown in scheme
M1** we observe adaptation responses which are not only
Figure 3. (a) Step responses for the first six Bi’s (identified by theiri’s) in scheme M2 when a step is applied on k2. Rate equations are
dA/dt ) k1 - k2 · A, dB1/dt ) k2 · A - k3 · B1, and dBi/dt ) ki+1 · Bi ki+2 · Bi, where all rate constants are, for the sake of simplicity, set equal
to 1.0. Panels (b) and (c) show the corresponding Bode plot, that is,
the influence of frequency on the magnitude and phase of the six Bi’s,
respectively.

16756 J. Phys. Chem. B, Vol. 112, No. 51, 2008

Drengstig et al.

radically different from the isolated motif’s (M1*) responses
but which show a striking similarity to experimental observations when two photoreceptors of Limulus are mutually inhibiting each other13 (Figure 4c,d)!
This suggests that schemes M1* or M1**, despite theirsimplicity, could be considered as building blocks for thephotoadaptation responses of single and coupled optic nervesin Limulus.
Moreover, Brodie et al.40 studied the frequency response of
the Limulus lateral eye. Assuming a reaction model as scheme
M2, the experimental Bode plot in the Brodie et al. study(Figure 4 in ref 40) indicates a network length of n ≈ 16
(including models of experimental setup) downstream from thesite at which the signal from the receptor is perceived.
Influence of Feedback and Feedforward Loops. Feedback
and feedforward loops are common regulatory elements. When
including negative feedback or positive feedforward loops inconsecutive networks, we generally observe additional over- or
undershooting behaviors, but the distribution of the robust
perfectly adapted sites within the network remains mostly
unaffected. A more detailed description of these behaviors isgiven in the Supporting Information.
Cyclic Reactions. Turning to cyclic reactions, we foundnetworks that can show both rate-constant-independent (robust)and rate-constant-dependent (nonrobust) perfect adaptation. For
the sake of simplicity, we consider here cyclic schemes withthree intermediates. A selection of four schemes (1-4) is given
in M3. Although reaction scheme 1 has an irreversible input(k1), the two output reactions with rate constants k3 and k5
produce a transfer function matrix H(s) where none of its
elements Hji(s) have a zero in the origo independent of the rate
constants.

Figure 5. Different adaptation kinetics observed for scheme 3 in M3.
(a) No adaptation in A for a step perturbation in k5, (b) partial adaptation
in C for a step perturbation in k-6, (c) perfect (robust) adaptation in B
for a step perturbation in k-4, and (d) overadaptation in A for a step
perturbation in k-4. Rate equations: dA/dt ) k1 - (k2 + k-6) · A + k6 · C+ k-2 · B, dB/dt ) k2 · A - (k4 + k-2 + k3) · B + k-4 · C, and dC/dt ) k5
- (k-4 + k6) · C + k4 · B + k-6 · A, with rate constants k1 ) 1.1, k2 )1.2, k-2 ) 1.7, k3 ) 1.3, k4 ) 1.4, k-4 ) 1.8, k5 ) 1.5, k6 ) 0.1, andk-6 ) 1.9.

0. Reversing one of the output reactions in scheme 1 (k5, see
scheme 3), B shows robust perfect adaptation with respect toall rate constants in the ring. Increasing the number ofirreversible steps of scheme 3 leads, as shown in 4, to robust
perfect adaptation in all intermediates for stepwise perturbations
in different ki. This illustrates, as already seen for consecutive
reaction M1, that the presence of irreversibility and theplacement of input and output reactions are crucial for theoccurrence of robust perfect adaptation. We found it interestingthat in scheme 3 (and in larger cyclic networks) all types ofadaptation responses can be observed (Figure 5).
Schemes 2 and 4 of M3 show a close resemblance to the
Barkai-Leibler (BL) model22 of bacterial chemotaxis. The BL
model has been discussed41 using the (thermodynamically)closed form as shown in M4

Therefore, no robust perfect adaptation can be obtained forthis scheme. However, the transfer function from ∆k4(s) to ∆A(s)HA,k4(s) )Bss((-k-2 + k6)s - k5k-2 + k6k3)
n(s)
(7)
)
d(s)
d(s)

has the following solution to n(s) ) 0s)

k5k-2 - k6k3
-k-2 + k6(8)which is rate-constant-dependent. Therefore, a zero in origo for
HA,k4(s) is obtained with perfect adaptation in A(t) when a step
perturbation is applied to k4 if k5k-2 - k6k3 ) 0. Otherwise,
partial adaptation or overadaptation (dependent on the valuesof the rate constants) can be observed. Thus, scheme 1 candemonstrate rate-constant-dependent nonrobust perfect adaptation.
By increasing the number of irreversible steps in 1, intermediates B and C become robust perfect adapted (see scheme 2).
No robust perfect adaptation in A is observed even when k-2 )The presence of a zero-order reaction step from C to Ri
assures perfect adaptation in the active form of the receptor Ra/
when applying a stepwise perturbation to k1 or k2. In view of
the results for schemes 2 and 4 of M3, the presence of a zeroorder term in M4 is not absolutely required when describing
scheme M4 as an open far-from-equilibrium steady-statesystem42,43 by explicitly including system input and system
output reactions. It may also be noted that by representing
scheme M4, as shown, formally as a closed system,41 it violates
the Principle of Detailed Balance44 and may, for certain rate
constant values, lead to negative concentrations in C.

Predicting Perfect Adaptation Motifs

J. Phys. Chem. B, Vol. 112, No. 51, 2008 16757

TABLE 1: Overview over kj f ki Substitutions (scheme M5)
index isubstituted by index j2
45
number of perfect adaptations

1
3
6
0

6
1
3
3

3
6
1
3

6
3
1
2

3
1
6
2

1
6
3
2

Scheme M5 shows a model suggested by Levchenko andIglesias12 in order to describe perfect adaptation in eukaryotic
chemotaxis. In this model, an activator A drives the receptor to
its active form Ra/, while an inhibitor I converts the active
receptor into an inactive form Ri.

Levchenko and Iglesias observed that perfect adaptation is
obtained when activation and inhibition of the receptor are
mediated by a common signal which activates A and I by aperfixed proportion, that is, kper5 ) Rk1 , where R is any real positive
constant. This suggestion for getting perfect adaptation isindicated in scheme M5 by white circular tags at k5 and k1 withcorresponding white tags at Ra/ and Ri (indicating that both Ra/
and Ri show perfect adaptation). However, by combining themethod of zero control coefficients with an analysis of thecorresponding transfer functions, five additional possibilities for
perfect adaptation can be found, complementing the suggestion
by Levchenko and Iglesias. Assuming first-order kinetics,12 the*can besteady-state concentration of active receptor Ra,ss
expressed in terms of the total receptor concentration Rtot )*+ Ri,ss as (see Supporting Information)
Ra,ss
/
Ra,ss
)

k1k3k6R
k1k3k6 + k2k4k5 tot(9)From eq 9, we see that for any substitution ki ) Rkj with i ∈/{2, 4, 5} and j ∈ {1, 3, 6}, the control coefficient C kRja,ss ) ∂ ln*
Ra,ss/∂ ln kj becomes zero, suggesting that a common stepwise
activation of ki and kj (kper) Rkperij ) may lead to perfectadaptation. Inspecting the corresponding transfer functions (seeSupporting Information), it was found that perfect adaptation
is observed in six cases (indicated in scheme M5) out of thenine possibilities that are suggested by a zero control coefficient/
C Rkja,ss (data not shown). Table 1 shows all possible rate constant
substitution combinations and the resulting number of perfectadaptation possibilities. In Table 1, substitutions which lead toperfect adaptation are indicated in boldface, while those showing
no perfect adaptation are indicated in italic.
Discussion
Adaptation processes are important for providing robustness45,46
to biological systems. To our knowledge, no systematic attempts
have so far been made to investigate the possible adaptation
responses (Figure 1) in reaction kinetic networks. Here, we haveshown how robust perfect adaptation sites and their responses(Figure 2) can be identified by finding zero control coefficientsand by analyzing the corresponding transfer functions. Somewhat to our surprise, already small networks can show arelatively large number of robust perfect adaptation sites, which,during evolution, may have been of importance where in a
network the signal from a receptor was received. For certain
cyclic networks such as scheme 3 of M3, we can even observe
simultaneously all adaptation types (Figure 1) in the same
network (Figure 5).
A prerequisite to observe robust perfect adaptation is thepresence of irreversibility in the network (as seen in thetransition from scheme S1 to M1) and treating the network as
an open far-from-equilibrium system.43,42
While the identification of robust perfect adaptation sites canbe linked to the identification of (exact) zero control coefficientsand the analysis of the associated transfer functions, thedetermination of nonrobust (rate-constant-dependent) adaptationsites in a reaction network is less straightforward. It requires
an analysis on how control coefficients depend on the rate
constants. In addition, rate-constant-dependent adaptation may
even occur in the presence of nonzero control coefficients whenreceptor signals from the same environmental variable are acting
at different sites in the network.
To illustrate this latter case (which can be considered as an) Rkpersuggestion by Levchenko andextension of the kperij
12
Iglesias), consider a network where component Ij,ss is dependent on a subset of P rate constants. Each rate constant {ki}i)1,...,P
is associated with a receptor signal from the same environmental
variable Θ (light, a hormone, etc.) changing the values of the
ki’s through Θ. The influence of Θ on Ij,ss can be written asP
d ln Ij,ss∂ ln kiCkIj,ss)id ln Θ
∂ ln Θ
i)1

∑

(10)

When C kIj,ss* 0, perfect adaptation may still be possible, butipositive and negative contributions of the C Ikj,ss’s are needed suchithat eq 10 becomes zero. This type of balancing approachrequires a fine-tuning of the parameters ∂ ln ki/∂ ln Θ for agiven set of control coefficients. Equation 10 has been suggested
as a mechanism for temperature compensation47,48 of circadianrhythms30,47,49-52 and steady-state fluxes31 in ectotherms. An
analysis of the influence of temperature on reaction kinetic
networks in terms of control coefficients and transfer functions
will be given elsewhere.
The analysis presented here is based on linearizing a system
around a steady state. This has the advantage of making theanalysis more tractable.26 Although a linearized model is a good
approximation for a nonlinear system when applying small
perturbations,27,29,53,54 the linearized system may not describe
correctly system responses when using large perturbations. In
this respect, studies using linearized models may be complemented by numerical techniques54,56,57 or by trying to apply
methods such as feedback control.29
Adaptation mechanisms contribute to the robustness of biologicalsystems, and there are attempts to provide a mathematical description of robustness.45,46 In this respect, we think that the heredescribed methods are useful in identifying and characterizing
adaptation responses within reaction kinetic networks.
Acknowledgment. We thank Morten Tengesdal and Ludger
Rensing for comments on the manuscript. This research was
supported, in part, by the National Science Foundation underGrant No. PHY05-51164. The KITP preprint number of thismanuscript is NSF-KITP-07-213.

16758 J. Phys. Chem. B, Vol. 112, No. 51, 2008
Supporting Information Available: Derivations of thetransfer functions for scheme S1 and scheme M5, derivation of
the amount of released/absorbed B during adaptation in scheme
M1*, and a description of the influences of negative feedback
and positive feedforward loops on adaptation behaviors inconsecutive reactions. This material is available free of charge
via the Internet at http://pubs.acs.org.
References and Notes(1) Hochachka, P. W.; Somero, G. N. Biochemical Adaptation.
Mechanism and Process in Physiological EVolution; Oxford University
Press: Oxford, U.K., 2002.
(2) Holland, J. H. Adaptation in Natural and Artificial Systems; MIT
Press: Cambridge, MA, 1992.
(3) Asthagiri, A. R.; Lauffenburger, D. A. Annu. ReV. Biomed. Eng.
2000, 2, 31–53.
(4) Koshland, D. E., Jr.; Goldbeter, A.; Stock, J. B. Science 1982, 217(4556), 220–5.
(5) Berg, H. C.; Tedesco, P. M. Proc. Natl. Acad. Sci. U.S.A. 1975,
72 (8), 3235–9.
(6) Alon, U.; Surette, M. G.; Barkai, N.; Leibler, S. Nature 1999, 397(6715), 168–71.
(7) Bray, D. Proc. Natl. Acad. Sci. U.S.A. 2002, 99 (1), 7–9.
(8) Mello, B. A.; Tu, Y. Biophys. J. 2003, 84 (5), 2943–56.
(9) Berg, H. C. E. coli in Motion; Springer-Verlag: New York, 2004.
(10) Mello, B. A.; Tu, Y. Biophys. J. 2007, 92 (7), 2329–37.
(11) Hansen, C. H.; Endres, R. G.; Wingreen, N. S. PLoS Comput. Biol.
2008, 4, 0014-0027.
(12) Levchenko, A.; Iglesias, P. A. Biophys. J. 2002, 82 (1 Pt 1), 50–
63.
(13) Ratliff, F.; Hartline, H. K.; Miller, W. H. J. Opt. Soc. Am. 1963,
53, 110–20.
(14) He, Q.; Liu, Y. Genes DeV. 2005, 19 (23), 2888–99.
(15) Asthagiri, A. R.; Nelson, C. M.; Horwitz, A. F.; Lauffenburger,
D. A. J. Biol. Chem. 1999, 274 (38), 27119–27.
(16) Hao, N.; Behar, M.; Elston, T. C.; Dohlman, H. G. Oncogene 2007,26 (22), 3254–66.
(17) Mettetal, J. T.; Muzzey, D.; Gómez-Uribe, C.; van Oudenaarden,
A. Science 2008, 319, 482–4.
(18) Segel, L. A.; Goldbeter, A.; Devreotes, P. N.; Knox, B. E. J. Theor.
Biol. 1986, 120 (2), 151–79.
(19) Hauri, D. C.; Ross, J. Biophys. J. 1995, 68, 708–722.
(20) Spiro, P. A.; Parkinson, J. S.; Othmer, H. G. Proc. Natl. Acad. Sci.
U.S.A. 1997, 94 (14), 7263–8.
(21) Mello, B. A.; Tu, Y. Proc. Natl. Acad. Sci. U.S.A. 2005, 102 (48),17354–9.
(22) Barkai, N.; Leibler, S. Nature 1997, 387 (6636), 913–7.
(23) Yi, T. M.; Huang, Y.; Simon, M. I.; Doyle, J. Proc. Natl. Acad.
Sci. U.S.A. 2000, 97 (9), 4649–53.
(24) Csikasz-Nagy, A.; Soyer, O. S. J. R. Soc. Interface 2008, 5, S41–
S47.
(25) Wilkie, J.; Johnson, M.; Reza, K. Control Engineering. An
Introductory Course; Palgrave: New York, 2002.
(26) Milsum, J. H. Biological Control Systems Analysis; McGraw-Hill:
New York, 1966.
(27) Ingalls, B. P.; Yi, T.-M.; Iglesias, P. A. In System Modeling in
Cellular Biology; Szallasi, Z., Stelling, J., Periwal, V., Eds.; MIT Press:
Cambridge, MA, 2006.

Drengstig et al.
(28) El-Samad, H.; Goff, J. P.; Khammash, M. J. Theor. Biol. 2002,
214 (1), 17–29.
(29) Sontag, E. D. Mathematical Control Theory. Deterministic Finite
Dimensional Systems, 2nd ed.; Springer-Verlag: New York, 1998.
(30) Ruoff, P.; Loros, J. J.; Dunlap, J. C. Proc. Natl. Acad. Sci. U.S.A.
2005, 102 (49), 17681–6.
(31) Ruoff, P.; Zakhartsev, M.; Westerhoff, H. V. FEBS J. 2007, 274(4), 940–50.
(32) Savageau, M. A. Nature 1971, 229 (5286), 542–4.
(33) Savageau, M. A. Biochemical Systems Analysis. A Study of Function
and Design in Molecular Biology; Addison-Wesley: Reading, MA, 1976.
(34) Fell, D. Understanding the Control of Metabolism; Portland Press:
London and Miami, FL, 1997.
(35) Heinrich, R.; Schuster, S. The Regulation of Cellular Systems;
Chapman and Hall: New York, 1996.
(36) Ruoff, P. Perfect Adaptation with Zero Control Coefficients, KITP.
http://online.itp.ucsb.edu/online/bioclocks07/ruoff1/.
(37) Radhakrishnan, K.; Hindmarsh, A. C. Description and Use of
LSODE, the LiVermore SolVer for Ordinary Differential Equations. NASA
Reference Publication 1327, Lawrence Livermore National Laboratory
Report UCRL-ID-113855; National Aeronautics and Space Administration,
Lewis Research Center: Cleveland, OH, 1993.
(38) Burns, J. A.; Cornish-Bowden, A.; Groen, A. K.; Heinrich, R.;
Kacser, H.; Porteous, J. W.; Rapoport, S. M.; Rapoport, T. A.; Stucki, J. W.;
Tager, J. M.; Wanders, R. J. A.; Westerhoff, H. V. Trends Biochem. Sci.
1985, 19, 16.
(39) Ingalls, B. P. J. Phys. Chem. B 2004, 108, 1143–1152.
(40) Brodie, S. E.; Knight, B. W.; Ratliff, F. J. Gen. Physiol. 1978, 72(2), 167–202.
(41) Alon, U. An Introduction to Systems Biology: Design Principles
of Biological Circuits; Chapman & Hall: New York, 2006.
(42) Kondepudi, D.; Prigogine, I. Modern Thermodynamics. From HeatEngines to DissipatiVe Structures; John Wiley & Sons: Chichester, U.K.,
1998.
(43) von Bertalanffy, L. PerspectiVes on General System Theory; George
Braziller: New York, 1975.
(44) Moore, J. W.; Pearson, R. G. Kinetics and Mechanism, 3rd ed.;John Wiley & Sons: New York, 1981.
(45) Wagner, A. Robustness and EVolVability in LiVing Systems;
Princeton University Press: Princeton, NJ, 2005.
(46) Kitano, H. Mol. Syst. Biol. 2007, 3, 1–7.
(47) Hastings, J. W.; Sweeney, B. M. Proc. Natl. Acad. Sci. U.S.A. 1957,
43, 804–811.
(48) Hazel, J. R.; Prosser, C. L. Physiol. ReV. 1974, 54 (3), 620–677.
(49) Ruoff, P. J. Interdiscip. Cycle Res. 1992, 23, 92–99.
(50) Edwards, K. D.; Lynn, J. R.; Gyula, P.; Nagy, F.; Millar, A. J.
Genetics 2005, 170 (1), 387–400.
(51) Eckhardt, N. A. Plant Cell 2006, 18 (5), 1105–08.
(52) Gould, P. D.; Locke, J. C.; Larue, C.; Southern, M.; Davis, S. J.;
Hanano, S.; Moyle, R.; Milich, R.; Putterill, J.; Millar, A. J.; Hall, A. Plant
Cell 2006, 18 (5), 1177–87.
(53) Kholodenko, B. N.; Kiyatkin, A.; Bruggeman, F. J.; Sontag, E.;
Westerhoff, H. V.; Hoek, J. B. Proc. Natl. Acad. Sci. U.S.A. 2002, 99 (20),12841–6.
(54) Ross, J. J. Phys. Chem. A 2008, 112 (11), 2134–43.
(55) Ruoff, P.; Noyes, R. J. J. Phys. Chem. 1985, 89, 1339–1341.
(56) Ermentrout, B. Simulating, Analyzing, and Animating Dynamical
Systems. A Guide to XPPAUT for Researchers and Students; Society for
Industrial and Applied Mathematics: Philadelphia, PA, 2002.
(57) Szallasi, Z.; Stelling, J.; Periwal, V. System Modeling in Cellular
Biology; MIT Press: Cambridge, MA, 2006.

JP806818C

Supporting Information (10 pages)
Predicting perfect adaptation motifs in reaction kinetic
networks

Tormod Drengstig1 , Hiroki R. Ueda2 and Peter Ruoff3,∗
1
Department of Electrical Engineering and Computer Science, University of Stavanger,
N-4036 Stavanger, Norway; 2 Center for Developmental Biology, RIKEN, Kobe, Japan;
3
Centre for Organelle Research, University of Stavanger, N-4036 Stavanger, Norway

Derivation of the Transfer function HB,k2 (s) for motif S1, Eq. 5
Consider motif S1 (main paper):
k1
k-1

k2

A k B
-2

k3
k-3

(S1)

with the rate equations:dA(t)= Ȧ(t) = k1 (t) − k2 (t)A(t) + k−2 (t)B(t) − k−1 (t)A(t)dtdB(t)= Ḃ(t) = k2 (t)A(t) − k3 (t)B(t) − k−2 (t)B(t) + k−3 (t)
dt

(A1)
(A2)

and the rate constants k1 (t), k2 (t), k3 (t), k−1 (t), k−2 (t) and k−3 (t).
The steady state concentrations Ass and Bss are given by Eq. A3 and no perfect adaptation can exist, because both Ass and Bss depend on all rate constants such that for all
ki we have nonzero control coefficients, i.e., CkAiss 6= 0 and CkBiss 6= 0.
Ass =k−3 (k−1 + k2 ) + k1 k2
k−2 (k−3 + k1 ) + k1 k3; Bss =
k3 (k−1 + k2 ) + k−1 k−2
k3 (k−1 + k2 ) + k−1 k−2(A3)

However, by introducing an irreversible input to the system (as indicated in the main
paper), for example by setting k−1 = 0, robust perfect adaptation in B(t) with respect1

to a step-wise change in k2 (or k−2 ) can be observed. The reason for this is that B(t) is−3
now independent of k2 or k−2 , because Bss = k1 +k. However, B(t) is still connectedk3
to k2 and k−2 as can be seen by inspecting the transfer function HB,k2 (s) for the entire
scheme S1, Eq. 5 (see below).
Ass , however, still depends on the remaining rate constants and therefore shows noperfect adaptation. Introducing irreversibility in scheme S1 by setting k−2 = 0 or/and
k−3 = 0, will not lead to perfect adaptation for neither A nor B, because also in thesecases Ass and Bss still depend on all the other rate constants.
In order to find the transfer function matrix H(s), the system is linearized around thesteady states values Ass and Bss , and the rate constants k1 , k2 , k3 , k−1 , k−2 and k−3 ,which gives the following linear model∆Ȧ(t) = − (k2 + k−1 )∆A(t) + k−2 ∆B(t)+ ∆k1 (t) − Ass ∆k2 (t) − Ass ∆k−1 (t) + Bss ∆k−2 (t)(A4)∆Ḃ(t) =k2 ∆A(t) − (k3 + k−2 )∆B(t)+ Ass ∆k2 (t) − Bss ∆k3 (t) − Bss ∆k−2 (t) + ∆k−3 (t)(A5)

or, in matrix form

 


∆Ȧ(t)−(k2 + k−1 )
k−2
∆A(t)
=
k2
−(k3 + k−2 ) ∆B(t)∆Ḃ(t)

1 −Ass0−Ass Bss+
0 Ass −Bss
0
−Bss


∆k1 (t)  ∆k2 (t) 0  ∆k3 (t) 


∆k
(t)
1 
−1


∆k−2 (t)∆k−3 (t)(A6)

The transfer function matrix can be found from the relationship (Eq. 3, main paper):H(s) =∆I(s)= (sI − A)−1 B
∆k(s)

(A7)

where ∆k(s) = [∆k1 (s), ∆k2 (s), ∆k3 (s), ∆k−1 (s), ∆k−2 (s), ∆k−3 (s)]T , and∆I(s) = [∆A(s), ∆B(s)]T . Applying Eq. A7 to the linear model in Eq. A6 gives the2following transfer function matrix

−1 

1 −Ass0−Ass Bss 0
s + k2 + k−1−k−2
H(s) =
0 Ass −Bss0
−Bss 1−k2
s + k3 + k−2=

1·(s + k2 + k−1 )(s + k3 + k−2 ) − k2 k−2 s + k3 + k−2
k−2
1 −Ass0−Ass Bss 0·
k2
s + k2 + k−1
0 Ass −Bss0
−Bss 11·(s + k2 + k−1 )(s + k3 + k−2 ) − k2 k−2s + k3 + k−2 −Ass (s + k3 )−Bss k−2...
k2
Ass (s + k−1 ) −Bss (s + k2 + k−1 ) . . .

=. . . −Ass (s + k3 + k−2 ) Bss (s + k3 )k−2...
−Ass k2
Bss (s + k−1 ) s + k2 + k−1

from which the element HB,k2 (s) is found as
HB,k2 (s) =∆B(s)Ass (s + k−1 )=∆k2 (s)(s + k2 + k−1 )(s + k3 + k−2 ) − k2 k−2(A8)which is identical to Eq. 5 in the main paper.

Amount of Released/Absorbed B during Adaptation in Motif M1* (Fig. 4)
Setting Eqs. A1 and A2 (see above) to zero, with the additional condition that k−1 =k−2 = k−3 = 0, we get the steady state concentrations in A and B as Ass = k1 /k2 and
Bss = k1 /k3 , respectively. At t = 0 we assume that A(0) = Ass and B(0) = Bss , andthat k2 undergoes a step-wise change to f ·k2 with f > 0 and f 6= 1.
For t ≥ 0 the response kinetics of A and B are calculated as:Ass−f k2 tA(t) =
1 + (f − 1)ef
k1 (f − 1) −f k2 t
−k3 t
B(t) = Bss +e−ek3 − f k2

(A9)

(A10)

The amount of released or absorbed B during the robust perfect adaptation of B iscalculated by using the integralZ k1 (f − 1)  e−k3 t e−f k2 t −I(t) =
B(t) − Bss dt =+C
(A11)
k3 − f k2
k3
f k2
3

· kk2 k1 3 we get
Because limt→∞ I(t) = 0 and I(0) = − f −1
fZ ∞
f − 1 k1
B(t) − Bss dt = lim I(t) − I(0) =·t→∞
f k2 k3
0(A12)

Influence of Negative Feedback and Positive Feedforward on Adaptation
Feedback and feedforward loops are common regulatory elements. We illustrate herewith a few examples how negative feedback and positive feedforward can affect adaptation (a description of all possible combinations of negative/positive feedback/feedforward
loops is beyond the scope of this paper). First we consider an extension of motif M2
in the main paper by including a negative feedback from intermediate Bn to k1 :

KI
k1

A

k2

k

k

k

k

B1 3 B2 4 ... n+1 Bn n+2(M6)1
with the inhibition term (to be multiplied with k1 ) KI +B. KI is the inhibition constant.
n
The legends of Figs. S1a and b give the rate equations and numerical values for KI andrate constants. The transfer function from k2 to Bn (n ≥ 2), HBn ,k2 (s), is given by:∆Bn (s)HBn ,k2 (s) ==∆k2 (s)

n+1
Y

ki

ki +

n+2
Y

Ass ·si=31(KI +Bn,ss )2

n+1
Y
i=1

(A13)

(s + ki)i=2

From Eq. A13 we get s = 0 as the only solution to the numerator polynomial n(s) =0, i.e., a zero in origo independent of the rate constants. The transfer functions forintermediates B1 to Bn−1 show slightly different structures, but they all have a zero inorigo (data not shown). Hence, a step in k2 results in robust perfect adaption for each
Bi -intermediate. HBn ,k2 (s) has both real and complex-conjugated poles, which results
in an underdamped response (over- and undershooting) in the perfect adapted Bi ’s(Fig. S1a). For high KI values the response becomes overdamped (Fig. S1b), which isdue to the fact that the real pole dominates over the complex-conjugated poles as KI
increases. In case the negative feedback from Bn is acting downstream of k2 , at ki ’s
with i > 2, all Bi ’s show robust perfect adaptation. If the perturbing step is applied to
kn+2 , i.e., at the end of the network, none of the Bi ’s show perfect adaptation.

4

To illustrate the influence of positive feedforward, we consider two cases each shownin scheme M7.

Kact

Kact
k1

A

k2

k

B1 3 B 2

k4... Bn-1k Bnk

k1

n+2

n+1

A

k2

k

k

kn+2

B1 3 B2 4 ... Bn-1 n+1 Bnk(M7)
In the left scheme of M7 A acts positively on kn+1 , while in the right scheme B1 acts
positively on kn+1 . In each case the positive feedforward loop is realized by multiplying in the rate equations kn+1 with the factors Kact A or Kact B1 , respectively, where
Kact is an ”activation constant”. In the case A acts positively on kn+1 , concentration
Bn−1 is not perfectly adapted, because A does not show perfect adaptation. All other
Bi ’s show robust perfect adaptation. The transfer function from k2 to Bn , HBn ,k2 (s)for n ≥ 3, is given as:

HBn ,k2 (s) =∆Bn (s)=∆k2 (s)



n
n
Y
Y
−Kact ·Ass ·s·kn+1 Bn−1,ss (s + ki ) − Ass
ki
i=3

i=3

Y

n


(s + ki ) s + Kact ·Ass ·kn+1 s + kn+2
i=2

(A14)
HBn ,k2 (s) has several zeros. Based on the fact that there is a zero in origo (s = 0)
to the solution n(s) = 0, independent of the rate constants, Bn shows robust perfectadaptation. However, the total response in Bn will be influenced by the other zeros as
indicated in Fig. 2 of the main paper.
In case intermediate B1 acts positively on kn+1 (the right scheme in M7) all Bi ’s arerobust perfectly adapted. The transfer function from k2 to Bn , HBn ,k2 (s) for n ≥ 4, isgiven as:HBn ,k2 (s) =∆Bn (s)=∆k2 (s)



n
n
Y
Y
Kact ·Ass ·s·kn+1 Bn−1,ss ·s (s + ki ) + Bss
ki
i=4

i=3

Y

n


(s + ki ) s + Kact ·Bss ·kn+1 s + kn+2
i=2

Figs. S1c and d show response kinetics of the two cases in M7 for n = 4.

5

(A15)

0.5

0.08

a

negative feedback
KI = 00.4
1amplitude

amplitude

2
3
4

0.1

1

0.06

0.3
0.2

negative feedback
high KI

0.07

0.05

2
3

0.04
0.03

4

0.02

0

b

0.01
-0.1

0

-0.2
0

5

10

15

20

25

30

-0.01
0

35

5

10

0.5

1.2

c

1

3

0.3

amplitude

amplitude

0.4

positive feedforward
A acting on k5

0.8
0.6

15

time, a.u.

time, a.u.

0.4
2

0.2

d

positive feedforward
B1 acting on k5

1

2
4

0.1

0.2
1

0

0
4

-0.2
-0.4
0

2

-0.1

4

6

8

10

12

14

16

18

-0.2
0

20

3

5

10

15

20

25

time, a.u.

time, a.u.

Fig. S1. Perfect adaptation kinetics in B1 to B4 for negative feedback scheme M6 (n = 4). The
step perturbation is applied to k2 with (a) KI = 0 (strong inhibition) and (b) KI = 5.0 (weak1− k2 ·A, dB1 /dt = k2 ·A − k3 ·B1 , andinhibition). The rate equations are: dA/dt = KIk+B
4dBi /dt = ki+1 ·Bi−1 − ki+2 ·Bi , i = 2, . . . 4 with all rate constants equal to 1.0. (c) Adaptation
kinetics for left scheme in M7 (n = 4, and step-wise perturbation of k2 ) where A is acting
positively on k5 with rate equations: dA/dt = k1 − k2 ·A, dB1 /dt = k2 ·A − k3 ·B1 , dB2 /dt =k3 ·B1 − k4 ·B2 , dB3 /dt = k4 ·B2 − Kact ·A·k5 ·B3 , dB4 /dt = Kact ·A·k5 ·B3 − k6 ·B4 . All rate
constant values and Kact are set to 1.0. (d) Similar system as in (c) but B1 is acting positively
on k5 (right scheme in M7, n = 4). All rate constants and Kact are set to 1.0. Note that all Bi
show robust perfect adaptation, but B3 is the only intermediate which shows different adaptation kinetics (undershooting).

6Derivation of Eq. 9 for motif M5
From motif M5 the following relationship holds at steady state
k1
k2
k5
Iss =
k6

Ass =

(A16)
(A17)

Moreover, the following differential equation can be found for the dynamics of intermediate Ra∗ :dRa∗ (t)k5k1= −k4 Ra∗ (t) + k3 Ri (t)
(A18)
dt
k6
k2
Assuming steady state, i.e. dRdta (t) = 0 gives∗∗
k2 k4 k5 Ra,ss= k1 k3 k6 Ri,ss

(A19)

∗
Rtot = Ra,ss+ Ri,ss

(A20)

∗
∗
)
k2 k4 k5 Ra,ss= k1 k3 k6 (Rtot − Ra,ss(A21)

Using the fact that
gives

which again can be organized as
∗
Ra,ss
=

k1 k3 k6
Rtot
k1 k3 k6 + k2 k4 k5(A22)

Derivation of the Transfer functions for motif M5
According to Table 1 there are six kj → ki substitutions that give perfect adaptationin both Ra∗ and Ri . In the following the individual substitutions and correspondingtransfer functions for perfect adaptation sites are presented in more detail.
1. Substitution k2 → k3 where k2 = αk3 produces the following transfer function
from ∆k3 (s) to ∆Ra∗ (s)H(s) =Ass Ri,ss ·s∆Ra∗ (s)=∆k3 (s)(s + αk3)(s + k3 Ass + k4 Iss )As can be seen, this transfer function has a zero in origo.

7

(A23)

2. Substitution k2 → k6 where k2 = αk6 produces the following transfer function
from ∆k6 (s) to ∆Ra∗ (s)∆Ra∗ (s)∆k6 (s)(k3 αRi,ss Ass + k4 Iss (Ri,ss − Rtot ))·s=−(s + αk6)(s + k3 Ass + k4 Iss )(s + k6 )H(s) =

k3 αk6 Ri,ss Ass + k4 αk6 Iss (Ri,ss − Rtot )+(s + αk6 )(s + k3 Ass + k4 Iss )(s + k6 )



(A24)

In order to show that this transfer function actually has a zero in origo, we haveto use Eq. A19 and Eq. A20 such that the real part of the denominator of Eq.A24
can be written as (using k2 = αk6, Eq. A17 and Eq. A16)
k5
k1∗
+ k4 k2 (Ri,ss − (Ra,ss+ Ri,ss ))k2
k6∗= k1 k3 k6 Ri,ss − k2 k4 k5 Ra,ss(A25)

k3 αk6 Ri,ss Ass + k4 αk6Iss (Ri,ss − Rtot ) = k3 k2 Ri,ss

Hence, by inserting Eq. A19 into Eq. A25, it is shown that the real part of thedenominator is zero, and the solution to n(s) = 0 in Eq. A24 is s = 0. Therefore,
the transfer function in Eq. A24 can be writtenH(s) =∆Ra∗ (s)(k3 αRi,ss Ass + k4 Iss (Ri,ss − Rtot ))·s=−∆k6 (s)(s + αk6 )(s + k3 Ass + k4 Iss )(s + k6 )(A26)which has a zero in origo.
3. Substitution k4 → k1 where k4 = αk1 produces the following transfer function
from ∆k1 (s) to ∆Ra∗ (s)∆Ra∗ (s)∆k1 (s)αIss (Ri,ss − Rtot )·s + k3 Ri,ss + k2 αIss (Ri,ss − Rtot )=(s + k2 )(s + k3 Ass + αk1 Iss )H(s) =(A27)

In order to show that this transfer function actually has a zero in origo, we haveto use Eq. A19 and Eq. A20 such that the real part of the denominator of Eq.A27
can be written as (using k5 = αk1 and Eq. A17)∗
k3 Ri,ss + k2 αIss (Ri,ss − Rtot ) = k3 Ri,ss + k2 αIss (Ri,ss − (Ra,ss+ Ri,ss ))∗= k3 Ri,ss − k2 αIss Ra,ss∗= k1 k3 Ri,ss − αk1 k2 Iss Ra,ss
k5∗= k1 k3 k6 Ri,ss − k2 k4 k6 Ra,ss
k6∗= k1 k3 k6 Ri,ss − k2 k4 k5 Ra,ss

8

(A28)

Hence, by inserting Eq. A19 into Eq. A28, it is shown that the real part of thedenominator is zero, and the solution to n(s) = 0 in Eq. A27 is s = 0. Therefore,
the transfer function in Eq. A27 can be written∆Ra∗ (s)αIss (Ri,ss − Rtot )·s
H(s) ==∆k1 (s)(s + k2 )(s + k3 Ass + αk1Iss )(A29)

which has a zero in origo.
4. Substitution k4 → k6 where k4 = αk6 produces the following transfer function
from ∆k6 (s) to ∆Ra∗ (s)αIss (Ri,ss − Rtot )s∆Ra∗ (s)=H(s) =∆k6 (s)(s + k3 Ass + αk6 Iss )(s + k6 )(A30)

As can be seen, this transfer function has a zero in origo.
5. Substitution k5 → k1 where k5 = αk1 produces the following transfer function
from ∆k1 (s) to ∆Ra∗ (s)∆Ra∗ (s)∆k1 (s)(k3 Ri,ss + αk4 (Ri,ss − Rtot ))·s + k3 k6 Ri,ss + αk2 k4 (Ri,ss − Rtot )=(s + k2 )(s + k3 Ass + k4 Iss )(s + k6 )(A31)

H(s) =In order to show that this transfer function actually has a zero in origo, we haveto use Eq. A19 and Eq. A20 such that the real part of the denominator of Eq.A31
can be written as (using k5 = αk1)∗
k3 k6 Ri,ss + αk2 k4 (Ri,ss − Rtot ) = k3 k6 Ri,ss + αk2 k4 (Ri,ss − (Ra,ss+ Ri,ss ))∗= k3 k6 Ri,ss − αk2 k4 Ra,ss∗= k1 k3 k6 Ri,ss − αk1k2 k4 Ra,ss∗= k1 k3 k6 Ri,ss − k2 k4 k5 Ra,ss(A32)

Hence, by inserting Eq. A19 into Eq. A32, it is shown that the real part of thedenominator is zero, and the solution to n(s) = 0 in Eq. A31 is s = 0. Therefore,
the transfer function in Eq. A31 can be writtenH(s) =∆Ra∗ (s)(k3 Ri,ss + αk4 (Ri,ss − Rtot ))·s=∆k1 (s)(s + k2 )(s + k3 Ass + k4 Iss )(s + k6 )which has a zero in origo.

96. Substitution k5 → k3 where k5 = αk3 produces the following transfer function
from ∆k3 (s) to ∆Ra∗ (s)H(s) =Ass Ri,ss ·s + Ass Ri,ss k6 + k4 α(Ri,ss − Rtot )∆Ra∗ (s)=∆k3 (s)(s + k3 Ass + k4 Iss )(s + k6 )(A33)

In order to show that this transfer function actually has a zero in origo, we haveto use Eq. A19 and Eq. A20 such that the real part of the denominator of Eq.A33
can be written as (using k5 = αk3 and Eq.A16)∗
Ass Ri,ss k6 + k4 α(Ri,ss − Rtot ) = Ass Ri,ss k6 + k4 α(Ri,ss − (Ra,ss+ Ri,ss ))∗= Ass Ri,ss k6 − k4 αRa,ssk1∗= Ri,ss k6 − k4 αRa,ss
k2∗= k1 k3 k6 Ri,ss − k2 k4 αk3 Ra,ss∗= k1 k3 k6 Ri,ss − k2 k4 k5 Ra,ss(A34)

Hence, by inserting Eq. A19 into Eq. A34, it is shown that the real part of thedenominator is zero, and the solution to n(s) = 0 in Eq. A33 is s = 0. Therefore,
the transfer function in Eq. A33 can be writtenH(s) =Ass Ri,ss ·s∆Ra∗ (s)=∆k3 (s)(s + k3 Ass + k4 Iss )(s + k6 )which has a zero in origo.

10(A35)