Modeling for Circuit Design

Our goal of modeling is to simulate the reactions that occur in our cell-free system with parameters represented by numerical data to predict the colorimetric outcomes of our biosensor. In order to achieve this outcome, we implemented ordinary differential equations (ODEs) that capture the essence of the protein dynamics in our system.

1. Cycle 0

ODEs

In our design, the oxidation process of bioamines by the constitutively promoted diamine oxidase (rDAO) that renders hydrogen peroxide (H₂O₂) is one of the most significant processes to activate the oxyS and the katG promoter.

The oxidation process of bioamines can be simplified as follows.

$Diamine \to H_{2}O_{2} \ (Reaction \ Rate = k_{1})$

This simplification can be done as we assume that the concentration of rDAO is constant or has reached equilibrium state when we are using the chassis.

$https://latex.codecogs.com/svg.image?k\left [ rDAO \right ]=k_1{}$

In our ODEs, rather than setting an initial concentration of bioamines, we assume a continuous input (m1) of diamine from our sample.

$\frac{d[Diamine]}{dt} = m_{1} - k_{1}[Diamine] - d_{1}[Diamine]\ \ \ \ \ \cdots (1)$

OxyR, the hydrogen peroxide stressor, is represented as OxyR_N when it is inactivated (N for non-activation) and OxyR_A (A for activation) when it is activated by H₂O₂. It is noted that only oxidized OxyR binds to a conserved OxyR binding site in the oxyS and katG regulatory region. Moreover, no degradation constant is rendered to PoxyS and Pkat as the number of plasmid in a bacteria is constant.

$H_{2}O_{2} + OxyR_{N} \to OxyR_{A}\ \ (Reaction \ rate = k_{2})$

$\frac{d[H_{2}O_{2}]}{dt} = k_{1}[Diamine] - k_{2}[OxyR_{N}][H_{2}O_{2}] - d_{2}[H_{2}O_{2}] \ \ \ \ \ \cdots (2)$

$\frac{d[OxyR_{N}]}{dt} = m_{2} - k_{2}[OxyR_{N}][H_{2}O_{2}] - d_{3}[OxyR_{N}] \ \ \ \ \ \cdots (3)$

$\frac{d[OxyR_{A}]}{dt} = k_{2}[OxyR_{N}][H_{2}O_{2}] - d_{4}[OxyR_{A}]\ \ \ \ \ \cdots (4)$

As the association and dissociation of the operator (PoxyS) and the transcription factor protein (OxyR_A) event occurs on a much faster time-scale than gene expression^[1], the association of OxyR_A and PoxyS to express mRNA molecules for fluorescent proteins can be treated in equilibrium.

$[OxyR_{A}] + [Unbound\ PoxyS] \rightleftharpoons [OxyR_{A} - bound \ PoxyS]$
$\textrm{(Reaction rate = } k_{3}, \textrm{Reverse reaction rate = }k_{-3} )$

As the rate of transcription from a regulated gene depends on the promoter occupancy, the rate of activated transcription (k₃) is proportional to the occupancy and can be represented as follows, when K = k_-3/k₃.

$\textrm{Rate of activated transcription }= \alpha \frac{[OxyR_{A}]/K}{1+[OxyR_{A}]/K}$

$\frac{d[PoxyS]}{dt} = -k_{3}[OxyR_{A}][PoxyS]\ \ \ \ \ \cdots (5)$

$\frac{d\left [ mRNA_{GFP}{} \right ]}{dt}=\frac{\alpha \left [ OxyR_A{} \right ]}{K+\left [ OxyR_A{} \right ]}-d_5{}\left [ mRNA_GFP{} \right ]\; \cdots \left ( 6 \right )$

*In the ODEs, [PoxyS] represents the concentration of OxyR_A bound operators. Similarly, when L = k_-4∕k₄, the expression of RFP can be described as follows.

$[OxyR_{A}] + [Unbound\ Pkat] \rightleftharpoons [OxyR_{A} - bound \ Pkat]$
$\textrm{(Reaction rate = } k_{4}, \textrm{Reverse reaction rate = }k_{-4} )$

$\frac{d[Pkat]}{dt} = -k_{4}[OxyR_{A}][Pkat]\ \ \ \ \ \cdots (7)$

$\frac{d\left [ mRNA_{RFP}{} \right ]}{dt}=\frac{\beta \left [ OxyR_A{} \right ]}{L+\left [ OxyR_A{} \right ]}-d_6{}\left [ mRNA_{RFP}{}\right ]\:\cdots \, \left ( 8 \right )$

*In the ODEs, [Pkat] represents the concentration of OxyR_A bound operators.

$mRNA_{GFP}\ \to \ GFP \ (Reaction\ rate: k_{5})$
$mRNA_{RFP}\ \to \ RFP \ (Reaction\ rate: k_{6})$

$\frac{d[GFP]}{dt} = k_{5}[mRNA_{GFP}]\ -\ d_{7}[GFP]\ \ \ \ \ \cdots (9)$

$\frac{d[RFP]}{dt} = k_{6}[mRNA_{RFP}]\ -\ d_{8}[RFP]\ \ \ \ \ \cdots (10)$

