Using biochemical circuits to approximately compute log-likelihood ratio for detecting persistent signals

TL;DR

This paper explores how biochemical circuits can be designed to approximate the computation of log-likelihood ratios for detecting persistent signals, using stochastic modeling and gene circuit implementations.

Contribution

It introduces a method to approximate log-likelihood ratio computation for persistent signal detection using biochemical circuits, based on stochastic modeling and detection theory.

Findings

Biochemical circuits can approximate log-likelihood ratios for signal detection.

Time-scale separation allows simplified computation of detection statistics.

Coherent feedforward gene circuits can implement the detection algorithm.

Abstract

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.

Equations

$$\frac{d{x}_{\ast }(t)}{dt}=$$

$$\frac{dy(t)}{dt}=$$

$$\frac{dz(t)}{dt}=$$

$$L(t)=\log \left(\frac{\mathrm{P}[{\tilde{\mathcal{Z}}}_{X_{\ast }}(t)|{\mathcal{H}}_1]}{\mathrm{P}[{\tilde{\mathcal{Z}}}_{X_{\ast }}(t)|{\mathcal{H}}_0]}\right)$$

$$\frac{dL(t)}{dt}=$$

$$g_{+}(N-{\tilde{Z}}_{X_{\ast }}(t))(J_1(t)-J_0(t))$$

$$J_i(t)=$$

$$\mathrm{E}[X_{\ast }(t)|{\tilde{\mathcal{Z}}}_{X_{\ast }}(t),{\mathcal{H}}_i]$$

$$\frac{d{\hat{p}}_t(s,x_{\ast })}{dt}$$

$$=$$

$$-g_{+}(x_{\ast }-\mathrm{E}[X_{\ast }(t)|{\tilde{\mathcal{Z}}}_{X_{\ast }}(t),{\mathcal{H}}_i])(N-{\tilde{\mathcal{Z}}}_{X_{\ast }}(t))$$

$$[{\mathcal{L}}_i{\hat{p}}_t](s,x_{\ast })$$

$$=$$

$$f_{-}(s+1)\hat{p}(s+1,x_{\ast })-f_{-}s{\hat{p}}_t(s,x_{\ast })+$$

$$k_{+}s(M-x_{\ast }+1)\hat{p}(s,x_{\ast }-1)-$$

$$k_{+}s(M-x_{\ast }){\hat{p}}_t(s,x_{\ast })+$$

$$k_{-}(x_{\ast }+1)\hat{p}(s,x_{\ast }+1)-k_{-}{x}_{\ast }{\hat{p}}_t(s,x_{\ast })$$

$$\frac{d{p}_t(s,x_{\ast })}{dt}=[{\mathcal{L}}_i{p}_t](s,x_{\ast })$$

$$\frac{d\tilde{L}(t)}{dt}\approx$$

$$g_{+}(N-{\tilde{Z}}_{X_{\ast }}(t))(X_1-X_0)\pi (t)$$

$$\frac{d{L}_1(t)}{dt}=$$

$$\frac{d{L}_1(t)}{dt}\approx$$

$$\frac{d{L}_2(t)}{dt}=$$

$$L_2(t)=$$

$$L_2(t)\approx$$

$$g_{+}\mathrm{E}[(N-{\tilde{Z}}_{X_{\ast }})X_{\ast }]=$$

$$\frac{1}{\mathrm{E}[Q]}\approx$$

$$I(\frac{1}{q})=$$

$$\mathrm{E}[(N-{\tilde{Z}}_{X_{\ast }})X_{\ast }]\approx$$

$$\mathrm{E}[(N-{\tilde{Z}}_{X_{\ast }})X_{\ast }]=$$

$$\mathrm{E}[(N-{\tilde{Z}}_{X_{\ast }})I(\frac{1}{X_{\ast }})]\approx$$

$$g_{+}\mathrm{E}[(N-{\tilde{Z}}_{X_{\ast }})]\approx$$

$$L_2(t)\approx$$

$$\frac{d{L}_2(t)}{dt}\approx$$

Full text

Keywords:

Statistical signal processing, signal detection, molecular computing, analog computation, biochemical circuits.

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 [Teo:2019bk], concentration regulation [Pasotti:2019fi] and to counter the effect of loading [McBride:2019is]; signal processing theory has been used to design molecular circuits that can be used for communications [Kuscu:2019ey, Chou:2019gf, Awan:2018ig] and filtering [Zechner:2016fb]. Secondly, one can view the information processing carried out by biomolecular circuits as analog computation [Sarpeshkar:2014dg, Teo:2015kva]. This view point has helped engineers to design molecular circuits that can perform logarithmic sensing [Daniel:2013ke], parity check [Marcone:2018kp] and integral control [Briat:2016ds]. 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 [Mangan:2003ia]. In the synthetic world, cell-based therapy [Mimee:2015jo] can make use of persistence detection on biomarkers. In this paper, we will solve the persistence detection problem by using statistical detection theory [Kay_v2] 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 [ShenOrr:2002jo] 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) [Milo:2002cg, ShenOrr:2002jo]. This means that the C1-FFL carries out important functions in cells. The authors in [ShenOrr:2002jo] 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 [ShenOrr:2002jo] 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 [Libby:2007fw]. 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 [Gardiner]. 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 log-likelihood 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 log-likelihood 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 [Chou:2018jh]. In comparison to [Chou:2018jh], this paper makes two different assumptions: (i) This paper considers an inducer-TF-gene pathway but [Chou:2018jh] considered only an inducer-TF pathway; (ii) This paper assumes that the inducer signal is stochastic while [Chou:2018jh] 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-TF-gene pathway. Second, we need to derive a new method to approximately compute the log-likelihood ratio (Sec. LABEL:sec:approx_llr). In particular, this paper needs to approximate the solution to a Bayesian filtering problem [Bronstein:2018eh] but this is not required in [Chou:2018jh] 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. LABEL:sec:c1ffl). These three aspects are the new elements of this paper in comparison to [Chou:2018jh]. Furthermore, in comparison with our earlier conference paper [Chou:2019ur], this paper explains why a C1-FFL can be used to approximately compute the approximate log-likelihood ratio (Sec. LABEL:sec:c1ffl) and provide full derivation on the computation of exact and approximate log-likelihood ratios (Appendices LABEL:app:sol:dp and LABEL:app:ia).

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. LABEL:sec:approx_llr, we present a method to approximately compute the log-likelihood ratio and use this approximation in Sec. LABEL:sec:c1ffl to show that the C1-FFL can be used to approximately compute the the log-likelihood ratio when it is positive. Finally, Sec. LABEL:sec:final 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 \ceeS_y are inducers. Both X and Y are TFs, which are expressed by their corresponding gene. The inducer S (resp. \ceeS_y) turns the inactive form X (Y) into its active form \ceeX_* (\ceeY_). The activation of gene Z requires the binding of both \ceeX_ and \ceeY_* 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 with their corresponding time signals in Fig. 1, e.g. $x_{*}(t)$ is the concentration of \ceeX_* at time $t$ and so on. In this paper, we will assume that the inducer \ceeS_y is always present and its concentration is always above the threshold needed to activate Y. Furthermore, we assume the activation of Y by \ceeS_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, [Mangan:2003ia] presents an ODE model for the C1-FFL, as follows:

[TABLE]

where $k_{+}$, $k_{-}$, $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}$ 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. LABEL:sec:approx_llr. 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, \ceeX_* 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 \ceeX_*, 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$).

III-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 \ceeX and its corresponding active form \ceeX_, as well as inactive $\widetilde{\rm Z}$ and the complex $\widetilde{\rm Z}_{X_{*}}$ which is formed by the binding of \ceeX_ to $\widetilde{\rm Z}$. These five species take part in the following four chemical reactions:

[TABLE]

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 LABEL:re:v. With this assumption, note that $k_{+}$ and $k_{-}$ in (2a) and (2b) are equal to those in (1a).

In terms of molecular biology, \ceeS is an inducer and \ceeX is a TF. In Reaction (2a), the species \ceeS activates \ceeX to produce \ceeX_*. Reaction (2b) is a deactivation reaction. The reactions (2a) and (2b) are depicted in both Figs. 2 and 3.

The species $\widetilde{\rm Z}$ is a gene. In fact, $\widetilde{\rm Z}$ in Fig. 3 is the same as \ceeZ 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 \ceeZ. 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 $\widetilde{\rm Z}$ to denote the gene and use \ceeZ to denote the protein expressed by $\widetilde{\rm Z}$. In Reaction (2c), an active \ceeX_* binds with the promoter of $\widetilde{\rm Z}$ to produce the complex $\widetilde{\rm Z}_{X_{*}}$. Lastly, Reaction (LABEL:cr:off2) is an unbinding reaction.

