arXiv:1802.01806v1 [q-bio.MN] 6 Feb 2018

Detection of sustained signals and its relation to
coherent feedforward loops
Chun Tung Chou
School of Computer Science and Engineering,
University of New South Wales,
Sydney, Australia.
E-mail: c.t.chou@unsw.edu.au
February 7, 2018
Abstract
Many studies have shown that cells use temporal dynamics of signaling molecules to encode information. One particular class of temporal
dynamics is sustained and transient signals, i.e. signals of long and short
durations respectively. This paper formulates a detection problem for
distinguishing sustained signals from transient ones. The solution of the
detection problem is to compute the likelihood ratio of observing a sustained signal to a transient signal. We show that, if the logarithm of this
likelihood ratio is positive (i.e. when the signal is likely to be sustained),
then this log-likelihood ratio can be approximately computed by a coherent feedforward loop. Although the capability of coherent feedforward
loops to discriminate sustained signals is known, its statistical information
processing interpretation has not been pointed out before.

Living cells use many different strategies to encode information. A strategy
is to encode information using the temporal dynamics of signaling molecules
[9, 2]. One particular class of temporal dynamics is sustained (also known as
persistent) and transient signals, i.e. signals of long and short durations respectively. The study in [8] found that PC-12 cells proliferated after a transient
ERK activation but differentiated after a sustained ERK activation. This shows
that the duration of ERK signaling affects the cell fate. Ref. [7] shows that coherent feedforward loops (FFL) can act as persistence detectors to differentiate
sustained signals from transient ones. A particular type of FFL is the coherent
type-1 feedforward loop with an AND logic at the output, or C1FFL for short.
Fig. 1a depicts the structure of the C1FFL which consists of a short arm, a long
arm and an AND logic combining the signals of the two arms. By modelling
the reaction dynamics of C1FFL with ordinary differential equations (ODEs),
Ref. [7] shows that a short and a long pulse applied to the input of the C1FFL,
give, respectively a zero and a non-zero output. Instead of starting from the
structure of C1FFL and show that it can detect sustained signals, we ask the
1

Input s(t)

Input s(t)
x*(t)
X
x*(t)

X

x*(t)

Y

Loglikelihood
ratio
computation

y(t)
AND

Log-likelihood
ratio

Output z(t)

(a)

(b)

Figure 1: (a) The coherent type-1 feedforward loop with AND logic. (b) The
detection theory framework.
question in the opposite direction. In this Letter, we start from the requirement
of detecting a sustained signal and use detection theory (DT) [6] to determine
the detector. We then show that the resulting detector can be approximately
realized by a C1FFL.
DT is a branch of statistical signal processing. Its aim is to use the measured data to decide whether an event of interest has occurred. For example,
DT is used in radar signal processing to determine whether a target is present or
not. In the context of this Letter, the events are whether the signal is transient
or sustained. A detection problem is often formulated as a hypothesis testing
problem, where each hypothesis corresponds to a possible event. Let us consider
a detection problem with two hypotheses, denoted by H0 and H1 , which correspond to respectively, the events of transient and sustained signals. Our aim is
to decide which hypothesis is more likely to hold. We define the log-likelihood
ratio (LLR) R:


P[measured data|H1 ]
R = log
(1)
P[measured data|H0 ]
where P[measured data|Hi ] is the conditional probability that the measured
data is generated by the signal specified in hypothesis Hi . Intuitively, if the
LLR R is positive, then the measured data is more likely to have been generated
by a sustained signal or hypothesis H1 , and vice versa. Therefore, the key idea
of DT is to use the measured data to compute the LLR and then use it to make
a decision.
We will now present a big picture explanation of how we will connect DT
with C1FFL. The signal x∗ (t) in Fig. 1a is the output signal of Node X in the
C1FFL. We can view the C1FFL as a 2-stage signal processing engine. In the
first stage, the input signal is processed by Node X to obtained x∗ (t) and this
is the part within the dashed box in Fig. 1a. In the second stage, the signal
2