Kinetic rate constants

Rate constants for specific reactions are difficult to obtain and vary according to different environments including the availability of resources. In our model, reactions were categorized by either fast-occurring or slow-occurring processes, then the constants were adjusted to suit the desired outcomes.

Protein Degradation Rates

Protein degradation, an aspect of proteolysis, ensures maintenance and rapid regulation of individual protein concentration in a system^[2]. While the exact degradation of the specific protein in our cell-free system is unknown, the half-life of a protein (the time it takes for the concentration of the target protein to be reduced by 50%) provides a good reference for the mathematical model. The degradation rate of the protein is a key attributor in stabilizing the model while providing insight about the expiration of the colorimetric output biosensors, as the expression of fluorescent proteins will be degraded over time. Assuming the degradation rate of proteins is first order, then the degradation rate (k) is inversely correlated to the half life, where k can be represented as below.

$k = \frac{ln(2)}{t_{1/2}}$

Results and Interpretation

Figure 1.1.1 GFP and RFP concentration modeling outputs of Cycle 0 with Bioamines represented on a log-scale axis.

Figure 1.1.2 GFP and RFP concentration modeling outputs of Cycle 0 with Bioamines represented on a linear axis

Results of Cycle 0 prove the successful RFP expression at a higher bioamine concentration (Figure 1.1). However, the expression of both GFP and RFP render a color mixing issue as the final colorimetric output at a high bioamine is a combination of visible light wavelengths of 395 nm and 588 nm, which renders a yellow color.

The specific shade of color can be visualized via the RGB Color Addition simulation provided by the Physics Classroom^[3]. As this phenomenon does not meet our purpose of generating a red color at a higher bioamine concentration to warn the consumers, thus the circuit was improved to cleave the GFP under high concentrations of bioamines.

Figure 1.2 The simulation of the combination of GFP and RFP. The color in the middle represents the likely colorimetric output of Cycle 0.

Sensitivity Analysis

Introduction

As mentioned briefly in the introduction, the mathematical model is a schematic representation of a real-world process. In this section explaining the sensitivity analysis, we would like to emphasize the purpose of mathematical models that can be classified into three main categories. First, mathematical models are used to highlight the understanding of the mechanisms that govern the phenomena of interest. Second, mathematical models are used to predict the currently unknown state of a system. Third, mathematical models are used to produce a desirable condition by controlling components of a system. All three objectives deeply involve the implementation of input components of a real system, which are represented as ‘parameters’ in a mathematical model. The three sensitivity models implemented for Fisherly portray remarkable sensitivity to the amino acid sequences of proteins^[4], which facilitates experimental studies and the results allows us to reduce the cost or time associated with further experiments by narrowing down the significant parameters required to meet the expected performance criteria. Thus, the identification of the most influential parameters inside a mathematical model is a vital step to fulfill the purposes of the mathematical model. We decided to identify the influential parameters of our model using the following three global sensitivity analysis methods: Morris, Sobol, and eFAST.

The Morris Method

The Morris Method is the global sensitivity analysis developed by Max D. Morris as a preliminary screening measure for inputs of a computational model. The elemental structure of the Morris Method is designed as an individually randomized one-factor-at-a-time (OAT) test. As the name OAT suggests, the value of a single input parameter is differentiated to measure the corresponding elementary effect.

For a better illustration, let’s consider a computational model Y(x), that accepts the input vector X containing K inputs.

$X = [x_{1}, x_{2}, x_{3}, \cdots , x_{K-1}, x_{K} ]$

As the analysis is repeated for n times, the i-th run will start with the input vector X_i consisting of random input parameter values selected from the possible input space. Within each run, the value of the parameter x_k is changed to x_k + ∆ every step, while keeping the values of other K - 1 parameters the same. At the end of each step the elementary effect of k-th parameter, which is defined as below, is calculated.

$EE_{k} = \frac{y(x_{1}, x_{2}, \cdots , x_{k} + \Delta , \cdots , x_{K}) - y(x_{1}, x_{2}, \cdots , x_{K})}{\Delta }$

Even though the elementary effect alone is a local sensitivity analysis of a model, the Morris Method is still considered global sensitivity analysis as the final result is computed by averaging all elementary effects throughout the different points of the input space.

We chose the Morris Method as our first global sensitivity analysis tool because of its evident advantage as a preliminary screening measure. As Dr. Morris commented, the complicated computational models are often “time-consuming and hence expensive”. However, the Morris Method is easy-to-implement and easy-to-interpret global sensitivity analysis that requires relatively low computational power and time. Also, the Morris Method is independent of explicit assumptions of inputs having important effects, monotonicity of outputs with respect to inputs, or adequacy of a low-order polynomial model as an approximation to the computational model. The Morris Method was the most ideal choice to filter out insignificant parameters that can be neglected in the future engineering cycles.

The results of the Morris Method are shown below.

Figure 1.3.1 The elementary effects of the corresponding parameter at a specific running point

Figure 1.3.2 Morris Sensitivity analysis results

