Using Biochemical Circuits to Approximately Compute Log-Likelihood Ratio for Detecting Persistent Signals

Given that biochemical circuits can process information by using analog computation, a question is: What can biochemical circuits compute? This paper considers the problem of using biochemical circuits to distinguish persistent signals from transient ones. We define a statistical detection problem over a reaction pathway consisting of three species: an inducer, a transcription factor (TF) and a gene promoter, where the inducer can activate the TF and an active TF can bind to the gene promoter. We model the pathway using the chemical master equation so the counts of bound promoters over time is a stochastic signal. We consider the problem of using the continuous-time stochastic signal of the counts of bound promoters to infer whether the inducer signal is persistent or not. We use statistical detection theory to derive the solution to this detection problem, which is to compute the log-likelihood ratio of observing a persistent signal to a transient one. We then show, using time-scale separation and other assumptions, that this log-likelihood ratio can be approximately computed by using the continuous-time signals of the number of active TF molecules and the number of bound promoters when the input is persistent. Finally, we show that the coherent feedforward gene circuits can be used to approximately compute this log-likelihood ratio when the inducer signal is persistent.


I. INTRODUCTION
The fact that living organisms use biomolecular circuits for information processing has provided much inspiration for engineers to design new biomolecular circuits.Firstly, engineers have been inspired to use engineering theory to design synthetic circuits.E.g., control theory has been used to design circuits for homeostatic control [1], concentration regulation [2] and to counter the effect of loading [3]; signal processing theory has been used to design molecular circuits that can be used for communications [4], [5], [6] and filtering [7].Secondly, one can view the information processing carried out by biomolecular circuits as analog computation [8], [9].This view point has helped engineers to design molecular circuits that can perform logarithmic sensing [10], parity check [11] and integral control [12].This paper considers the problem of persistence detection, i.e. the problem of deciding whether a particular chemical species has been present in sufficient quantity for a long enough time.In the natural world, the bacteria Escherichia coli (E.coli) perform persistence detection [13].In the synthetic world, cell-based therapy [14] can make use of persistence detection on biomarkers.In this paper, we will solve the persistence detection problem by using statistical detection theory [15] and show that a gene circuit (which is specific type of biochemical circuit) can be used to approximately compute the solution of this detection problem.
A gene circuit that can perform persistence detection in E. coli is the Coherent Type-1 Feedforward Loop with an AND logic at the output [16] or the C1-FFL for short.The C1-FFL is a network motif and is a frequently found circuit in both E. coli and Saccharomyces cerevisiae (yeast) [17], [16].This means that the C1-FFL carries out important functions in cells.The authors in [16] showed that the C1-FFL can act as a persistence detector.They did this by modelling the gene expression in the C1-FFL by using ordinary differential equations (ODEs) and show that a persistent (resp.transient) input to the C1-FFL will result in a high (zero) output.
The paper [16] took a deterministic approach to understand persistence detection.Given that the biochemical environment is stochastic, it is therefore necessary to understand how cells can infer information on the environment from a stochastic point of view [18].This paper considers a reaction pathway consisting of three chemical species: an inducer, a transcription factor (TF) and a gene promoter.In this reaction pathway, the inducer can activate the TF and the activated TF can bind with the gene promoter.In order to model a stochastic biochemical environment, we model the reaction pathway using the chemical master equation [19].We consider a detection problem whose aim is to infer whether the inducer signal is persistent or not by using the signal of the number of bound promoters over time.According to detection theory, the solution to this detection problem is to compute a loglikelihood ratio and we derive an ODE which describes the evolution of this log-likelihood ratio over time.In order to connect this ODE to the C1-FFL, we use time-scale separation and other assumptions to derive an intermediate approximation which is an ODE that can approximately compute the loglikelihood ratio for persistent signals.We then show that this intermediate approximation can be realised by using a C1-FFL.The key contribution of this paper is to show that it is possible to find the parameters of a C1-FFL so that its output is approximately equal to the log-likelihood ratio of statistical detection problem when the inducer signal is persistent.More specifically, the aim of this statistical detection problem is to detect whether the inducer signal is persistent in an inducer-TF-gene pathway.In addition, the methodology in this paper can be useful for designing synthetic molecular circuits for performing other signal processing tasks.
This paper makes advances compared to our previous work [20].In comparison to [20], this paper makes two different assumptions: (i) This paper considers an inducer-TF-gene pathway but [20] considered only an inducer-TF pathway; (ii) This paper assumes that the inducer signal is stochastic while [20] assumed that the inducer signal is deterministic.These two different assumptions mean that new methodologies are needed to show that the C1-FFL can be used to approximately compute the log-likelihood ratio of a persistence detection problem.First, we need to show how the log-likelihood ratio can be computed exactly (Sec.III-B) in an inducer-TFgene pathway.Second, we need to derive a new method to approximately compute the log-likelihood ratio (Sec.IV).In particular, this paper needs to approximate the solution to a Bayesian filtering problem [21] but this is not required in [20] as it considered a deterministic inducer signal.Third, we need to derive a method to show how the parameters of the approximate log-likelihood ratio computation can be mapped to the C1-FFL parameters (Sec.V).These three aspects are the new elements of this paper in comparison to [20].Furthermore, in comparison with our earlier conference paper [22], this paper explains why a C1-FFL can be used to approximately compute the approximate log-likelihood ratio (Sec.V) and provide full derivation on the computation of exact and approximate log-likelihood ratios (Appendices A and B).
The rest of this paper is organised as follows.Sec.II presents background information on the C1-FFL.We then define the detection problem and present its solution in Sec.III.After that in Sec.IV, we present a method to approximately compute the log-likelihood ratio and use this approximation in Sec.V to show that the C1-FFL can be used to approximately compute the the log-likelihood ratio when it is positive.Finally, Sec.VI presents a discussion and concludes the paper.