x∗ (t) is processed by the rest of the C1FFL to produce the output signal z(t).
We will now make a connection to DT. We apply DT to the dashed box in
Fig. 1a. We consider x∗ (t) as the measured data and use them to determine
whether the input signal is transient or sustained. DT tells us that we should
use x∗ (t) to compute the LLR. This means that we can consider the 2-stage
signal processing depicted in Fig. 1b where the input signal generates x∗ (t) and
the measured data x∗ (t) are used to calculate the LLR. If we can identify the
LLR calculation in Fig. 1b with the processing by the part of C1FFL outside of
the dashed box, then we can identify the signal z(t) with the LLR.
We will now define the problem for detecting a sustained signal using DT.
Our first step is to specify the signaling pathway in Node X, which consists of
three chemical species: signaling molecule S, molecular type X in inactive form
and its active form X* . The activation and inactivation reactions are:
k+

S + X −−→ S + X∗
k−

X∗ −−→ X

(2a)
(2b)

where k+ and k− are reaction rate constants. Let x(t) and x∗ (t) denote, respectively, the number of X and X* molecules at time t. Note that both x(t)
and x∗ (t) are piecewise constant because they are molecular counts. We assume
that x(t) + x∗ (t) is a constant for all t and we denote this constant by M .
We assume that the input signal s(t), which is the concentration of the
signaling molecules S at time t, is a deterministic signal. We also assumed that
the signal s(t) cannot be observed, so any characteristics of s(t) can only be
inferred.
We model the dynamics of the chemical reactions by using chemical master
equation [4]. This means that x∗ (t) is a realisation of a continuous-time Markov
chain. This also means that the same input signal s(t) can result in different
x∗ (t).
The measured datum at time t is x∗ (t). However, in the formulation of
the detection problem, we will assume that at time t, the data available to
the detection problem are x∗ (τ ) for all τ ∈ [0, t]; in other words, the data are
continuous in time and are the history of the counts of X* up to time t. We will
use X∗ (t) to denote the continuous-time history of x∗ (t) up to time t. Note that
even though we assume that the entire history X∗ (t) is available for detection,
we will show that the calculation of the LLR at time t does not require the
storage of the past history.
The last step in defining the detection problem is to specify the hypotheses
Hi (i = 0, 1). Later on, we will identify H0 and H1 with, respectively, transient
and sustained signals. However, at this stage, we want to solve the detection
problem in a general way. We assume that the hypothesis H0 (resp. H1 ) is
that the input signal s(t) is the signal c0 (t) (c1 (t)) where c0 (t) and c1 (t) are
two different deterministic signals. Intuitively, the aim of the detection problem
is to use the history X∗ (t) to decide which of the two signals c0 (t) and c1 (t) is
more likely to have produced the observed history.

3

Based on the definition of the detection problem, the LLR L(t) at time t is
given by:


P[X∗ (t)|H1 ]
L(t) = log
(3)
P[X∗ (t)|H0 ]
where P[X∗ (t)|Hi ] is the conditional probability of observing the history X∗ (t)
given hypothesis Hi . We show in the Appendix that the L(t) obeys the following
ODE:




dL(t)
dx∗ (t)
c1 (t)
=
log
−
dt
dt
c0 (t)
+
k+ (M − x∗ (t))(c1 (t) − c0 (t))

(4)