However, the outputs generated are ranks that simply indicate which of the parameters serve as the influential parameter of our model. Additionally, the results are prone to change depending on the number of simulations run. This can be attributed to the non-linear characteristic of our circuit, with multiple parameters influencing a single colorimetric output. In our trial, the number of maximum simulations was set to 70,000 times, and the most significant parameters were shown to be 16( $d_3{}$ ), 13( $m_1{}$ ), 14( $m_2{}$ ) whereas the moderate significant parameters were shown to be 7( $\beta$ ), 9( $L$ ), 8( $k$ ), and 6( $\alpha$ ). While $d_2{}$ and $m_1{}$ (for respective representations of these parameters, please refer back to the ODEs of Cycle 0) the finding of $d_2{}$ and $m_2{}$ is significant in that we are able to vary the expression of OxyR in our circuit by using a more or less sensitive promoter to render changes in the colorimetric output as we wish.

However, the Morris Method fails to give quantitative influence of each parameter on the output variation nor the insights about the non-linear effects of parameters on the output and interactions among different parameters. Therefore, we decided to move onto the next sensitivity analysis: the Sobol Method.

The Sobol Method

The Sobol Method is a form of a global sensitivity analysis that utilizes the Monte Carlo method as its basis for computation. Unlike the Morris method, the Sobol Method has the advantage of being able to provide multiple orders of Sensitivity Indices (also known as second order sensitivity indices), instead of a single interaction between two selected parameters. Moreover, in comparison with the eFAST which will be elaborated in the later section, it is more computationally efficient than the latter^[5].

Treat the model as a function where X is the vector of p model inputs, which are orthogonal to each other, and Y is the chosen model output^[6].

$https://latex.codecogs.com/svg.image?Y = f(X) = f(x_{1}, x_{2}, \cdots , x_{p})$
$https://latex.codecogs.com/svg.image?X = (x_{1}, x_{2}, \cdots , x_{p})$

Accordingly, the model output Y or f(X) can then be decomposed hierarchically as a sum of different functions when f₀, f_i(x_i), and f_i,j are defined as follows.

$f(X) = f_{0} + \sum_{i=1}^{p} f_{i}(x_{i}) + \sum_{1\leq i< j< p}^{p} f_{i,j}(x_{i}, x_{j}) + \cdots + f_{1, 2, \cdots , p} (x_{1}, x_{2}, \cdots , x_{p})\ \ \ \ \ \cdots (11)$

Then, we can take the square of either side of equation (11) to obtain the following equation.

$\int f^{2}(X)dX - f_{0}^{2} = \sum_{i=1}^{p} \sum_{k_{1}< \cdots < k_{i}}^{p} \int f_{k_{1}, \cdots , k_{i}}^{2} d(x_{1})\cdots d(x_{k_{i}})\ \ \ \ \ \cdots (12)$

Let’s define variance as follows.

$Var(z) = \int f^{2}(z)dz - \mathbb{P}(z)\ \ \ \ \ \cdots (13)$

Both sides of equation (12) can then be converted into probabilistic parameters to obtain the total variance as follows.

$Var(f(X)) = \sum_{i=1}^{p} Var(f_{i}(x_{i})) + \sum_{1\leq i< j< p}^{p} Var(f_{i, j}(x_{i}, x_{j})) + \cdots$
$+ Var (f_{1, 2, \cdots , p}(x_{1}, x_{2}, \cdots , x_{p})) \ \ \ \ \ \cdots (14)$

Upon obtaining equation (14), we can then define Sobol Indices (SI) as

$\textrm{First Order SI}\ \ \ S_{i} = \frac{Var(f(x_{i}))}{Var(Y)}$

$\textrm{Second Order SI} \ \ \ \ \ S_{i, j} = \frac{Var(x_{i, j})}{Var(Y)}$

$\textrm{Total Order SI} \ \ \ \ \ S_{i, j} = \frac{Var(x_{i, j, \cdots , p})}{Var(Y)}$

First Order SI is defined as the fraction of the variance of the output contributed only by parameter x_i over the total variance. To put it simply, S_i measures the effect of varying only x_i towards the model output Y, while it is normalized via the total output variance Var(Y). On the other hand, Total Order SI measures the effect of x_i and all of the interaction between parameters of any order with x_i towards Y and is normalled with Var(Y). A higher value of SI of a certain parameter implies that a minor change in this parameter may lead to a relatively bigger change in Y. In other words, Sobol Indices (SI) value indicates the sensitivity of a certain parameter that allows us to define the more influential parameters.

As the initial circuit has 2 outputs, 2 independent sobol sensitivity analysis for GFP and RFP respectively are done. We used the SALib package in Python to run our sensitivity analysis to obtain these results.

As the initial circuit has 2 outputs, 2 independent sobol sensitivity analysis for GFP and RFP respectively are done. We use SALib package in Python to run our sensitivity analysis to obtain these results:

Figure 1.4.1, 1.4.2 Sobol Analysis Result for cycle 0 ODE.
* is defined as SA values that are negative.