II. BACKGROUND ON C1-FFL
The C1-FFL can be depicted as a network where each link is associated with a signal and each node transforms the input signal(s) into an output signal.Fig. 1 shows the network of the C1-FFL.The input signal is s(t) and output signal is z(t).Both x * (t) and y * (t) are intermediate signals, and s y (t) is an external signal.
The C1-FFL in Fig. 1 is an abstraction of the molecular interactions which are depicted in Fig. 2. In the figure, both S and S y are inducers.Both X and Y are TFs, which are expressed by their corresponding gene.The inducer S (resp.S y ) turns the inactive form X (Y) into its active form X * (Y * ).The activation of gene Z requires the binding of both X * and Y * to the promoter of Z, i.e. the AND gate in Fig. 1.
Note that there is a one-to-one correspondence between the chemical species in Fig. 2  t and so on.In this paper, we will assume that the inducer S y is always present and its concentration is always above the threshold needed to activate Y. Furthermore, we assume the activation of Y by S y is fast, this allows us to write y * (t) = y(t) and we will use y(t) for y * (t) from now on.By using Hill function to model the gene expression, [13] presents an ODE model for the C1-FFL, as follows: on the right-hand side (RHS) of (1c) implements the AND gate in Fig. 1.With suitably chosen parameter values, the C1-FFL in (1) acts as a persistence detector in the sense that if the input signal s(t) is a persistent (resp.transient), then the output z(t) has a high (low) value.

III. STATISTICAL DETECTION ON A REACTION PATHWAY
Our aim is to consider a statistical detection problem to determine whether the input signal is persistent or not.However, in this section, we will consider a more general detection problem because it can readily be solved and we will specialise it to persistence detection in Sec.IV.This section is divided into two parts.We define the detection problem in Sec.III-A and present its solution in Sec.III-B.
Convention: In this paper, we use upper case letters to denote a chemical species, e.g.S, X * etc.For each chemical species, there are two corresponding continuous-time signals based on its concentration and molecular counts.E.g. for the chemical species X * , we denote its concentration over time as x * (t) (note: lower case x) and its molecular counts over time is X * (t) (note: upper case X).

A. Detection problem
In order that we can connect the detection problem to the C1-FFL later on, we will define the detection problem using a reaction pathway which is a subset of the C1-FFL species and reactions in Fig. 2. We have depicted the reaction pathway used in the detection problem in Fig. 3.The reaction pathway consists of five chemical species: S, inactive X and its corresponding active form X * , as well as inactive Z and the complex Z X * which is formed by the binding of X * to Z.These five species take part in the following four chemical reactions: where k + , k − , g + and g − are reaction propensity constants.
For the time being, we will make the simplifying assumption that the volume scaling needed to convert between propensity and reaction rate constants is 1.This simplification allows us to equate propensity constants with reaction rate constants.We will explain how non-unit volume can be dealt with in Remark 3.With this assumption, note that k + and k − in (2a) and (2b) are equal to those in (1a).
In terms of molecular biology, S is an inducer and X is a TF.In Reaction (2a), the species S activates X to produce X * .Reaction (2b) is a deactivation reaction.The reactions (2a) and (2b) are depicted in both Figs. 2 and 3.
The species Z is a gene.In fact, Z in Fig. 3 is the same as Z in Fig. 2. Note that Fig. 2 follows the standard convention in molecular biology where a gene and the protein that it expresses are given the same symbol Z.However, in this paper, we need different symbols for the gene and the protein that the gene expresses so that we can clearly distinguish their corresponding time signals.Therefore, we have chosen to use Z to denote the gene and use Z to denote the protein expressed by Z.In Reaction (2c), an active X * binds with the promoter of Z to produce the complex Z X * .Lastly, Reaction (2d) is an unbinding reaction.
Let S(t), X(t), X * (t), Z(t) and Z X * (t) denote, respectively, the number of S, X, X * , Z and Z X * molecules at time t.Note these signals are piecewise constant because they are molecular counts.We assume that X(t) + X * (t) (resp.Z(t) + Z X * (t)) is a constant for all t and we denote this constant by M (N ).
We will refer to S(t) as the input signal.The goal of the detection is to use the signal Z X * (t) to determine whether the input is persistent or not.We assume that the signal S(t) is generated by some species and chemical reactions upstream of (2).We will model these upstream reactions and (2) using the chemical master equation [19] which is a specific type of continuous-time Markov chain.We further assume that the upstream species do not react with X, X * , Z and Z X * .Intuitively, this means we can predict the behaviour of X, X * , Z and Z X * from that of S.
We have now defined the reaction pathway and its model.The next task is to specify the measured data and the hypotheses for the detection problem.The measured datum at time t is Z 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 Z X * (τ ) for all τ ∈ [0, t]; in other words, the data are continuous in time and are the history of the counts of Z X * up to time t inclusively.We will use Z X * (t) to denote the continuous-time history of Z X * (t) up to time t inclusively.
We now specify the hypotheses H i (i = 0, 1) for the detection problem.Later on, we will identify H 0 and H 1 with, respectively, transient and persistent signals.However, at this stage, we want to solve the detection problem in a general way.We assume that H 0 and H 1 are two distinct subsets of the set of all possible S(t).Intuitively, the aim of the detection problem is to decide which signal class H 0 or H 1 is more likely to have produced the observed history.
We remark that in the definition of the detection problem, the input signal S(t) is not directly observable.Since S reacts with the molecules in the reaction pathway in Fig. 3, the downstream signal Z X * (t) contains information on S(t).The aim of the detection problem is to infer the information on S(t) from this downstream signal.Given that we model the chemical system with the chemical master equation, both signals S(t) and Z X * (t) are noisy.