where [a]+ = max(a, 0). We also assume that the two hypotheses are equally
likely, so L(0) = 0. Since x∗ (t) is a piecewise constant function counting the
number of X* molecules, its derivative is a sequence of Dirac deltas at the time
instants that X is activated or X* is deactivated. Note that the Dirac deltas
corresponding to the activation of X carries a positive sign and the [ ]+ operator
keeps only these. Note that a special case of Eq. (4) with constant ci (t) and
M = 1 appeared in [10]. A more general form of Eq. (4) which includes the
diffusion of signaling molecules can be found in [3].
The importance of Eq. (4) is that, given the measured data x∗ (t), we can use
Eq. (4) together with ci (t) to compute the LLR L(t). If s(t) is similar to c1 (t)
(resp. c0 (t)), then the LLR L(t) generally increases (decreases) over time and
becomes more positive (negative). If our aim is to distinguish a sustained signal
from a transient one accurately, then we want the sustained signal to produce
a large positive L(t). Since the positive contribution of L(t) comes from the
first term in Eq. (4), we can get a large positive L(t) by making sure that a
sustained signal will produce many activations. This occurs when a sustained
signal has a duration which is long compared to time scale of the activation and
deactivation reactions (2) and we will make use of this condition later.
The LLR L(t) can certainly be computed numerically by a computer but
a more interesting question is whether it can be computed by using a set of
chemical reactions. We will address this next.
In order to build a connection with C1FFL, we need to specify the form of
ci (t). We now assume that the transient signal c0 (t) and sustained signal c1 (t)
are of the form of a rectangular pulse. The durations for c0 (t) and c1 (t) are,
respectively, d0 and d1 with d1 > d0 . We define two concentration levels for the
signaling molecules: basal concentration a0 and reference concentration a1 , with
a1 > a0 . The temporal profile of ci (t) is: for 0 ≤ t < di , ci (t) = a1 , otherwise
ci (t) = a0 . Furthermore, we assume the input s(t) is a pulse of duration d. The
temporal profile of s(t) is: for 0 ≤ t < d, s(t) = a, otherwise s(t) = a0 .
There are two difficulties why the computation in Eq. (4) cannot be carried
out by chemical reactions. These are: (1) The LLR can take any real value but
chemical concentration can only be non-negative; (2) It is difficult to calculate
derivatives and M − x∗ (t) using chemical reactions. Instead of using L(t), we
4

will derive an approximation L̂(t) which has the properties: L̂(t) ≈ L(t) for
sustained signals and L̂(t) = 0 for transient signals.
We first consider the case when the input s(t) is a sustained signal which
results in a positive L(t). We assume that d is long compared with the time
scale of the activation and deactivation reactions. Since these reactions are fast,
we can assume quasi-steady state:
x∗ (t) ≈

M k+ s(t)
k+ s(t) + k−

(5)

By using the form of ci (t) and s(t) in Eq. (4), we see that, for t ∈ [d0 , min{d, d1 }],
the contribution of the first term on the right-hand side (RHS) of (4) to L(t)
equals to the number of times that X is activated in the time interval t ∈ [d0 , t].
If t is sufficiently large and if the activation and deactivation reactions are fast,
then a large number of activations take place. In this case, we can get approximately the same L(t) by approximating the positive derivative of x∗ (t) by its
mean value:


dx∗ (t)
≈ k+ (M − x∗ (t))s(t)
(6)
dt
+
After substituting Eqs. (5) and (6) into Eq. (4), and after some manipulation,
we arrive at:
dL̂(t)
= x∗ (t) × {k− π(t) [φ(s(t))]+ } where
dt
 
a1
a1 − a0
φ(u) = log
,
−
a0
u

(7)
(8)

and π(t) = 1 for d0 ≤ t < d1 and zero otherwise. We first verify that L̂(t) ≈ L(t)
for t ∈ [0, min{d, d1 }) for long d and d1 . We use the Stochastic Simulation
Algorithm (SSA) [5] to obtain 100 realisations of x∗ (t). We then use these time
series and Eq. (4) to compute the mean and standard deviation of L(t). We
numerically integrate Eq. (7) to obtain L̂(t) for comparison. We use k+ = 0.02,
k− = 0.5, d0 = 5, d1 = 60, a0 = 0.25, a1 = a = 10.7. The top and bottom
plots in Fig. 2a show the results for d = 70 and d = 40. We can see that the
approximation is good. In general, the approximation is good if d and d1 are
long compared to the time scale of the activation and deactivation reactions.
We note from Fig. 2a that when the approximation holds, the maximum value
of L̂(t) and L(t) match very well.
There are a few noteworthy features on the RHS of Eq. (7). The RHS
consists of the multiplication of x∗ (t) with some other functions. We will map
x∗ (t) to the short arm of the C1FFL and the other functions to the long arm.
The multiplication is to implement the AND logic. Next, we note that the
function π(t) is zero in the time interval [0, d0 ) which means the RHS of Eq. (7)
is zero t < d0 independent of the input. This shows how the detector is making
use of the prior knowledge of c0 (t) by not making use of any observations before
time d0 .
5