Let $S(t)$, $X(t)$, $X_{*}(t)$, $\widetilde{Z}(t)$ and $\widetilde{Z}_{X_{*}}(t)$ denote, respectively, the number of \ceeS, \ceeX, \ceeX_*, $\widetilde{\cee Z}$ and $\widetilde{\rm 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. $\widetilde{Z}(t)+\widetilde{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 $\widetilde{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 [Gardiner] which is a specific type of continuous-time Markov chain. We further assume that the upstream species do not react with \ceeX, \ceeX_, $\widetilde{\cee Z}$ and $\widetilde{\rm Z}_{X_{*}}$. Intuitively, this means we can predict the behaviour of \ceeX, \ceeX_, $\widetilde{\cee Z}$ and $\widetilde{\rm Z}_{X_{*}}$ from that of \ceeS.

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 $\widetilde{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 $\widetilde{Z}_{X_{*}}(\tau)$ for all $\tau\in[0,t]$; in other words, the data are continuous in time and are the history of the counts of $\widetilde{\rm Z}_{X_{*}}$ up to time $t$ inclusively. We will use ${\widetilde{\cal Z}}_{X_{*}}(t)$ to denote the continuous-time history of $\widetilde{Z}_{X_{*}}(t)$ up to time $t$ inclusively.

We now specify the hypotheses ${\cal H}_{i}$ $(i=0,1)$ for the detection problem. Later on, we will identify ${\cal H}_{0}$ and ${\cal 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 ${\cal H}_{0}$ and ${\cal 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 ${\cal H}_{0}$ or ${\cal 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 \ceeS reacts with the molecules in the reaction pathway in Fig. 3, the downstream signal $\widetilde{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 $\widetilde{Z}_{X_{*}}(t)$ are noisy.

III-B Solution to the detection problem

The aim of the detection problem is to decide which hypothesis ${\cal H}_{i}$ $(i=0,1)$ is more likely to have generated the observed history ${\widetilde{\cal Z}}_{X_{*}}(t)$.

We assume that the two hypotheses are equally likely, therefore the log-ratio of posteriori probabilities $\frac{{\rm P}[{\cal H}_{1}|{\widetilde{\cal Z}}_{X_{*}}(t)]}{{\rm P}[{\cal H}_{0}|{\widetilde{\cal Z}}_{X_{*}}(t)]}$ is equal to the log-likelihood ratio $L(t)$:

[TABLE]

where ${\rm P}[{\widetilde{\cal Z}}_{X_{*}}(t)|{\cal H}_{i}]$ is the conditional probability of observing the history ${\widetilde{\cal Z}}_{X_{*}}(t)$ given hypothesis ${\cal H}_{i}$.

In Appendix LABEL:app:sol:dp, we show that the time evolution of $L(t)$ is given by the following ODE:

[TABLE]

where $[w]_{+}=\max(w,0)$, ${\rm E}[\;]$ denotes the expectation and ${\rm E}[X_{*}(t)|{\widetilde{\cal Z}}_{X_{*}}(t),{\cal H}_{i}]$ is the conditional expectation of $X_{*}(t)$ given the history and ${\cal H}_{i}$. Note that in deriving (2d), 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\left(\frac{J_{1}(t)}{J_{0}(t)}\right)$ is well defined.

Since $\widetilde{Z}_{X_{*}}(t)$ is a piecewise constant function counting the number of $\widetilde{\rm Z}_{X_{*}}$ molecules, its derivative is a sequence of Dirac deltas at the time instants that $\widetilde{\rm Z}_{X_{*}}$ forms or unbinds. Note that the Dirac deltas corresponding to the formation of $\widetilde{\rm Z}_{X_{*}}$ carries a positive sign and the $[\;]_{+}$ operator keeps only these. Fig. 4 shows an example $\widetilde{Z}_{X_{*}}(t)$ and its corresponding $\left[\frac{d\widetilde{Z}_{X_{*}}(t)}{dt}\right]_{+}$.

At the time instant $t$ that $\widetilde{\rm 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 (2d) 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 (2d) and to explain what information is important for persistence detection. This example will also provide some intuition on how we will approximately compute the log-likelihood in Sec. LABEL:sec:approx_llr.

III-B1 Numerical example

The aim of this example is to show that we can use the log-likelihood ratio computed by (2d) 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:

[TABLE]

(Step 1) The aim of this step is to approximate ${\rm E}[X_{*}(t)|{\widetilde{\cal Z}}_{X_{*}}(t),{\cal H}_{i}]$ in (2d). 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:

[TABLE]

where $i=0,1$, and $X_{0}$ and $X_{1}$ are defined in (LABEL:eq:ia:Xi). In other words, we approximate ${\rm E}[X_{*}(t)|{\widetilde{\cal Z}}_{X_{*}}(t),{\cal H}_{i}]$ by a rectangular pulse.

In order to argue for (2fat), we start by deriving the solution of the Bayesian filtering problem of using the history ${\widetilde{\cal Z}}_{X_{*}}(t)$ to determine the posteriori probability distribution $\hat{p}_{t}(s,x_{*})\triangleq{\rm P}[S(t)=x,X_{*}(t)=x_{*}|{\widetilde{\cal Z}}_{X_{*}}(t),{\cal H}_{i}]$. The typical temporal behaviour of $\hat{p}_{t}(s,x_{*})$ is depicted in Fig. 11. We can see that at the time instants where $\widetilde{\rm Z}$ binds or $\widetilde{\rm Z}_{X_{*}}$ unbinds, the probability $\hat{p}_{t}(s,x_{*})$ has a discrete jump; at other times, i.e. between two consecutive jumps, the probability $\hat{p}_{t}(s,x_{*})$ varies continuously. By using [Bronstein:2018eh, Eq. (21)], we can show that the time evolution of the continuously varying part of $\hat{p}_{t}(s,x_{*})$ is governed by the following ODE:

[TABLE]

where

[TABLE]

If $g_{+}$ is small, then the last term in (2fau) can be viewed as a perturbation to the ODE:

[TABLE]

This ODE is in fact the chemical master equation which governs the reactions (2f), (2) with $g_{+}=g_{-}=0$. For small $g_{+}$, we have $\hat{p}_{t}(s,x_{*})\approx p_{t}(s,x_{*})$. In addition, if both $g_{+}$ and $g_{-}$ are small, then the time between two jumps in $\hat{p}_{t}(s,x_{*})$ is long; this further implies that $\hat{p}_{t}(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 ${\rm E}[X_{*}(t)|{\widetilde{\cal Z}}_{X_{*}}(t),{\cal H}_{i}]$ by using the steady state mean number of \ceeX_* molecules assuming $g_{+}=g_{-}=0$. When $g_{+}=g_{-}=0$, the steady state mean number of \ceeX_* molecules is approximately given by the expression in (LABEL:eq:ia:Xi). Hence (2fat).

After using the approximation in Step 1, (2d) becomes:

[TABLE]

where $\pi(t)$ is defined in (LABEL:eq:pi).

We now demonstrate the accuracy of (2fat) and (2fax) using numerical examples. We use the same parameter values as in Sec. LABEL:sec:int_approx:numer. Fig. 12(a) plots the two sides of (2fat) for the two hypotheses. It can seen that, other than the transients, the approximation is fairly accurate. Next, we use Fig. 12(b) to demonstrate the accuracy of (2fax). By using one realisation of $\widetilde{Z}_{X_{*}}(t)$, we calculate the exact log-likelihood ratio $L(t)$ (2d) and the approximation $\tilde{L}(t)$ (2fax). The figure shows that the approximation is accurate. In addition, the figure also shows the RMS error between $L(t)-\tilde{L}(t)$ for 100 realisations of $\widetilde{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_{\rm ref}^{\rm ON,ref}$ and plot the result in Fig. 12(c).

(Step 2) The aim of this step is to replace $\left[\frac{d\widetilde{Z}_{X_{*}}(t)}{dt}\right]_{+}$ and $\widetilde{Z}_{X_{*}}(t)$ in (2fax) by alternative expressions. Since $\pi(t)$ in (2fax) is zero for $t<d_{0}$, we only have to consider input signals whose duration $d\geq d_{0}$. Recall that we assume in Sec. LABEL:sec:assumptions 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 \ceeX_* and $\widetilde{Z}_{X_{*}}$ are in steady state in the time interval $[d_{0},d]$

Note that the RHS of (2fax) 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 (2fax) to $L(t)$, which we will call $L_{1}(t)$:

[TABLE]

By using a method similar to that in [Chou:2019gf], which is based on the renewal theorem [Grimmett], we can show that (2fay) can be approximated by:

[TABLE]

Next, we consider the contribution of the second term on the RHS of (2fax) to $L(t)$, which we will call $L_{2}(t)$:

[TABLE]

By integrating the above equation, we have $L_{2}(t)=0$ for $t\in[0,d_{0})$; for $t\in[d_{0},\min(d,d_{1})]$, we have:

[TABLE]

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 (2fbb) by its ensemble average. In this part, we will overload the symbol $\widetilde{Z}_{X_{*}}$ to use it to refer to the random variable of the number of $\widetilde{\rm 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 the same way. With this overloading, the mean number of \ceeX_* and $\widetilde{\rm Z}_{X_{*}}$ molecules at steady state are denoted by ${\rm E}[X_{*}]$ and ${\rm E}[\widetilde{Z}_{X_{*}}]$ respectively. We can now rewrite (2fbb) as:

[TABLE]

In order to be able to connect to the C1-FFL, we will need to replace the expression $(N-{\rm E}[\widetilde{Z}_{X_{*}}])$ by a different expression. The derivation of this replacement expression requires a few auxiliary results.

(Auxiliary Result 1) By considering the global balance of the steady state of the reaction pathway (2), we have:

[TABLE]

(Auxiliary Result 2) If the amplitude of the input is sufficiently high, then the number of \ceeX_* molecules can be approximately modelled as a binomial distribution with $M$ trials and a success probability of $\frac{k_{+}\alpha}{k_{+}\alpha+k_{-}}$.

Consider a binomial distribution $B(Q;m,f)$ with parameters $m$ (number of trials) and $f$ (success probability), then for sufficiently large $m$ and $f$, we have

[TABLE]

where

[TABLE]

This result essentially says that the mean of the reciprocal of a binomial random variable (with $\frac{1}{0}$ excluded) is approximately equal to the reciprocal of the mean of the binomial random variable. If $f=1$ and $m\geq 1$, the binomial distribution has a single outcome with a non-zero probability so (2fbe) is exact. Intuitively, if a probability has a single modal distribution with a narrow spread, then (2fbe) holds approximately. For $f=0.1$, the relative error of using (2fbe) is 3.21% for $m=300$ and drops to 1.87% for $m=500$. In general, the approximation is better for large $m$ and $f$.

(Auxiliary Result 3) Since the inducer-TF reactions are faster than the TF-gene reactions and $M\gg N$, we can show that:

[TABLE]

We will argue that the above approximation holds by using time average to compute ${\rm E}[(N-\widetilde{Z}_{X_{*}})X_{*}]$. Let $t_{0}$, $t_{1}$, $\ldots$ be a sequence of time instants at which $\widetilde{\rm Z}_{X_{*}}$ changes its value. Since the continuous-time Markov chain associated with the chemical system is ergodic, we have:

[TABLE]

Since the TF-gene reactions are slow in comparison, $\widetilde{\rm Z}_{X_{*}}$ is a slow species while \ceeX_* is a fast species. This means the time interval $[t_{i},t_{i+1})$ (during which the count of $\widetilde{\rm Z}_{X_{*}}$ is a constant) is likely to be long compared to the time-scale of the fast species \ceeX_*. This allows us to approximate the integral on the RHS of (2fbj) by ${\rm E}[X_{*}](t_{i+1}-t_{i})$. Hence (2fbi). Note that the above argument is identical to the one used in [Cao:2005gj] to derive the slow-scale tau-leaping simulation algorithm.

(Auxiliary Result 4) By using the same argument as in Auxiliary Result 3, we can show that:

[TABLE]

We will now use the above auxiliary results and (2fbc) to derive the replacement expression. By using Auxiliary Results 1 and 3, we have:

[TABLE]

We then apply Auxiliary Result 2 to the RHS of (2fbl) to obtain:

[TABLE]

By applying Auxiliary Result 4 to the RHS of (2fbm), we have:

[TABLE]

By substituting (2fbn) into (2fbc), we have:

[TABLE]

By turning the above equation into the differential form, we have:

[TABLE]

Next, by combining (2faz) and (2fbp), we have:

[TABLE]

(Step 3) Since a set of chemical reactions can be modelled by a set of ODEs, we want to turn the ODE in (2fbq) into a form that can be implemented by a set of chemical reactions. However, (2fbq) cannot be directly implemented by chemical reactions because log-likelihood ratio can take both positive and negative values but chemical concentration is always non-negative. Although [Oishi:2011ig] has derived a chemical computation system that can have both positive and negative numbers, it requires double the number of species and reactions. As in our previous work [Chou:2019gf, Chou:2018jh], we choose to compute only the log-likelihood ratio when it is positive. We do that by applying $[\;]_{+}$ to the RHS of (2fbq); we have:

[TABLE]

We now replace $I(\frac{1}{X_{*}(t)})$ in (2fbr) by $\frac{1}{X_{*}(t)}$ to obtain:

[TABLE]

The removal of $I(\;)$ will not make much difference because the probability of having $X_{*}(t)$ equals to 0 is small when the input signal is persistent. Note that (2fbs) is the same as (LABEL:eq:Lfinal). This completes the derivation for (LABEL:eq:Lfinal).

In order to derive (LABEL:eq:Lfinal_mean), we start from (2fbs) and take expectation on both sides. If the amplitude $\alpha$ is sufficiently high, then there is a high probability that $X_{*}(t)$ is large. This means we can take the expectation operator to the inside of the $[\;]_{+}$ operator. After that we apply Auxiliary Results 2, 3 and 4 to obtain (LABEL:eq:Lfinal_mean).

Appendix C Utility maximisation

In this section, we will show that the utility maximisation problem (LABEL:eq:max_util) leads to a detection criterion based on the likelihood ratio. The proof here uses the same method as Appendix 3A in [Kay_v2] for proving the Neyman-Pearson lemma.

By using Lagrangian multiplier $\lambda\geq 0$, we rewrite (LABEL:eq:max_util) as the maximisation of:

[TABLE]

Let ${\cal O}$ be the observations that are available to the detection problem. Also let ${\cal A}_{0}$ and ${\cal A}_{1}$ be two disjoint sets which are to be used as the decision regions for the detection problem. In particular, the detection problem will decide for hypothesis ${\cal H}_{1}$ (resp. ${\cal H}_{0}$) if ${\cal O}\in{\cal A}_{1}$ (${\cal O}\in{\cal H}_{0}$). Given these definitions, we can write $P_{TP}$ and $P_{FP}$ as:

[TABLE]

By substituting (2fbv) and (2fbw) into (2fbu), we have the objective to be maximised is:

[TABLE]

In order to maximise (2fbx), we should include ${\cal O}$ in ${\cal A}_{1}$ if the integrand in (2fbx) is positive. In other words, ${\cal A}_{1}$ should be the set of all ${\cal O}$’s such that:

[TABLE]

If $U_{1}-\lambda C_{1}>0$ then (2fby) is equivalent to:

[TABLE]

This shows that we can use a criterion based on the likelihood ratio to maximise the utility. Note that if the requirement $U_{1}-\lambda C_{1}>0$ holds, then a non-empty ${\cal A}_{1}$ may be found because the benefit gained (or utility) from deciding for ${\cal A}_{1}$ outweighs the cost required. Otherwise, if $U_{1}-\lambda C_{1}\leq 0$, (2fby) suggests that ${\cal A}_{1}$ should be an empty set.

Note that (2fbx) can be interpreted as the maximisation of $P_{TP}$ subject to an upper bound on $P_{FP}$, which is the same class of optimisation formulation that the Neyman-Pearson lemma considers.

We remark that we can see from the above derivation that if both the mean utility and mean cost are linear in the probabilities $P_{TP}$, $P_{FP}$, $P_{TN}$ and $P_{FN}$, then the utility maximisation problem can be solved by using a criterion based on the likelihood ratio. We further remark that it is straightforward to generalise the above proof to the case with non-zero utility and cost.

Appendix D Necessity of delay

We use the method of contradiction to argue that there must be a delay in the cloud if $d_{0}>t_{\rm ss}$. Let us assume that there is a “memoryless” function $\psi({\rm E}[X_{*}(t)])$ which can carry out the computation in the cloud, i.e.

[TABLE]

Since we assume that the pathway (2) reaches the steady state by $d_{0}$, we can find time instants $t_{1}$ and $t_{2}$ where $t_{\rm ss}\leq t_{1}<d_{0}<t_{2}<d$ such that ${\rm E}[X_{*}(t)]$ is at steady state at both times $t_{1}$ and $t_{2}$, i.e. ${\rm E}[X_{*}(t_{1})]={\rm E}[X_{*}(t_{2})]$. Let us consider (2fca) at time instants $t_{1}$ and $t_{2}$, we have:

[TABLE]

which is a contradiction because ${\rm E}[X_{*}(t_{1})]={\rm E}[X_{*}(t_{2})]$. This establishes that, if $d_{0}>t_{\rm ss}$, there must be a delay in the cloud in Fig. LABEL:fig:compute:Lhat.