From these results, ignoring the degradation constants (d), we can obtain that the 6 parameters ( $https://latex.codecogs.com/svg.image?k_1{}$ , $https://latex.codecogs.com/svg.image?k_2{}$ , $https://latex.codecogs.com/svg.image?k_5{}$ , $https://latex.codecogs.com/svg.image?\alpha$ , $k$ , and $https://latex.codecogs.com/svg.image?m_2{}$ ) are sensitive towards GFP output, while 6 parameters ( $https://latex.codecogs.com/svg.image?k_1{}$ , $https://latex.codecogs.com/svg.image?k_2{}$ , $https://latex.codecogs.com/svg.image?k_6{}$ , $https://latex.codecogs.com/svg.image?\beta$ , $L$ , and $https://latex.codecogs.com/svg.image?m_2{}$ ) are significant for RFP output. From those parameters, there are three of them that are significant for both outputs: $https://latex.codecogs.com/svg.image?k_1{}$ , $https://latex.codecogs.com/svg.image?k_2{}$ and $https://latex.codecogs.com/svg.image?m_2{}$ .

After checking these three parameters (Figure 1.4.1, 1.4.2), we realized that any improvement towards them will not fix this color mixing issue. Changing those parameters produces similar graphs where both the GFP and RFP values increase as the parameter value increases. To solve this color mixing issue, we proposed to the circuit design team to either improve the activity of RBS that is upstream of RFP (increase $https://latex.codecogs.com/svg.image?k_6{}$ ), impeding the function of RBS that is upstream of GFP (decrease $https://latex.codecogs.com/svg.image?k_5{}$ ) or introduce a new system to solve the color mixing issue. After interviewing with Professor Yong Lai, the latter direction was chosen by the circuit design, and the TEV degradation system was proposed. Sobol analysis was done again for the improved circuit, Cycle 1, and the results are shown in the “Sobol Analysis of Cycle 1” section of this page.

Figure 1.5.1 Different $https://latex.codecogs.com/svg.image?k_1{}$ Parameters

Figure 1.5.2 Different $https://latex.codecogs.com/svg.image?k_2{}$ Parameters

Figure 1.5.3 Different $https://latex.codecogs.com/svg.image?m_2{}$ Parameters

Figure 1.5.4 Changing $https://latex.codecogs.com/svg.image?k_5{}$ and $https://latex.codecogs.com/svg.image?k_{6}{}$ Parameters

From these results, we can obtain that the 4 parameters ( $https://latex.codecogs.com/svg.image?k_4{}$ , $https://latex.codecogs.com/svg.image?k_5{}$ , $https://latex.codecogs.com/svg.image?d_1{}$ , and $https://latex.codecogs.com/svg.image?d_2{}$ ) are important in GFP output, while 5 parameters ( $https://latex.codecogs.com/svg.image?k_7{}$ , $https://latex.codecogs.com/svg.image?k_8{}$ , $https://latex.codecogs.com/svg.image?g_1{}$ , $m$ , $https://latex.codecogs.com/svg.image?d_3{}$ , and $https://latex.codecogs.com/svg.image?d_4{}$ ) with SI values above 0.05 being significant.

eFAST

(E)eFAST, or the extended Fourier Amplitude Sensitivity Test is in an updated version of the FAST or the Fourier Amplitude Sensitivity Test. Hence, the understanding of the FAST must come before looking into the details of eFAST. Just like the Sobol Method, the FAST is a variance-based sensitivity analysis method. However, a key difference from the Sobol method is where the computation for SA in eFAST is achieved using spectral analysis, instead of Monte Carlo method (5). Let’s say $y=f(X)$ to be the model of interest, where y represents the output of the model, $f$ is the function of the model, and $https://latex.codecogs.com/svg.image?X=(x_1{},x_2{}, ..., x_n{})$ is the input vector of the model composed of $n$ input factors. The input vector $X$ is assumed to be random with the probability density function of $https://latex.codecogs.com/svg.image?P(X)=P(x_1{},x_2{}, ..., x_n{})$ .

Furthermore, the domain of the input vector is regarded as a unit hypercube.

$https://latex.codecogs.com/svg.image?K_x{}^{n}=\left\{ x_1{},...,x_n{}|x_1{}\sim f_i{x_i{}}, i=1,...,n\right\}$

The input values of the FAST analysis are sampled from the input space via periodic search functions. Each input value is designated with a unique integer frequency, which is later used to determine the sensitivity indices of the input parameter. The general form of a periodic search function is defined as below.

$https://latex.codecogs.com/svg.image?x_i{(s)}= G_i{(sin\: w_i{s})}, \forall i=1,2,...,n$

In the general from, s is a scalar variable varying over the range − ∞ < s < + ∞. G_i are transformation function and {ω_i} is a set of different frequencies. For the initial FAST analysis, the search function was set as the following equation.

$https://latex.codecogs.com/svg.image?x_i{(s)}= \frac{1}{2}+\frac{1}{\pi }{arcsin(sin\: w_i{s}+\varphi _i{})}, \forall i\in [1,n]$

Both the scalar variable s and the random shift parameter $https://latex.codecogs.com/svg.image?\varphi _i{}$ have an uniform distribution over (0,2π]. Then, $https://latex.codecogs.com/svg.image?x_i{}\sim p_{x_i}(x_i)$ is obtained through the following transformation.

$https://latex.codecogs.com/svg.image?x_i{(s)}= F^{-1}(u_i{}), \forall i= 1,...,n$

The Fourier coefficients, derived by the evaluation of $f(x(s))$ , are represented as $https://latex.codecogs.com/svg.image?c_\omega {}$ .

$https://latex.codecogs.com/svg.image?c_\omega {}=\frac{1}{\pi}\int_{0}^{2\pi}f(x(s))e^{i\omega s}ds, \forall \omega \in N^{*}$

Lastly, the first-order sensitivity indices or the “main effect” of a parameter is computed according to the following equation.

$https://latex.codecogs.com/svg.image?S_{x_i}{}^{FAST}=\frac{\sum_{l=1}^{+\infty }\left| c_l_{\omega_i}{}\right|^{2}}{\sum_{\omega=1}^{+\infty }\left| c_{\omega_i}{}\right|^{2}}, \; \forall i=1,...,n$

Despite the advantage of FAST in terms of lighter computational burden, the classical FAST was only capable of computing the “main effect”, assuming the input parameters to be independent of each other. The lack of capacity to derive the “Total Sensitivity Index” was the major drawback of the classical FAST compared to the Sobol Method. The eFAST resolved the limitation of the classical FAST as the eFAST derives the total sensitivity index of $https://latex.codecogs.com/svg.image?x_i{}$ or $https://latex.codecogs.com/svg.image?ST_{x_i}{}$ by assigning higher frequency to $https://latex.codecogs.com/svg.image?x_i{}$ and low frequencies to other parameters $https://latex.codecogs.com/svg.image?x_{\sim i}{}$ . The total sensitivity index is defined as the following equation.

$https://latex.codecogs.com/svg.image?ST_{x_i}=\frac{\sum_{\omega =\frac{\omega _i{}}{2}+1}^{N/2 }\left| c_{\omega_i}{}\right|^{2}}{\sum_{\omega=1}^{N/2 }\left| c_{\omega_i}{}\right|^{2}}$

The results of the eFAST analysis are shown below.

Figure 1.6 eFAST Sensitivity Analysis results for Cycle 0

As shown in Figure 1.6, $m_2{}$ and $d_3{}$ were shown to be significant parameters, where the rest of the parameters had minimal effect. As $d_4{}$ is a parameter that represents the degradation rate of OxyR, we have disregarded the constant in our interpretation. Nevertheless, this is consistent with the findings of both the Sobol analysis and the Morris analysis method that $m_2{}$ is a significant parameter.

Figure 1.7 Varying $m_2{}$ parameters in Cycle 0

As seen in Figure 1.7, differing the parameters of $m_2{}$ by 10% brought about visible changes in the outputs. Combined with the results of our optimization results, we will be able to adjust our circuit components in Cycle 0 to obtain the most yield of the fluorescent proteins. However, we hope to serve as a reminder that all three sensitivity analyses were performed based on similar constructions of ODEs and parameter boundaries that may have led to concurrent results. Thus, the conclusions should further be clarified by experimental data to confirm the significant influences of varying specific element concentrations.

2. Cycle 1

The level of OxyR protein does not change in cells treated with H₂Osub>2, an indication that it is activated post-translationally. In vitro experiments demonstrated that OxyR can exist in either oxidized or reduced forms. Only the oxidized form binds DNA at the promoters of katG, dps, gorA, and so on, where activated OxyR stimulates transcription by interacting with the carboxyl-terminal domain of the α-subunit of RNA polymerase.

The key characteristic of Cycle 1 is the additional expression of the TEV protease which acts as the inactivation of the GFP. Both GFP and RFP are expressed by the constitutive promoter. Both the RFP and the GFP is constitutively promoted and the biosensor renders only the red color when there is high concentration of bioamines due to the cleavage of GFP by TEV protease. We designed GFP and OxyR to be expressed originating from the same promoter, and thus the same rate constant m3 was used for continuous input. Additionally, GFPt (visible) refers to the TEV cleavage site–tagged GFP, and the cleaved form is represented as GFPd. The RFPt refers to the TEV cleavage site–tagged RFP, and the cleaved form is represented as RFPd(visible).

$https://latex.codecogs.com/svg.image?\frac{d\left [ Diamine \right ]}{dt}=m_1{}-k_1{\left [ Diamine \right ]}$

$https://latex.codecogs.com/svg.image?\frac{d\left [H_2{}O_2{} \right ]}{dt}=k_1{\left [ Diamine \right ]}-k_2{\left [ OxyR_N{} \right ]\left [ H_2{}O_2{} \right ]}-d_1[H_2{}O_2{}]$

$https://latex.codecogs.com/svg.image?\frac{d\left [OxyR_N{} \right ]}{dt}=m_2-k_2{\left [ OxyR_N{} \right ]\left [ H_2{}O_2{} \right ]}-d_2[H_2{}O_2{}]$

$https://latex.codecogs.com/svg.image?\frac{d\left [OxyR_A{} \right ]}{dt}=k_2{\left [ OxyR_N{} \right ]\left [ H_2{}O_2{} \right ]}-d_3[OxyR_A]$

$\frac{d\left [mRNA_{TEV}{} \right ]}{dt}=\frac{\gamma \left [ OxyR_A{} \right ]}{S+\left [ OxyR_A{} \right ]}-d_4{}\left [ mRNA_{TEV}{} \right ] \; \; \left ( S=\frac{k_{-2}{}}{k_2{}} \right )$

$https://latex.codecogs.com/svg.image?\frac{d\left [TEV \right ]}{dt}=k_6{\left [ mRNA_{TEV}{} \right ]-d_5[TEV]$

$https://latex.codecogs.com/svg.image?\frac{d\left [GFPt \right ]}{dt}=m_3{}-k_7{\left [ GFPt \right ][TEV]-d_6[GFPt]$

$https://latex.codecogs.com/svg.image?\frac{d\left [GFPd \right ]}{dt}=k_8{\left [ GFPt \right ][TEV]-d_7[GFPd]$

$https://latex.codecogs.com/svg.image?\frac{d\left [RFPt \right ]}{dt}=m_4{}-k_9{\left [ RFPt \right ][TEV]-d_8[RFPt]$

$https://latex.codecogs.com/svg.image?\frac{d\left [RFPd \right ]}{dt}=k_{10}{\left [ RFPt \right ][TEV]-d_9[RFPd]$

Constant chart

Kinetic Reaction Rates

Reaction	Rate Constant	Estimated Value	Assumption
$https://latex.codecogs.com/svg.image?Diamine\rightarrow H_2{}O_2{}$	$https://latex.codecogs.com/svg.image?k_1{}$	$https://latex.codecogs.com/svg.image?9.63\times 10^{-5} (per\; sec)$	Half life: 2 hours^[7]
$https://latex.codecogs.com/svg.image?H_2{}O_2{}+OxyR_N{}\rightarrow OxyR_A{}$	$https://latex.codecogs.com/svg.image?k_2{}$	$https://latex.codecogs.com/svg.image?4.20\times 10^{-4} (per\; sec)$	Second order reaction^[8]
$https://latex.codecogs.com/svg.image?OxyR_A{}+katGp\leftrightarrow mRNA_{TEV_p}{}$	$https://latex.codecogs.com/svg.image?k_3{}$	No reference	N/A
$https://latex.codecogs.com/svg.image?mRNA_{TEV_p}{}\rightarrow TEV_p{}$	$https://latex.codecogs.com/svg.image?k_4{}$	No reference	N/A
$https://latex.codecogs.com/svg.image?GFPt + TEV_p{}\rightarrow GFPd+TEV_p{}$	$https://latex.codecogs.com/svg.image?k_5{}$	No reference	N/A
$https://latex.codecogs.com/svg.image?RFPt + TEV_p{}\rightarrow RFP+TEV_p{}$	$https://latex.codecogs.com/svg.image?k_6{}$	No reference	N/A

Reaction	Rate Constant	Estimated Value	Assumption
$https://latex.codecogs.com/svg.image?\phi\rightarrow Diamine$	$https://latex.codecogs.com/svg.image?m_1{}$	0	Active transport;First order reaction
$https://latex.codecogs.com/svg.image?\phi\rightarrow OxyR_N{}$	$https://latex.codecogs.com/svg.image?m_2{}$	0.1	Controllable
$https://latex.codecogs.com/svg.image?\phi\rightarrow GFPt$	$https://latex.codecogs.com/svg.image?m_3{}$	0.53	Controllable
$https://latex.codecogs.com/svg.image?\phi\rightarrow RFPt$	$https://latex.codecogs.com/svg.image?m_4{}$	0.10	Controllable

Degradation Rates

Protein	Constant	Value	Half-Life
$https://latex.codecogs.com/svg.image?H_2{}O_2{}\rightarrow \phi$	$https://latex.codecogs.com/svg.image?d_1{}$	$https://latex.codecogs.com/svg.image?8.02\times 10^{-6}\; (per\; sec)$	24 hours^[9]
$https://latex.codecogs.com/svg.image?OxyR_N{}\rightarrow \phi$	$https://latex.codecogs.com/svg.image?d_2{}$	$https://latex.codecogs.com/svg.image?1\; (per\; sec)$	No reference
$https://latex.codecogs.com/svg.image?OxyR_A{}\rightarrow \phi$	$https://latex.codecogs.com/svg.image?d_3{}$	$https://latex.codecogs.com/svg.image?1\; (per\; sec)$	No reference
$https://latex.codecogs.com/svg.image?mRNA_{TEV}{}\rightarrow \phi$	$https://latex.codecogs.com/svg.image?d_4{}$	$https://latex.codecogs.com/svg.image?2.10\times 10^{-3}\; (per\; sec)$	5.5 minutes^[10]
$https://latex.codecogs.com/svg.image?TEV\rightarrow \phi$	$https://latex.codecogs.com/svg.image?d_5{}$	$https://latex.codecogs.com/svg.image?6.42\times 10^{-5}\; (per\; sec)$	No reference
$https://latex.codecogs.com/svg.image?GFPt\rightarrow \phi$	$https://latex.codecogs.com/svg.image?d_6{}$	$https://latex.codecogs.com/svg.image?6.88\times 10^{-5}\; (per\; sec)$	2.8 hours^[11]
$https://latex.codecogs.com/svg.image?RFPt\rightarrow \phi$	$https://latex.codecogs.com/svg.image?d_7{}$	$https://latex.codecogs.com/svg.image?1.38\times 10^{-4}\; (per\; sec)$	No reference; same as RFPt
$https://latex.codecogs.com/svg.image?GFPd\rightarrow \phi$	$https://latex.codecogs.com/svg.image?d_8{}$	$https://latex.codecogs.com/svg.image?1.38\times 10^{-4}\; (per\; sec)$	No reference; sames as GFPd
$https://latex.codecogs.com/svg.image?RFPd\rightarrow \phi$	$https://latex.codecogs.com/svg.image?d_9{}$	$https://latex.codecogs.com/svg.image?9.63\times 10^{-6}\; (per\; sec)$	20 hours^[12]

Results and Interpretation

The outputs of Cycle 1 successfully proves the distinct color segregation under varying bioamine concentrations.

Figure 2.1.1 GFP and RFP concentration modeling outputs of Cycle 1 with Bioamines represented on a log-scale axis. $https://latex.codecogs.com/svg.image?k_1{}$ Parameters

Figure 2.1.2 GFP and RFP concentration modeling outputs of Cycle 1 with Bioamines represented on a linear axis. $https://latex.codecogs.com/svg.image?k_2{}$ Parameters

Comparison of high bioamine and low bioamine environments

In this section, we are going to implement realistic bioamine level thresholds and compare the colorimetric outputs accordingly. Our aim is to generate a green color when the fish is edible, and red when it is not safe to consume. Assuming the regulation for total bioamines resides around 1000 ppm, the following trend can be observed as follows.

A further analysis of the colorimetric output generated at different bioamine thresholds proves that the mechanism of generating different colored outputs functions. As seen from the outputs of our sensitivity analysis, we have good understanding of the significant parameters, implying the flexibility of our biosensor to be altered in certain parts to meet different requirements. While histamine is regulated to 200 ppm in many countries, the exact regulation standard for total bioamines (TBA) is unknown. Thus, our temporary TBA threshold of 1000 ppm is likely to be altered through further specification of local and international regulations in the future. Additionally, our biosensor holds the significance that it can be built to detect bioamines not limited to seafood, but other types of meat as well. The colorimetric outputs can be modified to take into account the various regulations of total bioamines via adding preliminary histamine into our solution buffer and increasing or decreasing the promoter’s sensitivity to bioamines.

Sobol Analysis of Cycle 1

Sobol analysis was performed again with the ODEs set for cycle 1 to provide insights not only on the improvement of cycle 1, but also on the future direction that the circuit design team can hold in order to further improve our circuit. The corresponding sensitivity analysis results are shown as follows.

Figure 2.2Sobol Analysis of Cycle 1 on GFP and RFP respectively

*is defined as SA values that are negative.

Interestingly, the Sobol Analysis of cycle 1 shows that there are no sensitive parameters that are commonly found in both GFP and RFP outputs. Assuming that degradation terms are ignored, this implies that changing GFP sensitive parameters (such as $https://latex.codecogs.com/svg.image?m_3{}$ ) would not affect the RFP performance, while changing RFP sensitive parameters (such as $https://latex.codecogs.com/svg.image?k_1{}$ , $https://latex.codecogs.com/svg.image?k_2{}$ , $https://latex.codecogs.com/svg.image?k_4{}$ , $https://latex.codecogs.com/svg.image?k_6{}$ , $https://latex.codecogs.com/svg.image?\alpha$ , $https://latex.codecogs.com/svg.image?k_7{}$ , $https://latex.codecogs.com/svg.image?m_2{}$ , $https://latex.codecogs.com/svg.image?m_4{}$ and $https://latex.codecogs.com/svg.image?d_4{}$ ) would not affect the GFP performance. These observations are further elaborated in Figure 2.3.1 and Figure 2.3.2 below.

Figure 2.3.1 Modeling of GFP with varying m3 Parameter Values

Figure 2.3.2 Modeling of RFP with varying k1 Parameter Values

Moreover, the results from Figure 2.4.2 are extremely significant for us. The parameter under observation, $https://latex.codecogs.com/svg.image?k_1{}$ is the rate constant for the reaction of the conversion of diamine to $H_2{}O_2{}$ , with $https://latex.codecogs.com/svg.image?k_1{}=k\left [rDAO \right ]$ . Thus, by increasing the rDAO concentration, we could increase $https://latex.codecogs.com/svg.image?k_1{}$ , hence increase the RFP output. We shared our observation with the circuit design, and they found out that pelB tag could increase the concentration of functional rDAO in our system. Refer to our design page for more information.

Circuit Design Experiment

Another method to increase the concentration of rDAO is to change our product to a cell free system. In CFS, rDAO can be pre-deposited in a certain amount by transforming it into BL21 cells, which serve as the source of cell lysates, and harvesting them in the log phase to ensure their active production of rDAO. Moreover, this also substantially reduces the variability of live cell productivity which would cause the ambiguity of k1. Based on these observations, the idea of increasing the concentration of rDAO was also shared with the chassis team to increase the further functionality of our biosensor.

Summary and Discussions

If we were asked to increase our GFP or positive control output while maintaining our RFP output (in other words, increasing only m3 as seen in Figure 2.4.1), this raises the color mixing issue again. Therefore, the dry lab team proposed another system in addition to our TEV degradation system, which is the post-translational inactivation system. The system that the dry lab team suggests is a transcriptional level deactivation system could be introduced in the future, such as the tetracycline-controlled Tet-Off-Tet-On gene expression systems.

Based on these results, mathematical modeling provides the optimal language for describing complex systems, but it only operates with the information fed into it. The assumptions and parameters are favored to meet the modeling expectations we initially envisioned. We hope that modeling provides a better understanding of the mechanism of our biosensor and that it provides a guideline for the experiments.

Supplementary Data and Codes

The supplementary data and codes can be found on the Github page of HKUST iGEM 2022.

References

[1] Ingalls, B. P. (2022). 7.1.2 Regulated gene expression. In Mathematical Modeling in Systems Biology: An introduction (pp. 192–193). essay, The MIT Press.

[2] Zhang, L., Gurskaya, N. G., Merzlyak, E. M., Staroverov, D. B., Mudrik, N. N., Samarkina, O. N., Vinokurov, L. M., Lukyanov, S., & Lukyanov, K. A. (2018). Method for real-time monitoring of protein degradation at the single cell level. BioTechniques, 42(4), 446–450. https://doi.org/10.2144/000112453

[3] Physics simulation: RGB color addition. The Physics Classroom. (n.d.). Retrieved September 26, 2022, from https://www.physicsclassroom.com/Physics-Interactives/Light-and-Color/RGB-Color-Addition/RGB-Color-Addition-Interactive

[4] Ratushny AV, Ramsey SA, Aitchison JD. Mathematical modeling of biomolecular network dynamics. Methods Mol Biol. 2011;781:415-33. doi: 10.1007/978-1-61779-276-2_21. PMID: 21877294; PMCID: PMC4482237.

[5] Vazquez-Cruz, M. A., Guzman-Cruz, R., Lopez-Cruz, I. L., Cornejo-Perez, O., Torres-Pacheco, I., & Guevara-Gonzalez, R. G. (2014). Global sensitivity analysis by means of EFAST and Sobol’ methods and calibration of reduced state-variable TOMGRO model using genetic algorithms. Computers and Electronics in Agriculture, 100, 1–12. https://doi.org/https://doi.org/10.1016/j.compag.2013.10.006

[6] Zhang, X. Y., Trame, M. N., Lesko, L. J., & Schmidt, S. (2015). Sobol Sensitivity Analysis: A Tool to Guide the Development and Evaluation of Systems Pharmacology Models. CPT: pharmacometrics & systems pharmacology, 4(2), 69–79. https://doi.org/10.1002/psp4.6

[7] Naila, A., Flint, S., Fletcher, G., Bremer, P., Meerdink, G., & Morton, R. (2012). Prediction of the amount and rate of histamine degradation by diamine oxidase (DAO). Food Chemistry, 135(4), 2650-2660. https://doi.org/10.1016/j.foodchem.2012.07.022

[8] Pedre, B., Young, D., Charlier, D., Mourenza, A., Rosado, L., Marcos-Pascual, L., Wahni, K., Martens, E., Rubia A., Belousov, V., Mateos, L. (2018). Structural snapshots of OxyR reveal the peroxidatic mechanism of H2O2 sensing. Proceedings of the National Academy of Sciences. https://doi.org/10.1073/pnas.180795411

[9] Breithaupt, J. (2007). Summary Review of Available Literature for Hydrogen Peroxide and Peroxyacetic Acid for New Use to Treat Wastewater. United States Environmental Protection Agency. https://www3.epa.gov/pesticides/chem_search/cleared_reviews/csr_PC-000595_12-Jul-07_a.pdf

[10] Bernstein, J., Khodursky, A., Lin, P., Lin-Chao, S., Cohen, S. (2002). Global analysis of mRNA decay and abundance in Escherichia coli at single-gene resolution using two-color fluorescent DNA microarrays. PNAS. https://doi.org/10.1073/pnas.112318199

[11] Halter, M., Tona, A., Bhadriraju, K., Plant, A. L., & Elliott, J. T. (2007). Automated live cell imaging of green fluorescent protein degradation in individual fibroblasts. Cytometry. Part A : the journal of the International Society for Analytical Cytology, 71(10), 827–834. https://doi.org/10.1002/cyto.a.20461

[12] He, L., Binari, R., Huang, J., Falo-Sanjuan, J., & Perrimon, N. (2019). In vivo study of gene expression with an enhanced dual-color fluorescent transcriptional timer. eLife, 8, e46181. https://doi.org/10.7554/eLife.46181