The form of Eq. (7) appears to be robust to how we define Hi . Instead of
assuming that Hi is a pulse of a specific duration di , we consider an alternative where we define hypothesis H0 (resp. H1 ) to be pulses whose duration is
uniformly distributed in [0, d0 ] ([d0 , d1 ]). These are composite hypotheses and
can be handled by using the Bayesian approach in [6, p.198]. We can show that
these hypotheses give an ODE in the same form as Eq. (7) except that π(t) is
slightly different. Thus the mapping to C1FFL also holds for these alternative
hypotheses.
Having shown that L̂(t) ≈ L(t) for long pulses, we now show that L̂(t) = 0
for short pulses. For a small d, the LLR L(t) becomes negative. In Eq. (7), the
[ ]+ operator is included so that L̂(t) ≥ 0 for all t. In particular, for transient
inputs with duration d < d0 , L̂(t) = 0 ∀t ≥ 0. This is consistent with the
behaviour of an ideal C1FFL which gives a zero output for short pulses [1].
Our next step is to show that Eq. (7) can be approximately realized by
the following reaction system, which models Node Y and the AND logic in the
C1FFL in Fig. 1a:
dy(t)
dt

=

hy x∗ (t)ny
ny
Ky + x∗ (t)ny
|

dz(t)
dt

=

{z

Hy (x∗(t))

x∗ (t) ×

−dy y(t)

(9a)

}

hz y(t)nz
Kznz + y(t)nz
{z
}
|

(9b)

Hz (y(t))

where hy , ny etc. are coefficients of the Hill functions. Our aim is to make
z(t) in Eq. (9b) to be approximately equal to L̂(t), and this can be achieved by
approximating k− π(t)[φ(s(t))]+ (= η(t)) in Eq. (7) by Hz (y(t)) in Eq. (9b). If
the input s(t) is a long pulse, then the time profiles of both η(t) and y(t) contain
a period of time that they plateau. The plateau in η(t) contributes to the ramplike increase in L̂(t) in Fig. 2a. This means we need to match the values of η(t)
and Hz (y(t)) at their plateau. For a pulse s(t) with an amplitude a, the heights
of the plateau of η(t) and Hz (y(t)) are, respectively, k− [φ(a)]+ (= f1 (a)) and
k+ a
Hz ( d1y Hy ( kM
))(= f2 (a)), and we want f1 (a) ≈ f2 (a) for as large a range of
+ a+k−
a as possible. Note that for all a such that f1 (a) > 0, both functions f1 (a) and
f2 (a) are strictly increasing, strictly concave and both f1 (∞) and f2 (∞) are
constants. Therefore, we can choose the Hill function coefficients to fit f2 (a) to
f1 (a). This argument takes care of the case when s(t) is a long pulse. For short
pulses, we need to realise π(t) whose purpose is to make the initial part of L̂(t)
zero. This is a feature shared by the ideal C1FFL model in [1]. Ref. [1] shows
that this can be realized by choosing a big enough Kz in Eq. (9b) so that the
production rate of z(t) is small initially.
We now present two numerical examples. We use the same k+ , k− , M ,
a0 and a1 values as before. We choose d0 = 10, d1 = 80. We use parameter
estimation to determine the reaction constants in Eq. (9) so that the C1FFL
output z(t) matches L̂(t) for a range of a. Fig. 2b compares L̂(t) and z(t) for
6

input s(t) with d = 10, 30 and 70 for a = 10.7. When d = 10, the output of the
C1FFL is small. For d = 30 and 70, the C1FFL output matches well with the
approximate LLR L̂(t). In the second example, we use inputs of a fixed duration
of 70 but vary the amplitude a from 2.7 to 85.7. Fig. 2b compares L̂(t) and
z(t) at t = 70. It can be seen that the C1FFL approximation works for a wide
range of a. These examples show that we can match the C1FFL output z(t)
to approximate LLR L̂(t) for input s(t) of different durations and amplitudes.
Although we have only shown that z(t) matches L̂(t) for pulse inputs, the match
also extends to other slowly-varying inputs.
Conclusions: This Letter considers the problem of detecting sustained signals
from transient ones. It makes use of detection theory and turns the problem into
one of computing the log-likelihood ratio. It further shows that the log-likelihood
ratio can be approximately computed by a coherent type-1 feedforward loop.
Although the capability of coherent feedforward loops to discriminate sustained
signals is known, its statistical information processing interpretation has not
been pointed out before.

A