B. Solution to the detection problem
The aim of the detection problem is to decide which hypothesis H i (i = 0, 1) is more likely to have generated the observed history Z X * (t).We assume that the two hypotheses .
are equally likely, therefore the log-ratio of posteriori probabilities ] is equal to the log-likelihood ratio L(t): where P[ Z X * (t)|H i ] is the conditional probability of observing the history Z X * (t) given hypothesis H i .
In Appendix A, we show that the time evolution of L(t) is given by the following ODE: where [w] + = max(w, 0), E[ ] denotes the expectation and E[X * (t)| Z X * (t), H i ] is the conditional expectation of X * (t) given the history and H i .Note that in deriving (4), we assume that the hypotheses have been properly chosen so that J i (t) > 0 (for t > 0 and i = 0, 1) which in turn implies that log J1(t) J0(t) is well defined.
Since Z X * (t) is a piecewise constant function counting the number of Z X * molecules, its derivative is a sequence of Dirac deltas at the time instants that Z X * forms or unbinds.Note that the Dirac deltas corresponding to the formation of Z X * carries a positive sign and the [ ] + operator keeps only these.Fig. 4 shows an example Z X * (t) and its corresponding At the time instant t that Z X * is formed or unbind, the expectation J i (t) also has a jump in value at time t.The term J i (t − ) in (4) refers to the value of J i (t) just before the jump.Lastly, we assume that the two hypotheses are a priori equally likely, so L(0) = 0.
We next present a numerical example to illustrate the properties of (4) and to explain what information is important for persistence detection.This example will also provide some intuition on how we will approximately compute the loglikelihood in Sec.IV.
1) Numerical example: The aim of this example is to show that we can use the log-likelihood ratio computed by (4) to distinguish a persistent signal from a transient one.In order to conduct the numerical study, we will assume that the inducer S is produced and degraded by the following chemical reactions: where S pre is a precursor that produces S, and f + and f − are reaction rate constants.The reason why we choose to use these reactions is that they allow us to use the time profile S pre (t), which is the number of precursor molecules a time t, to control the amplitude and duration of S(t).The chemical master equation is used to model the reactions ( 6) and ( 2).We assume that both hypotheses H 0 and H 1 include the knowledge of the reactions ( 6) and ( 2) because this knowledge is required for computing E[X * (t)| Z X * (t), H i ] in (5).In order to define the hypotheses, we first define reference signals R 0 (t) and R 1 (t) for S pre (t).Both reference signals are deterministic ON-OFF pulses and their time profile has the form: where d i is the ON-duration for R i (t); and, S ON,ref pre and S OFF,ref pre , are, respectively, the ON and OFF amplitudes.We assume that S ON,ref pre S OFF,ref pre so that one may consider S OFF,ref pre as a reference value for basal concentration.Furthermore, we assume that the duration of R 1 (t) is longer than R 0 (t), i.e. d 1 > d 0 .We define H 0 (resp.H 1 ) as the set of all S(t) that are generated by using R 0 (t) (R 1 (t)) as the time profile for S pre (t).Since d 1 > d 0 , we will identify H 0 (resp.H 1 ) as, respectively, the set of transient (persistent) signals.
For this example, the kinetic parameters for the reaction pathway (2) are: k + = 0.02, k − = 0.5, g + = 0.002 and g − = 0.05.The total number of TFs M is 100 and the number of genes N is 1.Furthermore, for reactions (6), f + = 0.37 and f − = 0.1.For the reference signals, we choose S ON,ref pre and S OFF,ref pre to be, respectively, 100 and 3 molecules; and d 1 = 600 and d 0 = 100.
We will now use the above set-up to detect a persistent signal.For simplicity, we assume that the actual S pre (t) is the reference signal R 1 (t) which results in a persistent input S(t); see the top plot of Fig. 5a for a sample of persistent S(t).We first use the chosen S pre (t) and the Stochastic Simulation Algorithm (SSA) [23] to produce a realisation of Z X * (t).The simulation assumes that there are zero S, X * and Z X * molecules initially.We then use optimal Bayesian filtering [21] to numerically compute E[X * (t)| Z X * (t), H i ] for both i = 0, 1.The numerical solution to this optimal Bayesian filtering problem requires us to solve a master equation given the hypothesis.In order to solve this master equation exactly, we have purposely chosen the problem parameter values so that the number of S and X molecules are limited.After solving the Bayesian filtering problem, we numerically integrate (4) to obtain the log-likelihood ratio L(t), which is plotted as the solid blue line in Fig. 5b.We see that the log-likelihood ratio is zero for t ≤ 100, ramps up in the time interval [100, 600] and plateaus after t ≥ 600.The log-likelihood ratio reaches and J 1 (t) − J 0 (t) in ( 4).
a positive value at the end, which means correct detection because it says the Z X * (t) signal is more likely to have been generated by a persistent signal.
Next, we use a transient input.We assume that the actual S pre (t) is the reference signal R 0 (t); see the bottom plot of Fig. 5a for a sample of transient S(t).We perform the same steps as before, namely SSA simulation, optimal Bayesian filtering and numerical integration to obtain the log-likelihood ratio for this transient input.The resulting log-likelihood ratio L(t) is plotted as red dashed lines in Fig. 5b.This L(t) becomes negative which means correct detection.
Fig. 5c shows the weighting factors log J1(t) J0(t) and J 1 (t)− J 0 (t) in ( 4) for the case when the S(t) is persistent.(The curves are similar when S(t) is transient.)It shows that these two weighting factors are mostly positive in the time interval [100,600] but are zero outside.This means the contribution to the log-likelihood ratio comes from the signal within [100,600].This can also be seen from Fig. 5b where the log-likelihood ratio does not change outside of [100,600] but increases (resp.decreases) for persistent (transient) signal within [100,600].This makes intuitive sense because the persistent input is different from the transient input within this time interval, so the signal in this time interval is useful for discriminating persistent signals from transient ones.
If a persistent signal can give a large positive log-likelihood ratio, then the probability of correctly detecting the persistent signal is higher.For this example, the positive contribution to log-likelihood ratio comes from the first term on the RHS of (4) because the weighting factors are non-negative, see Fig. 5c.In fact, each time when a X * binds to a Z in the time interval [100,600], it creates a positive jump in the magnitude of the log-likelihood ratio, which can be seen in Fig. 5b.This means that a persistent signal becomes easier to detect if X * binds to Z many times when the signal is ON.This can be achieved if the ON duration of the persistent signal has a longer time-scale compared to those of the binding and unbinding reactions of Z X * (i.e.reactions (2c) and (2d)) so that these reactions occur many times when the input is ON.
Although this example shows that the solution (4) can distinguish a persistent signal from a transient one by computing the log-likelihood ratio, the solution is hard to implement in a biochemical environment because Bayesian filtering requires extensive computation and a model of the pathway.In the next section, we will explore how we can compute the loglikelihood ratio approximately without using Bayesian filtering.

APPROXIMATELY
Our ultimate goal is to show that the computation of the loglikelihood ratio L(t) in ( 4) can be carried out by a C1-FFL in (1), i.e. there exists a set of parameters for the C1-FFL such that z(t) in ( 1) is approximately equal to L(t) in (4).It is not obvious from the expression of (4) that this can be done.The aim of this section is to derive an ODE, which will be referred to as the intermediate approximation, such that the output of this ODE is approximately equal to L(t) when the input is persistent.We will then use this intermediate approximation in Sec.V to relate to the C1-FFL.

