In our project, in normal state, HucR inhibits the promoter P_hucO, which results in low fluorescence intensity. However, when uric acid exists, it binds to HucR and the inhibition of transcription will be weakened, and RFP will be expressed so that there will be fluorescence.

One of the demands of our project is to plot a curve describing the relationship between the concentration of uric acid (denoted as $[UA]$ below) and the intensity of fluorescence (which is related to the concentration of RFP, and is denoted as $[RFP]$ below). According to our experiment, our result can be fitted well using quadratic function (see below), but quadratic function cannot justify the dynamics, and irrational negative results will happen for high $[UA]$. Therefore, we want to apply a more reasonable dynamic model to fit our results.

The Hill equation describes the fraction of a macromolecule saturated by a ligand as a function of the concentration of the ligand; it is used to determine the degree of contractibility of the receptor bound to the enzyme or receptor. This equation was first elucidated by Archibald Hill in 1910. Hill's equation as a description of the relationship between the concentration of compound adsorbed to the binding site, and the occupied fraction of the binding site is given by $$\theta = \dfrac{[LR]}{[R]+[LR]},$$ where $[L]$ denotes the concentration of unbound ligand, $[R]$ denotes the concentration of receptor, and $[LR]$ denotes the concentration of receptor with bound ligand. If we introduce the apparent dissociation $K_d$, then at equilibrium state, we have $$\theta = \dfrac{[L]^n}{K_d+[L]^n},$$ where $n$ is Hill coefficient.

For $n > 1$, a positive synergism is considered, i.e. once a ligand molecule has bound to a receptor, the affinity of that receptor for other ligands increases.

For $n = 1$, no synergism is considered, i.e., the affinity of the receptor for a ligand molecule does not depend on whether or not a ligand molecule has bound to it.

For $0 < n < 1$, a negative synergism is considered, i.e. once a ligand molecule is bound to the receptor, the affinity of that receptor for other ligands decreases.

We can rewrite the equation to the form below: $$\theta = \dfrac{[L]^n}{K_a^n+[L]^n}.$$ In this form, $K_a$ means the value of $[L]$ to make $\theta=\dfrac{1}{2}$, a.k.a., the microscopic dissociation constant.

The activity of P_hucO is regulated by uric acid indirectly, but to simplify our model, we treat it as direct regulation. According to [1], this effect can be described by Hill equation with $n=1$, therefore, for some constant $A$, the ratio of "available" promoter is $$\theta_1 = \dfrac{[UA]}{A+[UA]}.$$

Since the mRNA concentration is low, transcription and translation can be considered as a linear process, that is, $$\dfrac{d[mRNA]}{dt}=v_{p-mRNA}\theta_1-v_{d-mRNA}[mRNA].$$ Therefore, in steady state, let $\dfrac{d[mRNA]}{dt}=0$, we have $[mRNA]=B\theta_1$ for some constant $B$.

For translation of hrpR and hrpS, we have $$\dfrac{d[hrpR]}{dt}=k_2\theta_1-k_{d1}[hrpR]-k_+[hrpR]^2+k_-[hrpRS],$$ and for the binding, as it always holds $[hrpR]=[hrpS]$, we have $$\dfrac{d[hrpRS]}{dt}=k_+[hrpR]^2-k_-[hrpRS]-k_{d2}[hrpRS].$$ If $[hrpRS]$ is much lower than $[hrpR]$, then we can neglect $-k_+[hrpR]^2+k_-[hrpRS]$ in the expression of $\dfrac{d[hrpR]}{dt}$, and from this, if we let $\dfrac{d[hrpR]}{dt}=\dfrac{d[hrpRS]}{dt}=0$, we have $[hrpRS]=B'\theta_1^2$ for some constant $B'$.

The dynamics of RFP expression is similar, that is, for some constant $C$ and $D$, we have $$[RFP]=\dfrac{C[hrpRS]}{D+[hrpRS]},$$ and the combined result of our equation is $$[RFP]=\dfrac{C\theta_1^2}{D'+\theta_1^2}=\dfrac{C[UA]^2}{D'(A+[UA])^2+[UA]^2}=\left(\dfrac{E_1}{[UA]^2}+\dfrac{E_2}{[UA]}+E_3\right)^{-1}$$ for some constant $E_1$, $E_2$ and $E_3$.

We apply our model to experimental result. As our equation can be rewritten to $$\dfrac{1}{[RFP]}=\dfrac{E_1}{[UA]^2}+\dfrac{E_2}{[UA]}+E_3,$$ we can apply double reciprocal plot. Using polyfit in numpy, we obtained that $$E_1=13.66,\ E_2=0.1648,\ E_3=1.273\times 10^{-4}.$$ The curve is shown below.

The graph above demonstrates that the curve fitted by our model has a similar property to the quadratic function one: in both of models, OD600 value increase less quickly as $[UA]$ increases from $200\mu$mol/L to $1000\mu$mol/L. However, our model solves a drawback: as $E_1$, $E_2$ and $E_3$ are all positive, OD600 value will always increase when $[UA]$ increases, therefore, the irrational negative results can be avoided.

Moreover, if we let $\dfrac{d^2[RFP]}{d[UA]^2}=0$, we find that our system will reach the highest sensitivity to change of $[UA]$ when $[UA]=149.47\mu$mol/L, which is near the lower bound of our operating range. In addition, the sensitivity is good even when $[UA]=3000\mu$mol/L, indicating that our project will have a good performance in practice.

import numpy as np
import matplotlib .pyplot as plt
ua = # original data, omitted
rfp = # original data, omitted
inv_ua = 1.0 / ua
inv_rfp = 1.0 / rfp
poly = np.polyfit(inv_ua, inv_rfp, deg=2)
print(poly)
plt.scatter(ua, rfp, color="crimson", marker="+", label="experiment")
ua_fit = np.linspace(0.5, 3000.5, 3000)
rfp_fit = 1.0 / (poly[0]/ua_fit/ua_fit+poly[1]/ua_fit+poly[2])
plt.plot(ua_fit, rfp_fit, linestyle="--", color="darkgreen", label="fit")
plt.legend(loc="lower right")
plt.xlabel("[UA](μmol/L)")
plt.ylabel("[RFP]/OD600")
plt.title("Working curve about our device")
plt.show()


[1] Apparat, A., Jalihal, A. P., Law, J. N., Bharadwaj, A. & Murali, T. Benchmarking algorithms for gene regulatory network inference from single-cell transcriptomic data. Nature methods 17, 147–154 (2020).