Proof of (4)

In order to derive (4), we consider the history X∗ (t + ∆t) as a concatenation of
X∗ (t) and x∗ (t) in the time interval (t, t+∆t]. We assume that ∆t is chosen small
enough so that no more than one activation or deactivation reaction can take
place in (t, t + ∆t]. Given this assumption and right continuity of continuoustime Markov Chain, we can use x∗ (t + ∆t) to denote the history in (t, t + ∆t].
By using the Markov property, we can show that:


P[x∗ (t + ∆t)|H1 , x∗ (t)]
(10)
L(t + ∆t) = L(t) + log
P[x∗ (t + ∆t)|H0 , x∗ (t)]
By using the propensity functions of the reactions, we have:
P[x∗ (t + ∆t)|Hi , x∗ (t)] =
δx∗ (t+∆t),x∗ (t)+1 k+ (M − x∗ (t)) ∆t ci (t)(t)+
δx∗ (t+∆t),x∗ (t)−1 k− b(t) ∆t +
δx∗ (t+∆t),x∗ (t) (1 − k+ (M − x∗ (t))ci (t) ∆t − k− b(t) ∆t)

(11)

where δa,b is the Kronecker delta which is 1 when a = b. By substituting Eq. (11)
into Eq. (10) and taking the limit ∆t → 0, we have after some manipulations:


δx∗ (t+∆t),x∗ (t)+1
dL(t)
c1 (t)
= lim
log
−
∆t→0
dt
∆t
c0 (t)
δx∗ (t+∆t),x∗ (t) k+ (M − x∗ (t)) (c1 (t) − c0 (t))

(12)

In order to obtain Eq. (4), we use the following reasonings. First, the term
δ
∗ (t)+1
lim∆t→0 x∗ (t+∆t),x
is a Dirac delta at the time instant that an X molecule
∆t
7

is activated. Second, the term δx∗ (t+∆t),x∗ (t) is only zero when the number of X∗
molecule changes but the number of such changes is countable. In other words,
δx∗ (t+∆t),x∗ (t) = 1 with probability one. This allows us to drop δx∗ (t+∆t),x∗ (t) .

References
[1] U. Alon. An Introduction to Systems Biology: Design Principles of Biological Circuits. Chapman & Hall, 2006.
[2] Marcelo Behar and Alexander Hoffmann. Understanding the temporal
codes of intra-cellular signals. Current opinion in genetics & development,
20(6):684–693, December 2010.
[3] Chun Tung Chou. Maximum a-posteriori decoding for diffusion-based
molecular communication using analog filters. Nanotechnology, IEEE
Transactions on, 14(6):1054–1067, 2015.
[4] C. Gardiner. Stochastic methods. Springer, 2010.
[5] D Gillespie. Exact stochastic simulation of coupled chemical reactions. The
journal of physical chemistry, 1977.
[6] Steve M. Kay. Fundamentals of Statistical Signal Processing, Volume II:
Detection Theory. Prentice Hall, 1998.
[7] S Mangan and U Alon. Structure and function of the feed-forward loop
network motif. Proceedings of the National Academy of Sciences of the
United States of America, 100(21):11980–11985, October 2003.
[8] C J Marshall. Specificity of receptor tyrosine kinase signaling: transient versus sustained extracellular signal-regulated kinase activation. Cell, 80:179–
185, January 1995.
[9] Jeremy E Purvis and Galit Lahav. Encoding and Decoding Cellular Information through Signaling Dynamics. Cell, 152(5):945–956, February 2013.
[10] Eric D Siggia and Massimo Vergassola. Decisions on the fly in cellular sensory systems. Proceedings of the National Academy of Sciences,
110(39):E3704–12, September 2013.

8

LLR

2500

1250

0
0

20

40

60

80

60

80

time

LLR

1600

800

0
0

20

40

time

(a) Comparing L(t) to L̂(t).

Approx LLR / C1FFL output

1000

500

0

0

40

80

time

Approx LLR / C1FFL output

(b) Comparing L̂(t) to C1FFL output.

8000

6000

4000

2000

0

0

20

40

60

80

(c) Comparing L̂(t) and C1FFL output for
different pulse input amplitude a.

Figure 2: Numerical results. (Best view in color.)
9