A. Assumptions
The detection problem and its solution in Sec.III are general in the sense that they apply to any reaction pathways of the form (2) and hypotheses H i .In order to connect the detection problem to the C1-FFL model in (1), we will need to make specific assumptions to derive the intermediate approximation.We will specify these assumptions in this subsection.
We make the following two assumptions on the reaction pathway (2): • The time-scale of the inducer-TF reactions (2a) and ( 2b) is faster than that of the TF-gene promoter reactions (2c) and (2d).This translates to small g + and g − .• The number of TF molecules M is much higher than the number of genes N .We believe these are realistic assumptions.First, according to [24, Table 2.2], for E. coli, the time-scale for equilibrium binding of small molecules to protein is of the order of 1 ms and the time-scale for TF binding to gene promoter is of the order of 1s.Second, the copy number of most genes is either 1 or 2.
In order to make the derivation of the intermediate approximation analytically tractable, we assume that hypotheses H i are defined in the same way as in the numerical example in Sec.III-B1.This is so that we can approximate some signal, e.g. the weighting functions in Fig. 5c, by a piecewise constant function.
For analytical tractability, we further assume that the actual S pre (t) is a deterministic ON-OFF pulse of the form: where d is the pulse duration and the ON-amplitude f− f+ α will result in a steady state E[S(t)] of amplitude α when the pulse is ON.
We assume that the input signal S(t) is produced by using S pre (t) with the reactions (6), and the reaction rate constants f + and f − have been chosen such that the signal S(t) is slowly time varying compared to the time-scale of the inducer-TF reactions (2a) and (2b).This makes intuitive sense because a persistent or long S(t) needs to be "measured" by faster reactions.
We assume that if a persistent input is applied to the reaction pathway (2), the pathway is almost at steady state by d 0 where d 0 is the duration of the reference signal R 0 (t) in (7).This requirement can be met if the duration d 0 is long enough.It may be instructive to recall from the discussion in the numerical example in Sec.III-B1 that there is a time interval which is informative for persistence detection.For the assumptions in this section, the informative time interval can be shown to be

B. Intermediate approximation
An ideal persistence detector has the properties that a transient input will result in a zero output and a persistent input will result in a positive output [24].The C1-FFL, when acting as a persistence detector, can be considered to be an approximation of this ideal behaviour [24].However, it is not possible to map the log-likelihood ratio detector in (4) to the C1-FFL because the log-likelihood ratio becomes negative for transient signals but the concentration in the C1-FFL can only be non-negative.We will use the intermediate approximation as a bridge to connect (4) to the C1-FFL.The intermediate approximation has two key properties.First, if the input is transient, then the output of the intermediate approximation is zero.Second, if the input is persistent, then the output of the intermediate approximation is approximately equal to the loglikelihood ratio given by (4).Another purpose of the intermediate approximation is to replace the complex computation in (4), e.g.derivative and optimal Bayesian filtering, by simpler computation that can be implemented by chemical reactions.
The derivation of the intermediate approximation is given in Appendix B, using the assumptions stated in Sec.IV-A.The derivation shows that the time evolution of the intermediate approximation L(t) is given by the following ODE: where Furthermore, it can be shown that time evolution of E[ L(t)] obeys the following ODE: The behaviour of the intermediate approximation L(t) depends on two parameters of the input signal S(t): its mean ON-amplitude α and duration d (see (8) approximately equal to the log-likelihood ratio L(t) in the time interval 0 ≤ t < min{d, d 1 }.Note in particular that this approximation holds for a range of d and α values.From now on, we will choose d 1 to be ∞ so that L(t) ≈ L(t) holds for 0 ≤ t < d, i.e. when the persistent signal is ON.Note that an infinite For the first numerical experiment, we use α = 37, which means S pre (t) has an amplitude of S ON,ref pre when it is ON.This excitation means there is a mean probability of 0.6 that X is active.We use SSA simulation to generate 100 realisations of X * (t) and Z X * (t), and use them to compute the true loglikelihood ratio L(t) (which requires only Z X * (t)) and the intermediate approximation L(t) (which requires both X * (t) and Z X * (t)).Fig. 6a compares one realisation of L(t) (blue line) and L(t) (red line) for a persistent input.It can be seen that the intermediate approximation is fairly accurate.
We then use all 100 realisations to compute the mean of L(t), the mean of L(t), and the root-mean-square (RMS) error of L(t) − L(t).The results are plotted in Fig. 6b.This shows that L(t) is a good approximation of L(t) for many realisations.Next, we check the accuracy of using (16) to compute E[ L(t)].The black dashed lines in Fig. 6b plot E[ L(t)] computed from (16), which is almost the same as the mean obtained via SSA.This shows that ( 16) is an accurate method to compute E[ L(t)].
The discussion so far focused on persistent signals.The black line in Fig. 6a shows the intermediate approximation for a transient input with d = 100.It can be seen that the output is small.This agrees with our earlier prediction.
For a given set of hypotheses, the intermediate approximation holds for a range of α and d.We now change α to 99.9 (which means S pre (t) has an amplitude of 1.33S ON,ref   pre   when it is ON).By using 100 rounds of SSA simulations, we compute the the means of both L(t) and L(t), as well as the RMS error between them.The results are in Fig. 6c and they show that the intermediate approximation is accurate for different values of α.Fig. 6c also shows E[ L(t)] computed from ( 16) and we can see that the accuracy is good.
Remark 1: The reader may wonder why we do not define the hypotheses of the detection problem as: H 0 (resp.H 1 ) means the duration of the input signal is shorter (longer) than a given threshold.The reason is that these are composite hypotheses and the solution to the resulting detection problem is much harder, see [20,Remark 5.1] for a more in-depth discussion.We also want to point out that, even with our simpler formulation, the resulting detector gives a small output for any signal whose duration is shorter than d 0 .

C. Relevance to biological detection problems
In this section, we will discuss why the likelihood ratio is a relevant criterion for some biological detection problems.For this exposition, we will assume that the purpose of the detection is to determine if a certain nutrient is persistently present and if yes, then the organism wants to produce the enzyme to consume the nutrient.Our discussion is based on Bayesian decision theory which considers both utility and cost of actions, e.g. the successful detection and consumption of a nutrient gives a positive utility to the organism at the cost of producing the enzyme.
We will use the hypotheses H 0 and H 1 to refer to the conditions that the nutrient is respectively, absent and present, in the environment.We will also use negative and positive, respectively, to refer to the two hypotheses H 0 and H 1 .Given these two hypotheses, the detection problem may decide for either H 0 (negative) or H 1 (positive).Table I summarises the four possible combinations from the two environmental conditions and the two detection outcomes.These four combinations are labelled as True Negative (TN), False Negative (FN), False Positive (FP) and True Positive (TP) using the terminologies commonly used in statistics.For the time being, we assume that the utility can only take two values: U 1 and U 0 where U 1 > U 0 = 0.The living organism can only get the positive utility U 1 for TP as this is the only situation which the nutrient is present and detected.Similarly, for the cost, we assume C 1 > C 0 = 0.The living organism incurs a cost when FP or TP occurs as enzymes are made in these cases.
Let P 0 and P 1 denote the true probabilities that the nutrient is, respectively, absent and present.Let also P T P denote the probability of TP etc.The mean utility is P T P P 1 U 1 and the mean cost is P F P P 0 C 1 + P T P P 1 C 1 .Let C max be the maximum cost that the living organism can afford.From a Bayesian decision point of view, the goal is to maximise the mean utility subject to a constraint on the mean cost, i.e. max P T P P 1 U 1 (17) subject to P F P P 0 C 1 + P T We show in Appendix C that the solution of this utility maximisation problem is to choose a suitable positive threshold ρ such that the organism should decide for H 1 if the likelihood ratio P[data|H1] P[data|H0] ≥ ρ.In fact, the derivation shows that, if U 1 > λC 1 where λ is the Lagrangian multiplier of the above optimisation problem, then the above utility maximisation problem is equivalent to maximising P T P subject to an upper bound on P F P , which is in fact the scenario covered by the Neyman-Pearson lemma.This shows that the likelihood ratio is a suitable statistic to be used.Note that we have assumed for simplicity that U 0 = 0 and C 0 = 0 earlier, however it can be shown that, under some conditions, the result still holds for non-zero U 0 and C 0 .
1) Numerical example: In this section, we want to study the performance of using the approximate log-likelihood ratio L(t) for decision making.Since our interest is to detect persistent signals, we can do that by testing whether the log-likelihood ratio is greater than or equal to a positive threshold ρ.We use a numerical example to illustrate the impact of the positive decision threshold ρ on the TP and FP rates.For the TP rate, we assume the input is a long pulse of duration d = 3000 and generate 100 realisations of X * (t) and Z X * (t).We use the data to compute 100 values of L(d) and L(d) where L(d) and L(d) are, respectively, the log-likelihood ratio (4) and approximate log-likelihood ratio (9) at time t = d.We then compare these values against the threshold to obtain the TP rates.We vary the threshold between 10 and 80.The blue and red lines in Fig. 7 show the TP rates for L(t) and L(t) respectively.It shows that there is a wider range of thresholds that can be used to obtain a high TP rate for L(t).This is because L(t) has a lower variance in comparison.We will discuss the magenta curve in Fig. 7 in Sec.V-B.
the FP rate, we generate 100 realisations using a short pulse of duration d = 100(= d 0 ) and use the data to compute  16) to the C1-FFL.

L(d) and L(d).
For all the 100 realisations, the L(d)'s are negative and L(t) ≈ 0. Therefore for all the thresholds in 10 − 80, the FP rates are zero for both L(t) and L(t).This shows that we are able to find thresholds that give a large TF while keeping FP low.
The above numerical experiment on computing the FP rate also shows that it is not a problem to "round" the negative log-likelihood ratio in L(d) to a near zero approximation L(d).This is because, if the test is to check whether the loglikelihood ratio is above a sufficiently positive threshold, then a negative L(d) or an almost zero L(d) will lead to the same decision.

V. USING THE C1-FFL TO APPROXIMATELY COMPUTE LOG-LIKELIHOOD RATIO
We have shown in the previous section that the mean of the intermediate approximation E[ L(t)] is an accurate approximation of the mean log-likelihood ratio E[L(t)] when the input is persistent.The aim of this section is to show that we can use the C1-FFL in (1) to approximately compute E[ L(t)] in (16).

A. Relating E[ L(t)] to the C1-FFL
For the time being, we will assume d z in (1c) is zero and show that it is possible to find C1-FFL parameters such that z(t) in ( 1) is approximately equal to E[ L(t)] in (16).We will explain in Remark 2 how a non-zero d z can be handled.
We depict the calculations of ( 16) in Fig. 8.We split the computation on the RHS of ( 16) as the product of E[ Z X * (t)] and g − π(t) [φ(E[X * (t)])] + where the multiplication operation is depicted as an AND gate in the figure.Note that the two-branch structure of the computation in Fig. 8 has a direct resemblance with that of the C1-FFL in Fig. 2. We first consider the computation of the two branches separately.
We first consider the computation of E[ Z X * (t)] from E[X * (t)], which is the branch on the left in Fig. 8.By using the volume scaling assumption (Sec.III), we equate molecular count Z X * (t) with concentration z X * (t), and similarly, X * (t) with x * (t).We propose to compute z X * (t) from x * (t) by using: Note that the above approximation is obtained from assuming that the TF-gene reactions (2c) and (2d) are close to steady state.By using the parameters in the numerical example in Section IV-B1, we have plotted the two sides of ( 18) in Fig. 9a when the input is persistent with a duration d of 200.Note that there are two transients of z X * (t) in the figure, one after time 0 when the input turns ON and the other at time d when the input turns OFF.We can see that, other than these two transients, the two sides of ( 18) are almost equal.We will show later on these two transients have little effect on the accuracy of the overall computation.Note that the RHS of ( 18) has the form of a Hill function and we can identify it with which is depicted by a cloud in Fig. 8.We first argue that, for most of the admissible choices of d 0 , there must be a time delay element in the cloud.This implies that there must be some chemical reactions in the cloud in order to create this time delay.To understand which d 0 is admissible, we recall that we assume in Section IV-A that the pathway ( 2) is almost at steady state by the time d 0 .As an illustration of this assumption, consider Fig. 9a which shows the time profile of z X * (t), the vertical dashed line shows the time (which we will denote as t ss ) by which z X * (t) is sufficiently close to steady state, the assumption means an admissible d 0 must be greater than or equal to t ss .
In Appendix D, we show that there must be a delay in the cloud if d 0 is strictly greater than t ss .Intuitively, a delay is needed because π(t) is zero in [0, d 0 ) but X * almost reaches steady state before d 0 .Given that there must be a delay element in the cloud in Fig. 8, we can achieve that by inserting a transcription node (e.g.Node Y in Fig. 2) in the cloud.With this insertion, we can identify the branch on the right in Fig. 8 with the indirect branch in the C1-FFL in Fig. 2.
The next step is to show that we can find Hill functions H xy () and H yz () in (1) which will enable us to compute g − π(t) [φ(x * (t))] + .Since the inducer-TF pathway is fast, we can approximate x * (t) by its steady state value x * ,ss for t < d.We can use (1b) to show that Our aim is to achieve the approximation: We will turn this into two requirements.First, in order to imitate π(t), we require H yz (y(t)) to have a sharp transition from a low value to a high value around d 0 .Note that this requirement implies that y(t) should not have reached steady state by d 0 and can be achieved by suitable choice of d y .Second, when H yz (y(t)) reaches steady state, which happens some time after time d 0 , its amplitude is determined by g − [φ(x * )] + .We can consider these two requirements separately because the first requirement is concerned with a time event around time d 0 while the second requirement is concerned with the steady state amplitude which is reached at a time later than d 0 .
We first consider how H yz (y(t)) can be used to realised a transition around time d 0 .Let us recall that H yz (y(t)) has the form hyzy(t) nyz K nyz yz +y(t) nyz .Given that y(t) is an increasing function of t, if we choose K yz to be y(d 0 ), then we have y(t) < K yz for t < d 0 and y(t) > K yz for t > d 0 .Hence, if n yz is sufficiently large, then H yz (y(t)) will have a sharp transition around d 0 .
The argument in the last paragraph works for a particular choice of x * ,ss , which corresponds to a particular choice of input signal amplitude.We will now explain how H yz (y(t)) can be used to realised a transition at time d 0 for a range of input amplitudes.Since H xy (x * ,ss ) is a Hill function, then for sufficiently large x * ,ss , the value of H xy (x * ,ss ) does not change much because H xy (x * ,ss ) saturates.As a result, there is a range of x * ,ss such that y(t) does not change a lot.This means if we choose K yz to be the y(d 0 ) corresponding to a particular x * ,ss which saturates H xy (x * ,ss ), we can a achieve sharp transition around time d 0 for a range x * ,ss .
We will now give an illustration on the first requirement.Fig. 9b plots the two sides of (20) for a particular input amplitude.We see that we can use H yz to approximate g − π(t) [φ(x * ,ss )] + well except near the falling edge of g − π(t) [φ(x * ,ss )] + .We can now consider the two branches in Fig. 8 together.We now compare the whole RHS of ( 16), which is and the RHS of (1c), which is H xz (x * (t)) times H yz (y(t)).Fig. 9c compares these two expressions.An interesting observation is that the mismatch in some time intervals in Fig. 9a and 9b are cancelled out when the multiplication is made.
We now consider the second requirement which requires the steady state value of H yz (y(t)) to be g − [φ(x * )] + .According to (19), the steady state amplitude of y(t) is Hxy(x * ,ss) dy , therefore the requirement is equivalent to finding H xy and H yz such that: holds for a large range of x * ,ss .Let χ = 1 x * ,ss .By using the expressions of H xy () and H yz () in (1) as well as φ() in (10), we can rewrite (21) in terms of χ as: where η = Kyzdy hxy nyz .If we consider the expression on the RHS of ( 22) as a function of χ, it has a linearly decreasing segment followed by a constant segment at zero.Recall that we have chosen n yz to be large earlier, and if we furthermore choose n xy is to be 1, then the LHS may be approximated by three segments: a constant for small χ, linearly decreasing for intermediate χ and zero for large χ.This shows that it is possible to find Hill functions H xy and H yz such that the second requirement holds.We will demonstrate this with a numerical example.

B. Numerical example
This numerical example uses the first set of parameter values as in Sec.IV-B1.Our aim is to show that we can choose a set of C1-FFL parameters so that z(t) ≈ L(t) for a range of input amplitude α.We do this by fixing H xz () as N g+x * (t) g+x * (t)+g− , and by fitting the C1-FFL parameters in H xy () and H yz () in ( 1) by nonlinear optimisation.The data for fitting is obtained from varying α from 3.125 to 62.5.For each α, we compute E[ L(t)] using (16)  Fig. 10a compares E[L(t)] and the C1-FFL output z(t) for three different values of α: 7.5, 30 and 75.It can be seen that they match very well.Next, we compare the value of z(t) and L(t) at time t = 800 for α ∈ [1.9, 112].Fig. 10b shows that the match is good for a large range of α.
Although we use optimisation to obtain the parameters of the Hill functions H xy () and H yz ().We find that their values are compatible to our intuitive argument in Sec.V-A.First, we say that K yz should be chosen close to the value of y(d 0 ) for a range of input amplitudes.For input amplitudes of 5, 15 and 480 (which correspond to a mean probability of 0.4, 0.6 and 0.8 that the promoter Z X * is bound), the values of y(d 0 ) are respectively, 229.2, 247.3 and 252.9.These values are pretty close to K yz = 255.4obtained from using optimization.Second, the steady state value of y(t) is around 259 for a range of input amplitudes and y(d 0 ) has not yet reached the steady state.Third, we argue that H yz ( Hxy(x * ,ss) dy ) and g − [φ(x * ,ss )] + should be approximately equal, i.e. (21) or (22).Fig. 10c plots the two sides of ( 22) and shows that they match well except near the turning point.Fourth, the optimised value of n xy is 1.65 which is similar to the predicted value of 1.
We next study the behaviour of the fitted C1-FFL when s(t) in ( 1) is the stochastic persistent input used in Sec.IV-B1.We apply the tau-leaping simulation algorithm in [26] to simulate the stochastic behaviour of the C1-FFL from its ODE model.We perform 100 rounds of simulation to obtain the mean of the C1-FFL output z(t).We compare the mean of z(t) against the mean log-likelihood computed in Sec.IV-B1 in Fig. 10d and we can see that they match fairly well.Next, we use z(t) to study the TP (true positive) and FP (false positive) rates in the same way as in Sec.IV-C1.Both persistent and transient inputs are used in this study, and 100 simulation runs are performed for each type of input.The magenta line in Fig. 7 (Page 8) shows how the TP rates vary with the decision threshold.It can be seen that the TP curve for the C1-FFL is in between those for L(t) and L(t).We find that we can again find a range of positive thresholds such that the TP rate is high and FP rate is zero.
Remark 2: We have shown that on the condition that d z in (1c) is zero, then we can find parameters of the C1-FFL such that z(t) in ( 1) is approximately equal to E[ L(t)] in (16).Let us now assume d z is non-zero and we add the term −d z E[ L(t)] to the RHS of ( 16), then we can use the same C1-FFL parameters as before to make z(t) in (1) to be approximately equal the new E[ L(t)].
Remark 3: The derivations in Secs.III-V assume that the reactions take place in a volume V of 1.If the volume is not 1, ( 16) remains the same.Note that both g − and the term inside [ ] − in ( 16) are independent of V .If we divide both sides of ( 16) by V , we have Thus, the argument in Sec.V-A is still valid if we interpret the output of the C1-FFL as log-likelihood ratio per volume.
Remark 4: Our interpretation of the C1-FFL as a statistical persistence detector shows an interesting signal processing architecture involving parallel processing.The C1-FFL has two arms.The short arm produces the signals X * (t) and Z X * (t) which contain information on whether the inducer signal is persistent or not.The longer arm then makes use of the signals X * (t) and Z X * (t) to approximately compute, in a parallel manner, the approximate log-likelihood ratio.It is also interesting to see that the gene Z, with its two promoter sites, has the dual roles of generating the signal to be processed as well as processing the signal.

C. Discussion
Biochemical circuits are known to be impacted by both intrinsic and extrinsic noise [27].An open problem is to design a C1-FFL, which is based on the mass-reaction kinetics, for persistence detection using the probabilistic framework studied in this paper.In particular, we envisage using this design problem to study the impact intrinsic noise (which comes from reaction kinetics) and extrinsic noise (which comes from input signal properties such as how the signal is generated, amplitude, duration etc.) on the detection and false alarm rates.
Another open problem is to study how our work can be used in conjunction with the sequential probability ratio test (SPRT) [28].The SPRT uses the log-likelihood ratio L(t), such as that defined in (4) which uses data up till time t, with two thresholds θ 1 and θ 0 where θ 1 > θ 0 .At any point in time, if L(t) > θ 1 or L(t) < θ 0 is true, then a decision is reached (deciding for, respectively, H 1 or H 1 ); otherwise a decision has not been reached and data collection will continue.It is suggested in [29] that the SPRT is useful for biological detection because theoretically SPRT has a shorter decision  ).Note the horizontal axis is χ = 1 x * ,ss (d) Comparing the mean C1-FFL output from tau-leaping against mean L(t) from SSA.
time.We can use our approximate log-like ratio L(t) in the spirit of SPRT.In our case, only one threshold θ 1 (> 0) will be used.At any point in time, if L(t) > θ 1 then a decision is reached for H 1 ; otherwise a decision has not been reached.(Note the similarity to the study in Sec.IV-C1.)This set up will be useful for practical applications where an action is needed if H 1 is true but otherwise no actions are required.We envisage that this set up will have a shorter decision time as in SPRT.It will be interesting to study how the design parameters mentioned in the last paragraph will influence the decision time and accuracy.This study is important for practical implementation.

VI. CONCLUSIONS
This paper shows that the C1-FFL gene circuits can be used to approximately compute the log-likelihood ratio for statistical detection of persistent signals.In order to derive this result, we use new methods to show how log-likelihood ratio can be approximately computed and to show how this approximation can be mapped to the parameters of a C1-FFL.An open problem is to study the impact of intrinsic and extrinsic noise on the performance of the C1-FFL based statistical detector.Another open problem is to consider using the approximate log-likelihood ratio computation together with sequential probability ratio test to realise fast persistent detection.Finally, note that the methodology in this paper can be used to derive biochemical circuits that can be used to solve other statistical signal processing problems.APPENDIX A PROOF OF (4) Recalling that Z X * (t) is the history of Z X * (t) in the time interval [0, t].In order to derive (4), we consider the history Z X * (t + ∆t) as a concatenation of Z X * (t) and Z X * (t) in the time interval (t, t + ∆t].We assume that ∆t is chosen small enough so that no more than one reaction can take place in (t, t + ∆t].Given this assumption and right continuity of continuous-time Markov Chains, we can use Z X * (t + ∆t) to denote the history of Z X * (t) in (t, t + ∆t]. Consider the likelihood of observing the history Z X * (t + ∆t) given hypothesis H i : where we have expanded Z X * (t + ∆t) in ( 23) using concatenation.
By using (25) in the definition of log-likelihood ratio, we can show that: The condition probability P[ Z X * (t+∆t)|H i , Z X * (t)] is the prediction of the number of Z X * molecules at time t + ∆t based on its history up till time t.This conditional probability can be obtained by solving a Bayesian filtering problem over a continuous-time Markov chain which describes the dynamics of the chemical reactions in (2) and those that produce S [21].We considered how this conditional probability could be evaluated in our earlier work [30].The key result in [30] says that P[ Z X * (t + ∆t)|H i , Z X * (t)] can be expressed in terms of the predicted rate of the chemical reactions that Z X * are involved in.By using [30], [21], we have: where δ a,b is the Kronecker delta which is 1 when a equals to b and zero otherwise, and is the expected number of X * molecules at time t given Hypothesis i and the history Z X * (t).
Note that P[ Z X * (t + ∆t)|H i , Z X * (t)] in ( 27) is a sum of three terms with multipliers δ Z X * (t+∆t), Z X * (t)+1 , δ Z X * (t+∆t), Z X * (t)−1 and δ Z X * (t+∆t), Z X * (t) .Since these multipliers are mutually exclusive, we have: where we have used the approximation log(1 + f ∆t) ≈ f ∆t and have ignored terms of order (∆t) 2 or higher to obtain (28).Note also that the above derivation assumes that J1(t) J0(t) is strictly positive so its logarithm is well defined; this can be achieved by proper choice of the hypotheses.
By substituting ( 28) into ( 26), we have after some manipulations and after taking the limit ∆t → 0: In order to obtain (4), we use the following reasonings.First, the term lim ∆t→0 ∆t is a Dirac delta at the time instant that an X * molecule binds with Z to form a Z X * .Since the binding instants are also the times at which Z X * (t) jumps by +1, we can identify this term with δ Z X * (t+∆t), Z X * (t) = 1 with probability one.This allows us to drop δ Z X * (t+∆t), Z X * (t) .Hence (4).We remark that a more general framework of deriving log-likelihood ratio and logposteriori probability in the reaction-diffusion master equation framework can be found in [31], [30], [32].
APPENDIX B DERIVATION OF (9) The aim of this section is to show that the log-likelihood ratio computation in (4) can be approximated by the intermediate approximation in (9) for persistent signals in the time interval [0, min(d, d 1 )].The derivation is based on the assumptions stated in Sec.IV-A.We have divided the derivation into three steps.
(Step 1) The aim of this step is to approximate E[X * (t)| Z X * (t), H i ] in (4).The computation of this expectation requires the solution of a Bayesian filtering problem which is computationally intensive.We will argue that if g + and g − are small, then we can use the approximation: where i = 0, 1, and X 0 and X 1 are defined in (11).In other words, we approximate E[X * (t)| Z X * (t), H i ] by a rectangular pulse.
In order to argue for (30), we start by deriving the solution of the Bayesian filtering problem of using the history Z X * (t) to determine the posteriori probability distribution pt (s, x * ) P[S(t) = x, X * (t) = x * | Z X * (t), H i ].The typical temporal behaviour of pt (s, x * ) is depicted in Fig. 11.We can see that at the time instants where Z binds or Z X * unbinds, the probability pt (s, x * ) has a discrete jump; at other times, i.e. between two consecutive jumps, the probability pt (s, x * ) varies continuously.By using [21, Eq. ( 21)], we can show that the time evolution of the continuously varying part of pt (s, x * ) is governed by the following ODE: If g + is small, then the last term in (31) can be viewed as a perturbation to the ODE: This ODE is in fact the chemical master equation which governs the reactions ( 6), (2) with g + = g − = 0.For small g + , we have pt (s, x * ) ≈ p t (s, x * ).In addition, if both g + and g − are small, then the time between two jumps in pt (s, x * ) is long; this further implies that pt (s, x * ) is mostly around steady state, e.g.see Fig. 11 in the time interval [320,386].Overall, this means that we can approximate E[X * (t)| Z X * (t), H i ] by using the steady state mean number of X * molecules assuming g + = g − = 0.When g + = g − = 0, the steady state mean number of X * molecules is approximately given by the expression in (11).Hence (30).After using the approximation in Step 1, (4) becomes: where π(t) is defined in (14).
We now demonstrate the accuracy of ( 30) and (34) using numerical examples.We use the same parameter values as in Sec.IV-B1.Fig. 12a plots the two sides of (30) for the two hypotheses.It can seen that, other than the transients, the approximation is fairly accurate.Next, we use Fig. 12b to demonstrate the accuracy of (34).By using one realisation of Z X * (t), we calculate the exact log-likelihood ratio L(t) (4) and the approximation L(t) (34).The figure shows that the approximation is accurate.In addition, the figure also shows the RMS error between L(t) − L(t) for 100 realisations of Z X * (t).It can be seem that the RMS error is small.In order to show that the approximation holds for different parameter settings, we use a different value S Note that the RHS of ( 34) is the sum of two terms that do not depend on L(t).We can therefore consider the contribution of each term to L(t) separately.
First, we consider the contribution of the first term on the RHS of (34) to L(t), which we will call L 1 (t): By using a method similar to that in [5], which is based on the renewal theorem [33], we can show that (35) can be approximated by: Next, we consider the contribution of the second term on the RHS of (34) to L(t), which we will call L 2 (t): By integrating the above equation, we have L 2 (t) = 0 for t ∈ [0, d 0 ); for t ∈ [d 0 , min(d, d 1 )], we have: Since the reaction pathway is in steady state in the time interval [d 0 , min(d, d 1 )], we can replace the time average in (38) by its ensemble average.In this part, we will overload the symbol Z X * to use it to refer to the random variable of the number of Z X * molecules at steady state.This should not cause any confusion because the meaning should be clear from the context.In addition, we will overload the symbol X * in

Fig. 3 :
Fig. 3: The reaction pathway for the detection problem.

Fig. 5 :
Fig. 5: (a) Sample persistent S(t) (top plot) and transient S(t) (bottom plot).(b) The log-likelihood ratio for a long signal (or persistent signal) and a short signal.(c) The top and bottom plots show log J1(t) J0(t) [d 0 , min(d, d 1 )].Intuitively, this assumption allow us to use the steady state statistics in the time interval [d 0 , min(d, d 1 )] to replace d Z X * (t) dt + and g + (N − Z X * (t)) in (4) by, respectively, g − Z X * (t) and g − Z X * (t) X * (t) .

Fig. 6 :
Fig. 6: Numerical results for Sec.IV-B1.(a) Comparing one realisation of L(t) and L(t) for a persistent signal.Also shows a realisation of L(t) for a transient signal.(b) Comparing the mean of L(t) from SSA, mean of L(t) from SSA, mean of L(t) by (16).Also shows RMS error between L(t) and L(t) from SSA. (c) Same type of comparison as (b) but for different values of α and d.
The utilities and costs of the four combinations of environmental conditions and detection outcomes.

Fig. 7 :Fig. 8 :
Fig. 7: The impact of the decision threshold on the TP rates.
and use them to fit the parameters of the C1-FFL.The fitted values of the C1-FFL parameters are: k xy = 9.66, n xy = 1.65,K xy = 4.36, d y = 0.037, h yz = 0.135, n yz = 50.0and K yz = 255.4.

Fig. 12 :
Fig. 12: (a) Demonstrating the accuracy of (30).Top plot: H 0 ; Bottom plot: H 1 .(b) and (c) Comparing exact log-likelihood ratio L(t) (4) and the approximation L(t) (34).For (b), S ON,ref ref d y and d z are reaction rate constants; h xy , n xy , K xy , h xz , n xz , K xz , h yz , n yz and K yz are coefficients for Hill functions H xy (x * (t)), H xz (x * (t)) and H yz (y(t)).Lastly, x(t) + x * (t) is the constant M .The multiplication of H xz and H yz The first two properties are concerned with transient signals, which are those input signals whose duration is no longer than d 0 or whose amplitude α is small.The intermediate approximation says that transient signals give a small L(t).On the other hand, persistence signals have a duration longer than d 0 and have a sufficiently large amplitude α.For persistent signals, the intermediate approximation L(t) is ).Three important properties for L(t) are: 1) If d < d 0 , then for all t, we have L(t) and E[ L(t)] are zero or small.This is due to π(t), which is zero outside of [d 0 , d 1 ), and the fact that X * (t) is likely to be small for t ≥ d 0 .2) If the amplitude α is lower than a threshold, then for all t, we have L(t) is zero or small and E[ L(t)] is zero.We will explain this for E[ L(t)].Since E[X * (t)] is an increasing function of α, this means a small α will give a small E[X * (t)].If E[X * (t)] is less than X1−X0 log X 1 X 0 for all t, then [φ(E[X * (t)])] + on the RHS of (16) is zero and this implies E[ L(t)] is zero for all t.The explanation for L(t) is similar.3) If d is longer than d 0 and α is sufficiently large, then for 0 ≤ t < min{d, d 1 } we have L(t) ≈ L(t) where L(t) is given in (4).
ON,ref ref and plot the result in Fig. 12c.(Step 2) The aim of this step is to replace d Z X * (t) dt + and Z X * (t) in (34) by alternative expressions.Since π(t) in (34) is zero for t < d 0 , we only have to consider input signals whose duration d ≥ d 0 .Recall that we assume in Sec.IV-A that the pathway (2) reaches steady state by d 0 if the input has a duration of at least d 0 .This means that the probability distributions of X * and Z X * are in steady state in the time interval [d 0 , d]