Model

Overview

Computational modelling of "PERspectives" has been proven to be crucial in order to accelerate and improve the wet lab tasks, aiming to achieve significant and reproducible results. Moreover, by setting the basis of the computational aspect of PERspectives, we further establish a potential workflow for future applications both on our as well as on relevant fields.

Our computational approach was based on modelling the senders and the recievers populations seprately. The aim of the model is to predict the change in concentrations of key chemical substances, the rate of gene expression and to evaluate the regulatory capacity of our system. Also, based on the experimental work carried out by the wet lab, we determined some kinetic constants to create a customised kinetic model for our system.

But we did not limit ourselves to just creating a kinetic model. For our system it is also important to observe the effect of lactone diffusion as well as the stability analysis.

In the following sections, we describe an outline of the aforementioned models and, most importantly, the thinking process we followed during their development. We hope they will serve as a helpful companion to anyone who wishes to use "PERspectives" in the future. As everything else in iGEM, our work is open to recommendations and improvements from the SynBio community.

Senders
Receivers
Stability Analysis
Pattern Recognition
Diffusion

Senders System Explained




pTeT Promoter


The purpose of our team is to create 3x1 patterns of E.coli senders subpopulations that express the LuxI gene at a different, desirable rate.
To achieve this, we designed a regulatory system based on the pTet promoter. Given that there is no aTc in the system, the promoter is repressed by the TetR regulator. If aTc is present in the system, then it binds with the TetR, creating a complex, that is pumped out of the cell, allowing the mRNA of the following gene to be produced. For the production of the TetR our team used a constitutive promoter from the iGEM Registry.

Simulating the weights of the neural network with RBSs


RBS Library Generation

Our goal was to find a way to simulate the weights of a perceptron algorithm with a biological analogue. In cooperation with the wet lab members, it was decided that the weights of the neural network would be simulated with Ribosome Binding Sites (RBSs) that assign different translation rates. In order to accomplish this task, the Salis Lab RBS Calculator Tool was used. This tool allows the user to see the translation rate of an RBS, but also taking into account the following gene, since it is known that the translation rate of an RBS depends on it. There are two ways to use the RBS Calculator:

1. Upload an RNA sequence (RBS + gene) and receive as a result the translation rate of this RBS when followed by this specific gene.
2.Demand a specific translation rate of a specific gene and receive as a result RBS sequences that will provide it.

You can find out more on how to use RBS Calculator by visiting the documentation page of the tool, or by reading our Tutorial that you can find here. After studying on and consulting experts on the subject. (See Integrated Human Practices) Our team used the first way to calculate the translation rate of the following RBS Libraries:

1. Anderson RBS Library from iGEM Registry
2. RBS Library used in 'Precise and reliable gene expression via standard transcription and translation initiation elements' paper’.

The need for creating new synthetic RBS sequences that would fit exactly to the design of our project emerged. So the dry lab members enriched the library that was built based on existing RBSs, with new synthetic ones. The sequences of the new RBSs where found by using again Salis Lab RBS Calculator tool, but these time by requesting a specific translation rate and receiving the synthetic sequence by the tool.

Once the library with all the 166 RBSs was ready and all the translation rates were calculated, the RBSs were divided into three classes:

1. The strong intensity RBSs intensity RBSs, with translation rate more than 90000 au.
2. The medium intensity RBSs, with translation rate around 25000-40000 au.
3. The weak intensity RBSs, with translation rate less than 15000 au.

Practically, therefore, a library of RBSs emerged from which the strong ones would simulate hight weights, those with a moderate translation rate a moderate weight while the weak RBSs simulate the low weights for the perceptron.

Download our RBS Library



Finally, the RBSs that we decided to test in the lab are the following:

Name Translation Initiation Rate (au) Strength
BCD2 287840.97 Very Strong
S4 85298.9 Strong
S5 78686.2 Strong
S10 31126.5 Medium
S11 30580.8 Medium
S12 11058.3 Weak
S13 10752.1 Weak
BCD1 9455.2 Weak
BCD8 2190.04 Weak
BSD12 12782.75 Weak
BCD14 22019.63 Weak
Table 1: Final RBSs chosen for experimental use.


Using Quorum Sensing to stimulate the weighted sum of the neural network


In order to be able to simulate the weighted sum of the perceptron, we took advantage of the natural way bacteria communicate, the quorum sensing mechanism. Thus, our senders subpopulation must be able to express the LuxI gene. The expression of the gene is regulated by the aTc promoter mentioned above. In this way, if we wanted to stimulate a small weight, we would use the aTc promoter followed by the LuxI gene with a low intensity RBS. Moreover we would do the same thing respectively to simulate a medium or high perceptron weight. In this way every subpopulation of the senders would produce a specific amount of AHL which is determined by the RBS we used. By combining them together, we would have achieved the simulation of the weighted sum of the perceptron by adding up the total amount of lactone produced by each subpopulation.

Senders' Circuit Overview: TetR is expressed by a constitutive promoter (BBa_J23106 and PlacI were tested) and the BBa_B0034 RBS. It binds to specific operators in the Ptet promoter and blocks RNA polymerase binding and, thus, transcription. Please refer to Engineering Success page for more details.

Figure 1: Senders’ Circuit Overview.

Biochemical Reactions


For the kinetic modelling of the senders' system, we had to convert all the aforementioned design into biochemical equations, each of which is characterized by its own rate of consumption or production of the chemical species. The resulting equations are:

1. ∅ → TetR
2. TC + TetR→TetRTC
3. TetRTC + Ptetoff ↔ Pteton
4. Pteton → Pteton + mRNALuxI
5. TetR + TetO → Ptetoff
6. mRN ALuxI → LuxI
7. TC→ ∅
8. TetRTC → ∅
9. mRNALuxI → ∅


Equation (1) describes the production of TetR by the constitutive promoter. Equation (2) describes the formation of the TetR+aTc comlex when aTc is present in the system and this formation leads to equation (3) describing the activation of the Ptet promoter. Equations (4) and (6) describe the translation of the gene and finally the production of the LuxI enzyme. Equation (5) represents the deactivation of the Ptet promoter when aTc has been consumed. Finally, equations (7),(8) and (9) represent the degradation of the main chemical species of a aTc, the complex TetR+aTc and the mRNA respectively [1].

Once the LuxI gene has been expressed, the enzymatic reaction of the formation of the quorum sensing molecule (AHL) is taking place. The reaction is considered to be of a first order kinetic since the substrate of the enzyme reaction is produced naturally and abundantly within the cell.

Kinetic Model


To develop the kinetic model, our team worked on two levels. At first, the differential equations for the production of chemical species inside a single cell were developed. Then we scaled-up our approach by studying the production of the lactone as a whole from each subpopulation, as well as its diffusion in the system until it reaches the receivers.

Single cell model

Initially, to simplify the model, the change in the concentration of the basic chemical species in a single cell of each subpopulation was studied. For the creation of the differential equations, both production and consumption as well as the degradation rates were taken into account. The production and consumption rates of each chemical can be represented by the following differential equations:

\(\frac{d[ATC_{int}]}{dt} = k_{-comlex}[complex]- k_{complex}[ATC_{int}][TetR] - d_{ATC}[ATC_{int}] \ (1) \)



\(\frac{d[tetR]}{dt} = \frac{p_{J23105} * k_{tetr} * CN_{senders}}{d_{mRNATetR}} + k_{-complex} [complex] - k_{complex} [ATC_{int}] [TetR] - d_{TetR} [TetR] \ (2) \)



\(\frac{d[complex]}{dt} = k_{complex}[ATC_{int}][TetR] - k_{-complex}[complex] - d_{complex}[complex] \ (3)\)



\(\frac{d[LuxI]}{dt} = \frac{p_{tet} * k_{LuxI} * CN_{senders}}{d_{mRNALuxI}} * [\frac{Kd_{tet}^n}{Kd_{tet}^n + [TetR]^n} - d_{LuxI} [LuxI]] \ (4) \)



\(\frac{d[AHL_{int}]}{dt} = k_{A} [LuxI] - d_{AHL} [AHL] \ (5)\)



At this point it would be useful to mention the different assumptions made for the single cell kinetic model:

1. We assume that the expression of a gene can be fully characterised by a Hill Equation.
2. We assume that the production of the AHL by the LuxI enzyme is a first order kinetic reaction, so the ODE describing this reaction is formed respectively.

Equations (1),(2) and (3) describe the change of the concentration of TetR, the complex TetR+aTc and the LuxI enzyme respectively. Equation (4) is the representation of the first order enzymatic reaction that expresses the change of the concentration of the lactone that is the quorum sensing molecule. The total amount of this lactone simulates the weighted sum of the artificial neural network.

Initially, the kinetic constants for the one-cell model were found in the literature as we needed to analyze the system and have an indicative value for the amount of lactone produced over time and to choose whether the plasmids we would medium or high copy. We followed very specific criteria for the choice of kinetic constants. We always checked if the experiments in the research works from which we were inspired, used the same microorganisms as us and the same experimental conditions as those that we would develop our system. By keeping the desired amount of AHL produced constant, we ran our model and concluded that the copy number of the plasmid used at the senders' system should be low. The the matrix below consists of the initial values of our kinetic constants found in literature:

Kinetic Constant Description Value Source
CN Plasmid Copy Number (Senders) 10(Low) Wet Lab Construct Used
n LuxI Hill Coefficient (in aTc-TetR system) 3 Team:HSiTAIWAN 2016
d_comlex Degradation rate of aTcTetR 0.025(1/min) Literature
d_mRNATetR Degradation Rate of mRNA TetR ln(2)/30 = 0.023(1/sec) Literature
dTetR Degradation Rate of TetR 0.631 [1/sec] Literature
d_AHL Degradation Rate of AHL 0.01 min-1 Literature
d_mRNALuxI Degradation Rate of Luxi Mrna 0.347 (1/min) Literature
d_LuxI Degradation Rate of LuxI protein 0.00167(1/min) Team MIT 2010
d_aTc Degradation Rate of aTc ln(2)/20 (1hours) Literature
KA Synthesis Rate of AHL by LuxI 0.04(1/min) Literature
k-1_aTcTetR Dissociation Rate of aTcTetR 4.2*10^-4(1/min) Literature
k_aTcTetR Synthesis Rate of aTcTetR k-1/Kd
Kd_pTet Dissociation constant of tetR with the pTet Promoter 50 nM Team:HSiTAIWAN 2016
Kd_aTcTetR Dissociation constant of aTC with the tetR molecule 15 nM Team:HSiTAIWAN 2016
p_pTet Transcription Rate of ptet promoter 25.3*10^(-3) [1/sec] Calculated using "Salis Lab" Calculator
p_J23105 Transcription Rate of J23105 promoter 30*10^(-3) [1/sec] Registry
k_luxi Translation Rate of LuxI Taken from our RBS Library Calculated using "Salis Lab" Calculator
k_tetr Translation Rate of TetR 11*10^(-3)[1/sec] Calculated using "Salis0.129*10 Lab" Calculator
Nmin Initial value of E.coli 0.129*10^(7) Calculated from lab data about OD
Nmax Maximum value of E.coli 2.14*10^(7) Calculated from lab data about OD
Vbead Volume of culture 200 μL From Wet Lab Setup
Table 2: Kinetic constants obtained from literature.


In terms of the non-literature values, those are obtained by the In Silico Tool of Salislab that offered a proportional scale of strength for our promoter and rbs constructs.

Time Scale Convertion


Translation Rates

The maximum possible translation rate offered by salislab is around 250000 a.u. which corresponds to the maximum elongation rate, i.e 18 aminoacids/sec. Knowing that each amino acid has 3nucleotides, the following time-scale formula emerges:

k_gene=(4.14 proteins/mRNA/second)*[TIR_salis/250000]

Transcription Rates

In a similar manner, the maximum transcription elongation rate in E.Coli is 90nts/sec and the RNA polymerase II footprint was found to be 40 nucleotides [4]. Finding out that the maximum transcription rate of Salislab is 100000 a.u. ,we have the formula:

\( p_{gene}=(90/40)*(TIR_{salis}/100000) \), or

\( p_{gene}=2.25*(TIR_{salis}/100000) \)

The above formulas were used in order to receive the time-scale values of our constructs.

Whole cell model

Then, wanting to have a complete picture of our overall system at the level of bacterial populations, we upgraded our model for multiple cells. The scale-up had to include the study of phenomena such as the growth of the bacterial population and how this affects the amount of lactone produced as well as the diffusion of the lactone in and out of the cells. The assumptions made for the whole cell kinetic model are the following:

1. We assume that the expression of a gene can be fully characterised by a Hill Equation.
2. We assume that the production of the AHL by the LuxI enzyme is a first order kinetic reaction, so the ODE describing this reaction is formed respectively.
3. We assume that the diffusion of the lactone outside the cells is characterized by the same Fick's constant as the diffusion of the lactone inside the cells.
4. We assume that the growth of the bacterial population occurs logarithmically.

The production and consumption rates of each chemical can be represented by the following added differential equations:

\(\frac{dN}{dt} = \mu * [N] *(1-\frac {[N]}{[N_{max}]}) \ (1) \)



\(\frac{d[ATC_{ext}]}{dt} = D_{ATC} [N] ([ATC_{int}] - \frac{V_{E.Coli}}{V_{bead}} [ATC_{ext}]) - d_{ATC} [ATC_{ext}] \ (2) \)



\(\frac{d[ATC_{int}]}{dt} = D_{ATC} * (\frac{V_{E.Coli}}{V_{bead}} [ATC_{ext}] - [ATC_{int}]) + k_(-complex)[comlex] -k_{complex} [ATC_{int}][TetR] -d_{ATC} [ATC_{int}] \ (3) \)



\(\frac{d[tetR]}{dt} = \frac{p_{J23105} * k_{tetr} * CN_{senders}}{d_{mRNATetR}} + k_{-complex} [complex] - k_{complex} [ATC_{int}] [TetR] - d_{TetR} [TetR] \ (4) \)



\(\frac{d[complex]}{dt} = k_{complex}[ATC_{int}][TetR] - k_{-complex} [complex] - d_{complex} [complex] \ (5)\)



\(\frac{d[LuxI]}{dt} = \frac{p_{tet} * k_{LuxI} * CN_{senders}}{d_{mRNALuxI}} * (\frac{Kd_{tet}^n}{Kd_{tet}^n + [TetR]^n} - d_{LuxI} [LuxI]) \ (6) \)



\(\frac{d[AHL]}{dt} =k_{A} [LuxI] - d_{AHL} [AHL] \ (7) \)



\(\frac{d[AHL_{int}}{dt} = k_{A} [LuxI] + D_{AHL} ( \frac{V_{E.Coli}}{V_{bead}} [AHL_{ext}]) - d_{AHL} [AHL_{ext}] \ (8) \)



\(\frac{d[AHL_{ext}}{dt} = D_{AHL} [N] ( [AHL_{int}] \frac{V_{E.Coli}}{V_{bead}} [AHL_{ext}]) - d_{AHL} [AHL_{ext}] \ (9) \)



Equation (1) is a representation of how the number of cells is growing with time. Equations (2) and (3) explain how the concentration of the aTc inducer is changing inside and outside of a cell due to diffusion, and similarly equations (8) and (9) describe the diffusion phenomena for AHL-OC6.

Model Development

In order to develop our model, we used python programming language. More specifically, we used the the libraries scipy and lmfit. lmfit provides a Least-Squares Minimization routine and class with a simple, flexible approach to parameterizing a model for fitting to data. More information on the senders model as well as the results of our experiments can be found in our GitLab repository here.

Senders Model Code


Results


As mentioned above, the initial purpose of our kinetic model-both for the single cell and the whole cell versions-was to understand our system and its various parts, not only in the time domain, but also its behavior based on the different aTc inputs, even with a limited amount of experimental data. To that end initially, we performed some experiments based, for which the parmaeters of the system we tuned based on the current literature. In this way, we obtained some preliminary results on the systems behavior that were shared and discussed with the wet lab. More information on the Engineering Success page.

It should be noticed that due to the fact that the TetR anc kd_tet values are of similar scale, the repressible model script couldn't produce the expected results, especially when we tried to fit parameters via the lmfit() library. As a result, in terms of coding, we used an inducible system for the production of LuxI by the TetR-aTc complex, since it didn't lead to misfit errors. The results presented below refere to our kinetic model of the BCD2 RBS variant, at full induction (aTc=10**(-4) M ).

Single-Cell Model


Figure 2: Change in concentration of chemical species after solving the ODs for the single-cell model.

Whole-Cell Model


Figure 3: Change in concentration of chemical species after solving the ODs for the whole-cell model.

Discussion

From the graphs above, we can comprehend the following:
1. The concentration of external and internal aTc follows similar values, due to a significant decrease of volume from the external volume of the bead (at 200μL) to the Vecoli (at around 1000pL). This was also experimentally validated by [5] where a maximum 10-fold gap between extracellular and intracellular aTc was observed (see Supplementary Table 10).
2. The diffusion of aTc is immediate leading to an instant reduction of the TetR Repressor. However, at around the 3h mark the entire aTc concentration tends to zero, which subsequently brings back the TetR repressor and diminishes the amount of AHL produced.
3. In the whole cell model, we observe a decrease of the internal AHL molecules due to the fast diffusion phenomena throughout the E.coli membrane (described by the D_AHL constant). The entire AHL amount of the cell culture is estimated at 5*10**(11) molecules (for a 200μL Bead).

In terms of the transfer function of the senders, we observed a clear difference between the single-cell and the whole cell models:

Figure 4: Ρepresentation of the transfer function for single-cell and whole-cell model (RBS: BCD2)


As we observe, the lack of diffusion phenomena in the first modelling approach lead to a false estimate of the aTc value required for full induction of the senders population. As a result, our team decided to keep the whole cell model for parameter fitting and simulations, since it would correlate better with the complex conditions of our wets lab experiments.

Lab Data Comparison and Model Fit


At first, a straight-comparison of our literature whole-cell model with the lab induction data was required in order to compare the two models, detect the differences and fault estimates.

Figure 5: Comparison of the theoretical model with the experimental data (RBS: BCD2)


Despite its ability to predict the experimental value of LuxI measured at full induction, our literature model lacked the capability to accurately predict the initiation of induction. Due to that fact, a model fitting method-via error minimization-was required in order to obtain the correct model parameters. As metnioned previously, we used the python library lmfit and more specifically the routine minimize() with the integrated Nelder-Mead slgorithm in order to perform our parameter optimization.

Nelder-Mead algorithm is a numerical method used to find the minimum or maximum of an objective function in a multidimensional space. It is a direct search method (based on function comparison) and is often applied to nonlinear optimization problems for which derivatives may not be known. A simple flow-chart of the Nelder-Mead algorithm is presented below:

Figure 6: Flow-chart of the Nelder-Mead algorithm. [6]


On the table below, we show the parameters opted to be fitted to our experimental data as well as the fitting results for the BCD 2 RBS. For each parameter you can find: the parameter symbol, description, the value obtained from literature review and its source and finally the calculated value, based on the experimental data fit.:

Symbol Description Literature Value Source Fitted Value
mi Senders Dilution Rate 0.019 [1/min] Harvard Bionumbers 0.032 [1/min]
a_tetr Production Rate of TetR k_tetr*p_ptet=1.18[1/min] Calculated with "SalisLab" Calculator 10.22
b_ptet Leakiness of Ptet promoter 0 (%) assumed 20(%)
a_luxi Production Rate of LuxI (for BCD2) k_luxi*p_luxi=312[1/min] Calculated with "SalisLab" Calculator 128
n LuxI Hill Coefficient 3 Team:HSiTAIWAN 2016 2
kd_ptet Dissociation Constant of TetR with the ptet promoter 50nM Team:HSiTAIWAN 2016 2.5nM
Table 3: Kinetic constants overview.


Figure 7: Model fitting for senders parameters for RBS BCD2.


As noted in Table 3, the theoretical model that was created based on the parametes we obtained by literature review, deviates from the model that was created based on the fitting of the experimental data. More specifically, there is a big deviation between the known and predicted dissociation constant. The dissociation constant is a crucial value for our system since it determines the input aTc on which the induction will begin and on our case the literature value was a bit too high for our own system. Moreover, the parameters that were based on the proportional values for the RBSs obtained from the Salis RBS caluclator seem to not agree with our experimental data. Finally, the binding affinity is also a variable dependent on the experimental setup and measurements. Binding affinity corresponds to do with the steepness of the plotted curve.

The whole-cell fitted model of the senders was used for the implementaion of the pattern recognition (see more in the Pattern Recognition section), where the communication of the senders and the receivers cultures is simulated.

Receivers System Explained



The receivers population simulates the activation function of the perceptron algorithm. So it was really important for us to create constructs in a way that we would secure a steep activation of the production of the final fluorescent protein, mNeon Green. We had to design plasmids that would give the receivers cells the ability to express our protein very strongly only if the amount of lactone they receive from the senders exceeds a certain threshold. We decided to test two different constructs in order to check which one will present the best response.


Receivers' Constructs


Construct no.1: Open Loop (OpLo)

The open loop construct is the simplest. It consists of a Plux promoter followed by the luxR gene and a Plux promoter followed by the mNeonGreen gene. Basically, this is a system where LuxR feeds the Plux promoter to express both itself and the fluorescent protein. Thus, when a specific threshold of OC6 is exceeded, the production of the LuxR increases sharply, which implies a very fast production of mNeonGreen. In this way it succeeds the steep activation.

Figure 1: Representation of Open Loop construct.

Construct no.2: PFR Loop
The PFR construct aims to increase the steepness of our activation function via the following process:
The LuxR is a leaky promoter and thus its leakiness provides a baseline production of LuxR that after the addition of OC6 will activate the production of more LuxR and thus LuxR-AHL. The LuxR-AHL dimer produced will then induce the production of mNeonGreen (as in the OpLo case) and it will simultaneously repress the production of phlf, via a LuxR repressible promoter. The phlf is responsible for the repression of the mNeonGreen construct, so its repression , along with the induction by LuxR-AHL, will lead to a simultaneous increase in the output concentration. The simultaneous control of mNeonGreen by the LuxR-AHL dimer and the phlf repressor is accomplished via a pLux-phlo hybrid promoter. This produces a steeper activation function than the Open Loop circuit, as it is also noted in the results section.

Figure 2: Representation of PFR Loop construct.

Biochemical Reactions


For the kinetic modelling of the receivers system, we had to convert all the aforementioned design choices in the constructs into biochemical equations, each of which is characterized by its own rate of consumption or production of the chemical species. The resulting equations are:

1. ∅ →LuxRreg
2. LuxAHL + LuxRreg ↔LuxRcomplex
3. LuxRcomlex + LuxRcomplex ↔ DLuxR
4. DLuxR + P LuxRof f ↔P LuxRon
5. P LuxRon → P LuxRon +mRNAmNeon
6. mRNAmN eon→mN eon
7. LuxAHL→ ∅
8. LuxRreg→ ∅
9. LuxRcomplex→ ∅
10. DRLux→ ∅
11. mRN AmN eon→ ∅
12. mN eon→ ∅


Equation (1) describes the production of the LuxR regulator. Equation (2) describes the formation of the AHL and the regulator complex when aTc is present in the system and this formation leads to equation (3) describing the formation of a dimer. This dimer leads to the activation of the LuxR promoter in equation (4). Equations (5) and (6) describe the translation of the gene and finally the production of the LuxI enzyme. Finally, equations (7),(8),(9),(10),(11) and (12) represent the degradation of the main chemical species.

Kinetic Model


To develop the kinetic model of the receivers, our team worked on two levels as for the senders' kinetic model. At first, the differential equations for the production of chemical species inside a single cell were developed. Then we scaled-up our approach by studying the production of the fluorescent protein from the whole receivers population. The difference lies in the fact that for this model our team analyzed and also tested three constructs experimentally. The purpose of this effort was to see which of the three systems we had designed would work best as an activation function of a perceptron and would have the sharpest response. You can find more information in the Engineering Success page. Overall, we created one single-cell and one whole-cell model for each construct.

Construct no.1: Open Loop


Single cell model

Initially, to simplify the model, the change in the concentration of the basic chemical species in a single cell of each subpopulation was studied. For the creation the differential equations, both production and consumption as well as the degradation rates were taken into account The production and consumption rates of each chemical can be represented by the following differential equations:

\(\frac{d[AHL_{int}]}{dt} = k_{-1}*[Mon]+D_{AHL}*(\frac{V_{eColi}}{V_{bead}} *AHL_{ext} - AHL_{int}) -k_{1}*[AHL_{int}][LuxR] - d_{AHL} [AHL_{int}] \ (1)\)



\(\frac{d[LuxR]}{dt}= \frac{p_{J23105}*k_{LuxR}*CN_{rec}}{d_{mRNALuxR}} + k_{-1} [Mon] - k_{1}[LuxR][AHL_{int}] - d_{LuxR}[LuxR] \ (2)\)



\(\frac{Mon}{dt} = k{1}[LuxR][AHL_{int}] + 2k_{-2}[Dimer]-k_{-1}[Mon]-2k_{2}[Mon]^2 - d_{Mon}[Mon] \ (3)\)



\(\frac{Dimer}{dt}= k_{2}[Mon]^2 - k_{-2}[Dimer]-d_{Dimer}[Dimer] \ (4) \)



\(\frac{d_{mng}}{dt} = \frac{p_{Lux}*k_{mng}*CN_{rec}}{d_{mRNAmng}} * [\beta_{Lux}+(1-\beta_{Lux})*\frac{[Dimer]^n}{Kd_{Lux}^n+[Dimer]^n}] -d_{mng}[mng] \ (5) \)



At this point it would be useful to mention the different assumptions made for the single cell kinetic model:

1. We assume that the expression of a gene can be fully characterised by a Hill Equation.
2. The LuxR-AHL dimerization is described via a second-order kinetic reaction, so the ODE describing this reaction is formed respectively.
3. The degradation rate of AHL inside and outside the membrane remains the same.


Equation (1) represents the change of concentration of the AHL, in order to complete the induction inside the cell. Equation (2) describes the change of the concentration of the LuxR from the lactone reaching the receiver's cell. Equations (3) and (4) describe the change of the concentration of the monomer complex and the dimer, practically leading to the activation of the promoter. Equation (5) is the representation of the production of the fluorescent protein, the mNeonGreen.

Whole cell model

In order to have a complete picture of our overall system at the level of bacterial populations, we upgraded our model for multiple cells. The scale-up had to include the study of phenomena such as the growth of the bacterial population and how this affects the amount of lactone produced, as well as the diffusion of the lactone in and out of the cells.

The assumptions made for the whole cell kinetic model are the following:

1. We assume that the expression of a gene can be fully characterised by a Hill Equation.
2. We assume that the diffusion of the lactone outside the cells is characterized by the same Fick’s constant as the diffusion of the lactone inside the cells.
3. We assume that the growth of the bacterial population is occurs logarithmically.

The production and consumption rates of each chemical can be represented by the following added differential equations:

\(\frac{dN}{dt} = \mu * [N] *(1-\frac {[N]}{[N_{max}]}) \ (1)\)



\(\frac{d_{AHL_{ext}}}{dt} = D_{AHL}*[N]*(AHL_{int}- \frac{V_{eColi}}{V_{bead}} *AHL_{ext}) - d_{AHL}[AHL_{ext}] \ (2)\)



\(\frac{d[AHL_{int}]}{dt} = k_{-1}*[Mon]+D_{AHL}*(\frac{V_{eColi}}{V_{bead}} *AHL_{ext} - AHL_{int}) -k_{1}*[AHL_{int}][LuxR] - d_{AHL} [AHL_{int}] \ (3) \)



\(\frac{d[LuxR]}{dt}= \frac{p_{J23105}*k_{LuxR}*CN_{rec}}{d_{mRNALuxR}} + k_{-1} [Mon] - k_{1}[LuxR][AHL_{int}] - d_{LuxR}[LuxR] \ (4)\)



\(\frac{d[Mon]}{dt}= k_{1} [LuxR][AHL_{int}]+2 k_{-2}[Dimer]-k_{-1}[Mon]^2 - d_{Mon}[Mon] \ (5)\)



\(\frac{d[Dimer]}{dt} = k_{2}[Mon]^2 - k_{-2} [Dimer] - d_{Dimer} [Dimer] \ (6) \)



\(\frac{d[mng]}{dt} = \frac{p_{Lux}*k_{mng}*CN_{rec}}{d_{mRNAmng}} + [\beta_{Lux} +(1-\beta_{Lux})* \frac{[Dimer]^n}{Kd_{Lux}^n+[Dimer]^n}] - d_{mng} [mng] \ (7)\)



Equation (1) is a representation of how the number of cells is growing with time. Equations (2) and (3) explain how the concentration of the lactone AHL is changing inside and outside of a cell due to diffusion.

Construct no.2: PFR Loop


Single cell model

> In a similar manner, we create the single and whole cell model of the PFR construct:

\(\frac{d[AHL_{int}]}{dt} = k_{-1}[Mon] - k_{1} * [AHL_{int}][LuxR] - d_{AHL}[AHL_{int}] \ (1) \)



\(\frac{d[LuxR]}{dt} = \frac{p_{Lux}*k_{LuxR}*CN_{rec}}{d_{mRNALuxR}} * [(1-\beta_{Lux})+\beta_{Lux} *\frac{[Dimer]^n}{Kd_{Lux}^n+[Dimer]^n}] + k_{-1}[Mon]-k_{1}[LuxR][AHL_{int}] - d_{LuxR}[LuxR] \ (2) \)



\(\frac{d[Mon]}{dt} = k_{2} [Mon]^2 - k_{-2}[Dimer] - k_{-1}[Mon]-2k_{2}[Mon]^2 -d_{Mon}[Mon] \ (3) \)



\(\frac{d[Dimer]}{dt} = k_{2} [Mon]^2 - k_{-2}[Dimer] - d_{Dimer}[Dimer] \ (4)\)



\(\frac{d[phlf]}{dt} = \frac{p_{Lux}*k_{phlf}*CN_{rec}}{d_{mRNAphlf}} * [\beta_{Lux} + (1-\beta_{Lux}) * \frac{Kd_{Lux}^n}{Kd_{Lux}^n+[Dimer]^n} ]- d_{phlf}[phlf] \ (5)\)



\(Hill_{Lux}= \beta_{Lux} + (1-\beta_{Lux}) *\frac{[Dimer]^n}{Kd_{Lux}^n+[Dimer]^n} \)



\(Hill_{phlf}= \beta_{phlf} + (1-\beta_{Lux}) *\frac{Kd_{Lux}^n}{Kd_{Lux}^n+[Dimer]^n} \)



\(\frac{d[mng]}{dt} = \frac{p_{hybrid}*k_{mng}*CN_{rec}}{d_{mRNAmng}} *Hill_{Lux} * Hill_{phlf} - d_{mng} [mng] \ (6)\)



Equations (1) through (4) remain the same as in the Open Loop System and describe the creation of the dimer from LuxR and AHL. Equation (5) has to do with the repression of the phlf via the LuxR-AHL dimer. Equation (6) describes the induction of mneongreen via the hybrid promoter. The hybrid promoter has been modelled as a multiplication of two separate hill equations, one for the induction of the pLux promoter and one for the phlo-box,where phlf is being hitched in order to cause the repression of mneon.

Whole-cell model:

The diffusion phenomena of the OC6 through the cell membrane have been described via the equations (1)-(3) and they are exactly the same as in the open loop.

\(\frac{dN}{dt} = \mu * [N] *(1-\frac{[N]}{[N_{max}]}) \ (1)\)



\(\frac{d[AHL_{ext}]}{dt} = D_{AHL} *[N] *( AHL_{int} - \frac{V_{eColi}}{V_{bead}} *AHL_{ext}) - d_{AHL}[AHL_{ext}] \ (2)\)



\(\frac{d[AHL_{int}]}{dt} = k_{-1}*[Mon]+D_{AHL}*(\frac{V_{eColi}}{V_{bead}}*AHL_{ext} - AHL_{int}) -k_{1}*[AHL_{int}][LuxR] - d_{AHL} [AHL_{int}] \ (3) \)



\(\frac{d[LuxR]}{dt}= \frac{p_{Lux}*k_{LuxR}*CN_{rec}}{d_{mRNALuxR}} [\beta_{Lux} + (1-\beta_{Lux}) *\frac{[Dimer]^n}{Kd_{Lux}^n+[Dimer]^n}] + k_{-1} [Mon] - k_{1}[LuxR][AHL_{int}] - d_{LuxR}[LuxR] \ (4)\)



\(\frac{d[Mon]}{dt} = k_{1}[LuxR][AHL_{int}] + 2 k_{-2}[Dimer] - k_{-1} [Mon] -2k_{2}[Mon]^2 - d_{Mon} [Mon] \ (5)\)



\(\frac{d[Dimer]}{dt} = k_{2} [Mon]^2 - k_{-2} [Dimer] \ (6)\)



\(\frac{d[phlf]}{dt} = \frac{p_{Lux}*k_{phlf}*CN_{rec}}{d_{mRNAphlf}} * [\beta_{Lux} + (1-\beta_{Lux}) * \frac{Kd_{Lux}^n}{Kd_{Lux}^n+[Dimer]^n} ]- d_{phlf}[phlf] \ (7)\)



\(Hill_{Lux}= \beta_{Lux} + (1-\beta_{Lux}) *\frac{[Dimer]^n}{Kd_{Lux}^n+[Dimer]^n} \)



\(Hill_{phlf}= \beta_{phlf} + (1-\beta_{Lux}) *\frac{Kd_{Lux}^n}{Kd_{Lux}^n+[Dimer]^n} \)



\(\frac{d[mng]}{dt} = \frac{p_{hybrid}*k_{mng}*CN_{rec}}{d_{mRNAmng}} *Hill_{Lux} * Hill_{phlf} - d_{mng} [mng] \ (8)\)



Tables 1 and 2 below contain the values of our kinetic constants obtained by reviewing the current literature for the OpLo and PFR constructs respectively:

Kinetic Constant Description Value Source
CN Plasmid Copy Number (Receivers) 17(Medium) Wet Lab Constructs
n Hill Coefficient 1.5 Estimated (before fit)
k1 Synthesis Rate of LuxRAHL k-1/kd1
k2 Synthesis Rate of LuxRAHL2 k-2/kd2
k-1 Dissociation Rate of Monomer (LuxR-AHL) 10(1/min) Literature
k-2 Dissociation Rate of of dimer (LuxR-AHL)2 1(1/min) Literature
kd2 Dissociation Constant of Dimer (LuxR-AHL)2 20 Nm Literature
kd1 Dissociation Constant of Monomer (LuxR-AHL) 100 nM Literature
kdlux Dissociation constant of (LuxR-AHL)2 to the plux promoter 200nM Literature
d_RA2 Degradation Rate of LuxRAHL2 0.017(1/min) Team Valencia 2018
d_RA Degradation Rate of LuxRAHL 0.156(1/min) Team Valencia 2018
d_Mng Degradation Rate of Mng 0.010(1/min) Literature
d_mRNAmNGt Degradation Rate of Mrna_Mng 0.039(1/min) Literature
d_LuxR Degradation Rate of LuxR 0.002(1/min) Literature
d_mRNALuxR Degradation Rate of LuxR Mrna 0.347(1/min) Literature
β_Lux Basal Expression of pLux Promoter (as min/max) 4.8%=(1043/21290) Registry
D_AHL Diffusion Rate of External AHL through the membrane 10(1/min) Literature
V_Ecoli E.Coli Cell Volume 1μ (m**3)=10**(-15) L Literature
μ dilution rate 0.019(1/min) Literature
Transc_Plux Transcription Rate of LuxR promoter 0.79/60(1/sec) Team Valencia 2018
Vbead Volume of culture 200 μL From Wet Lab Setup
Nmin Minimum value of cells in Bead 0.29*10**(7) From OD experimental Data
Nmax Maximum value of cells in Bead 2.142*10^(7) From OD experimental Data
d_AHL Degradation Rate of AHL 0.01 min^-1 Literature
p_J23105 Transcription Rate of constitutive promoter for LuxR 58*10^(-3) (1/sec) Registry
p_Lux Transcription Rate of LuxR Promoter for mneongreen 0.79(1/min) "Salis Lab" Calculator
k_luxr Translation Rate of LuxR (RBS BBa_B0034) 18*10^(-3) (1/sec) "Salis Lab" Calculator
k_mneongreen Translation Rate of mneongreen(RBS BBa_B0034) 27*10^(-3) (1/sec) "Salis Lab" Calculator
Table 1: Kinetic constants found in literature for the Open Loop construct.


Kinetic Constant Description Value Source
p_lux Transcription Rate of LuxR promoter for LuxR 155*60*10^(-3) (1/sec) "Salis Lab" Calculator
k_luxr Translation Rate of LuxR (RBS BBa_B0034) 27*10^(-3) (1/sec) "Salis Lab" Calculator
p_lux_repr Transcription Rate of LuxR_Repressible promoter for phlf 3.69*10^(-3) (1/sec) "Salis Lab" Calculator
k_phlf Translation Rate of Phlf (RBS BBa_B0034) 4*10^(-3) (1/sec) "Salis Lab" Calculator
p_hybrid Transcription Rate of Hybrid Promoter for mneongreen 10*10^(-3) (1/sec) "Salis Lab" Calculator
k_mng Translation Rate of mneongreen (RBS BBa_B0034) 18.3*10^(-3) (1/sec) "Salis Lab" Calculator
d_mrnaphlf Degradation Rate of mrna phlf 0.02 (1/sec) "Salis Lab" Calculator
d_mrnaphlf Degradation Rate of mrna phlf 0.02 (1/sec) Team SUSTech_Shenzhen 2021
d_phlf Degradation Rate of phlf 0.042(1/min) Team Hong_Kong_HKUST 2017
Kd_phlf Dissociation Constant for the phlf repressible promoter 2*10^(-7) nM Team Hong_Kong_HKUST 2017
Table 2 : Additional kinetic constants found in literature for the PFR Loop construct.


OpLo Model Code



PFR Model Code


Results


Similar to pipeline we followed for modelling our senders system, our team also compared the theoretical single and whole cell models for the receivers, both in time-scale and for different induction concentrations of AHL. However, in this section, we are going to present only a part of our results, since the full analysis and codes, with their respective results, are available on our GitLab repository, here.

Below we present the time-scale results of our whole-cell PFR model, at full induction of 10^(-6) M OC6.

Figure 2: Change in concentration of chemical species after solving the ODs.

From the graphs above we can clearly notice that, as in the case of aTc, the diffusion of AHL through the membrane is rapid, and it leads to a similar concentration of internal and external AHL [8]. This leads to the creation of the dimer inducer within the 100 minutes, which not only induces the production of mNeonGreen, but also reduces the repressor phlf. This observation lasts until the 600 minute mark, where-since there is no AHL inside the cell-the dimer starts significantly decreasing, thus stopping the repression of phlf. The increase of phlf and reduction of dimer mean that the produced fluorescence starts dropping (around the 600 minute mark).

Similar results were shown by the rest single-and-whole cell models of OpLo and PFR.

Transfer Function of OpLo


First, we represent the literature transfer function both for the single and the whole cell models.

Figure 3: Transfer function for single-cell and whole-cell model.

Similar to the senders model, the whole-cell description of our system gives more consistent results and is closer to the experimental data, since the induction in that case starts for lower concentrations of AHL. A similar difference was observed in the PFR single-and-whole cell models as well.

Comparison to Experimental Data and Model Fit of OpLo


First, we aimed to compare our theoretical whole-cell model to the experimental data. By reviewing the results below, we assume that-despite its approximately valid representation of our system-our literature model lacks the accuracy needed in order to be used for the pattern recognition task. Thus we again decided to perform model fit using the python library lmfit.

Figure 4: Comparison of the theoretical model with the experimental data for the Open Loop construct.

Symbol Description Initial Value Source Fitted Value
n Hill coefficient for OpLo 1.5 estimated 0.5
b_plux Leakiness of pLux promoter 4.9% Registry-BBa_R0062 10%
a_luxr Production Rate of LuxR p_J23105*k_luxr=5.68 SalisLab 49.78
a_mng Production Rate of mneongreen p_plux*k_mng=0.85 Salis Lab 0.44
kd_lux Dissociation Constant of LuxR Promoter 200nM Literature 10nM
d_mng Degradation Rate of mneongreen 0.01(1/min) Literature 0.023(1/min)
d_mrna_mng Degradation Rate of mrna mneongreen 0.039(1/min) Literature 0.019(1/min)
mi Dilution Rate 0.019(1/min) Literature 0.0028(1/min)
Table 3: Fitted values for Open Loop construct.

Below the fitted model for OpLo is presented:

Figure 5: Model fitting for Open Loop construct

The model fit process gave us some important variables, which could not be accurately estimated otherwise. These include: the curve steepness, the production rates of LuxR and mNeonGreen as well as the dissociation constant of the PLux promoter. Some of those values were used for the characterization of the PFR construct below.

PFR: Transfer Function and model fit


Using the fitted values of d_mng,d_mrnamng, kd_lux and b_plux that would remain unchanged for the PFR and with the rest of the values from literature, we extracted the same comparison of theoretical and and experimental transfer functions:

Figure 6: Comparison of the theoretical model with the experimental data for the PFR loop construct.


As we can clearly observe, the OpLo results used for the characterization of PFR resulted in a more accurate prediction of the later construct. However, there are still some inaccuracies, in terms of the curve steepness,the leakiness of the hybrid promoter and the characterization of the phlf molecule. After the fitting, we obtained a more accurate reconstruction of the experimental data via our model:

Figure 7: Model fitting for PFR Loop construct.


Symbol Description Initial Value Source Fitted Value
a1 Production Rate of LuxR k_luxr*p_lux=15 Salis Lab 5.67
a2 Production Rate of Phlf k_phlf*p_lux_repr=0.053 Salis Lab 0.1
a3 Production Rate of mneongreen p_hybrid*k_mng=0.6588 Salis Lab 0.62
n Hill Coefficient for PFR 2 Estimated 2.31
b_phlf Leakiness of phlf promoter 10% Estimated 0.36
kd_phlf Dissociation Constant of phlf with phlf promoter 2*10^(-7)) Team Hong_Kong_HKUST 2017 7.95*10^(-7)
d_phlf Degradation rate of phlf 0.042(1/min) Team Hong_Kong_HKUST 2017 0.03
d_mrna_phlf Degradation rate of mrna_phld 1.2(1/min) Team SUSTech_Shenzhen 2021 0.06
Table 4: Fitted values for PFR Loop construct.


After the successful characterization of the transfer function circuits, we aimed to directly compare them while simulating a pattern recognition task and decide which implementation is more accurate in terms of expected and actual output concentrations.

Stability Analysis


As it has been analysed in the Impementation page, our team is inspired to utilize its biological perceptron for recognizing complex patterns in a variety of different circumstances, leading to the creation of a novel “thinking” biosensor of E.Coli consortia. However, it is of paramount importance to gain in-depth knowledge about more complex aspects of our system, apart from the concentration analysis that our main model offers. One such aspect are the diffusion phenomena through the cell membrane, that try to demonstrate the change of AHL concentration as the molecules get from the senders to the receivers consortia and the various effects on both populations.On the other hand, the stability analysis that will follow offers context on a different issue regarding our biosensor; the frequency that it could receive new input values for classification.

Below we summarize the fundamental theory from control systems. At first we should validate that our system is stable and will lead to a steady state concentration of mNeonGreen. For that purpose and for the entire analysis, our initial model equations have been simplified to just two equations, one for the senders and one for the receivers, as shown below:

\( a_{send} = a_{1} + a_{2} +a_{3}\)



\(\frac{d[x_{1}]}{dt} = f_{1}(t)= (a_{send}) * \frac{u^{n_{1}}}{K_{1}{^n_{1}}+u^{n_{1}}} - d_{1}* x_{1}\)



\(\frac{d[x_{2}]}{dt} = f_{2}(t)= (a_{rec}) * \frac{x_{1}{^n_{2}}}{K_{2}{^n_{2}} + x_{1}{^n_{2}}} - d_{2} * x_{2} \)



,where the three senders subpopulations are being modelled as one, which produces the summation of the subpopulations OC6. This OC6 value, shown as x1 here, is used as induction input to the receivers population, whose fluorescence is shown as the x2 variable. Finally, both equations are using Hill functions to model the induction.

However, in order to properly study the system's frequency behavior, we should first perform a linearization method which brings our system to the state-space. After performing such a task, we received the following linearized control model:

\( x_{dot}(t) = A * x(t) + B* u(t) \)

\(y_{dot}(t) = C * x(t) +D * u(t)\)

\(x(t) = [x_{1}(t),x_{2}(t)]\)

\(u(t)=aTc\)

\(x_{1}(t)= AHL\)

\(x_{2}(t) = out\)


\( where \ A=[a_{11}, a_{12}][a_{21}, a_{22}] \)


\(B=[b_{1}],b_{2}]\)



\(a_{11}= \frac{df_{1}}{dx_{1}} = -d_{1} \)



\(\frac{d[x_{2}]}{dt} = f_{2}(t)= (a_{send}) * \frac{u^n2}{K_{1}^n2+u^n2} - d_{2}* x_{2}\)



\(b_{1}=\frac{df_{1}}{du}=a_{send} = \frac{a_{send}*n_{1}*K_{1}^{n_{1}} * u(t_{0})^{n_{1}-1}}{(u(t_{0})^n_{1} + K_{1}^{n_{1}})^2}\)



\(a_{21}=\frac{df_{2}}{dx_{1}}=a_{rec} = \frac{a_{rec}*n_{2}*K_{2}^{n_{2}} * x_{1}(t_{0})^{n_{2}-1}}{(x_{1}(t_{0})^{n_{2}} + K_{2}^{n_{2}})^2} \)



\(a_{22} = \frac{df_{2}}{dx_{2}}=-d_{2}\)



\(b_{2} = \frac{df_{2}}{du}=0\)



It should be assumed that the initial value of u and x1, i.e. u(t0) and x1(t0) should be different to zero, in order to have non-zero values and thus a meaningful analysis. Considering from our diffusion analysis that the diffusion is a rapid procedure that begins around the 6 second mark, it would be safe to assume that we would achieve non zero values for t0>1minute.

In terms of the stability of our system, it is common truth [9] that a stable system has all of its eigenvalues negative, which is true in our case since the λI-A=0 equation brings the following solutions:

\(λ_{1} = -d_{1} \ and \ λ_{2} = -d_{2}\)


,which are negative, assuming positive degradation rates of aTc and AHL.

Moving on, we quantified how often we could update the input states of our perceptron without leading to any instabilities (i.e. problematic output). And there also comes control theory to offer us a useful tool, the Bode plots and the meaning of phase margin. In particular, phase margin is a metric that corresponds to the difference between the phase lag of our system at a given frequency minus the value of -180degrees. This ia a crucial value for our system since phase values below -180degrees lead to unstable systems. In simpler terms, the higher the phase margin for a certain frequency the better for our system.

Furthermore, it is a common ground that, when designing electronic systems and amplifiers, we set a minimum value of phase margin (usually using the frequency where we have unity gain), and use a variety of theoretical and practical methods in order to accomplish that. More information on the electrical side offers the 6th chapter of the book [10]. However, in our case it is just enough to set as maximum frequency f_max the frequency where a minimum amount of phase margin is accomplished, which will be set at 35 degrees and corresponds to a phase lag of -145degrees. Our system should not receive input values at frequencies that reduce the phase margin below that threshold.

Initially, in order to receive the bode plots of our system, we converted the state-space system to a transfer function H(s), which was accomplished via matlab's module ss2tf(). We received:

\(H(s)= \frac{n1*n2*K_{1}^{n1}*K_{2}^{n2}*a_{rec}*a_{send}*x1(t0)^{n1}*u(t0)^{n1-1}}{x1(t0)*(s+d1)*(s+d2)} = \frac{A}{x1(t0)*(s+d1)*(s+d2)}\)



where A and x1(t0) are constant values, for known t0 initial time of analysis. Below the parameter table for stability analysis constants is presented.
Symbol Value
K1 50nM
K2 200nM
n1 2
n2 1.5
d1 0.01 [1/min]
d2 0.01 [1/min]
x1(t0) 2*10^(-9)
u(t0) 10^(-6)
a_rec 0.85
a_send 312
Table 1. Parameter Tables for stability analysis transfer function with the OpLo construct parameters. The parameters refer to the BCD2 literature construct.


Since the gain A/x1(t0) of our system is not the main purpose of it, it was not our main concern. On the other hand, we have to take a closer look at the denominator of the transfer function, and the values d1,d2 which are the poles of the system. Generally, a pole leads to a decrease of the phase lag of our system by 90 degrees. So in our case, the two poles of d1 and d2 will eventually lead to an unstable system after a certain frequency, which we would like to be the highest possible. The problem of decrease in the phase value is more serious in the case of two poles which are close to each other, since this leads to a more rapid decrease of our system's phase lag, almost instantly to -180 degrees.

For our system, we noticed that d1=d2=0.01 [1/min] as it has been documented in the parameters of the model's page. As a result, the phase bode plot of our initial system lead to a steep decrease in phase lag and thus in instability.

Figure 1 Bode plots for the initial transfer function (with d1=d2)

Furthermore, the wanted phase margin of 35degrees is accomplished with an operating period of 210 minutes meaning that our sensor should receive novel inputs every 210 minutes.

However, in the case where we would like to increase the operating frequency, we would have to increase at-least by 10-fold the one of the two poles. In that case, the instability of the system is dependent only by the new,larger pole. So, we propose the addition of a degradation tag which would increase the degradation rate of mneongreen on the receivers. A proposed degradation tag (ssRa) could be received from the paper [12].

As a result, the new bode plots are the following:

Figure 2. Bode plots for a 10-fold increase of d2, leading to d2=10*d1

, and we also received a phase margin of 35degrees for an operating period of 38minutes, a significant increase in allowed input frequency.

However, the above modelling effort comes with its assumptions. First and foremost the use of only two equations for the entire system would probably lead to quantitative inaccuracies about the exact frequency and phase values in each case. On the other hand, introducing a ‘pole splitting’ method doesn't come without its consequences, the main being an overall decrease in output value, which can be witnessed on the bode plots above. That is especially true in electronic components such as amplifiers which need to reach high gain values at the highest frequencies possible. In our case, it might not be of high significance, but it is possible that a 100-fold in degradation rates would lead to a stable system, where however the maximum fluorescence would be limited to < 1000 a.u.

So in order to effectively use such a modelling scheme to practice, a more analytical approach should be implemented. In our case, a lack of time meant that our team could never really use the biological perceptron technology as a real-time device, so such modelling wouldn’t be useful for further characterization of our current system. So the efforts of our dry lab modelling were mainly on the characterization of the sub-constructs and the effective prediction of the perceptron’s results for the optimization of the experimental situations (link for proof) However, the importance of introducing control theory principles and techniques of comprehending the frequency response of a design genetic system should not be understated. Based on our knowledge, only a few teams [13, 14] have used any type of control theory principles to their modelling efforts.

Pattern Recognition


The biological sensor our team aims to create, is intending to analyze complex stimuli. The system we propose was therefore originally designed to recognize 3x1 bits patterns. These patterns are created by on-off bits. We consider “on” a subpopulation of senders that is induced with aTc and “off” a senders subpopulation that is not induced with aTc, and as a result there will be no lactone production regardless of the RBS. The possible patterns are:

Figure 1: On-Off patterns. On is equivalent to an induced subpopulation of senders, Off equivalent to an uninduced subpopulation of senders.

We have created a library of 166 RBSs. The RBSs were divided into three families according to the rate of translation they show when followed by the LuxI gene. The RBSs that have high translation initiation rate are called “Strong (S)” RBSs, the ones with medium translation initiation rate are called “Medium (M)” RBSs and the ones with low translation initiation rate are called “Weak (W)”. So by creating different combinations of RBS triads we expect different amount of lactone produced.

Figure 2: Patterns by combining subpopulations of senders with RBS of different strength.


Thus combining RBSs with different strengths in populations that are induced or not, we create patterns that can be provided to the senders as inputs, which our system will be able to recognize. Theoretically, from each 3x1 bit pattern we expect a different amount of lactone production when all of the senders subpopulations are induced with aTc. Since the strength of each RBS is translated into an amount of lactone produced and finally the subpopulations of senders are mixed, the position of each bit does not matter but only what RBS each subpopulation has and if the population is induced or not. Furthermore, it is obvious that the maximum amount of the lactone will be produced in case no. 2 where the subpopulations of the senders all have very strong RBSs. Respectively, the lower amount of the lactone is expected in case no. 4 where the subpopulations of the senders all have very weak RBSs.



Our proposed implementation of the perceptron weights with the use of synthetic RBS sequences could become the template for the development of plug-and-play devices that enable the modification of the inputs' weights regarding the desired behaviour and application. We selected to use RBSs instead of promoters for the pattern recognition because of their versatility and it is quicker to and easier to predict the translation initiation rate of an RBS instead of directly evolutionize a promoter to have the desired strength. Overall, such a multicellular system has various applications in many scientific fields, from medicine to the environment and the thriving field of biocomputing. It could become the starting point for the development of “smart biosensors” or “smart drugs” and even perform complex computational tasks using much less energy than the amount needed by a computer, aiming to solve one of the most popular problems nowadays, the high cost of computational performances



Use of our model to simulate a pattern recognition task


After the characterization of our model, we wanted to simulate both the senders and the receivers populations and the way they would interact with each other. In this way, not only would we get a validity of our characterizations, but also make a direct comparison of the OpLo and PFR receivers. But at first we need to set up the pattern to be recognized as well as the RBS’s to be used.


For that purpose we decided to use only the sender’s constructs that were tested experimentally and accurately determine the a_luxi parameter of our model which depicts the strength of the AHL produced.

RBS Value Wet Value Dry a_luxi_fitted
BCD2 43704 43680 128
BCD12 16981 16209 47.5
S11 34943 34125 100
S13 19903 20475 60
S4 5926 5972 17.5
S5 597 597 1.75
S10 17106 17062 50
S12 2887 2900 8.5
BCD8 4050 4095 12
BCD14 23202 23205 68
Table 1. Model fit on the different sender’s RBS according to the experimental data

Also, the pattern simulated is presented below:

Index Senders 1 Senders 2 Senders 3 Desired Output
0 0 0 0 0
1 0 0 1 0
2 0 1 0 0
3 0 1 1 1
4 1 0 0 1
5 1 0 1 1
6 1 1 0 1
7 1 1 1 1
Table 2. Pattern to be simulated on our model


Via training from our software for the pattern above, it was concluded that the sequences to be used should be the following:



Figure 3. Received weight values and RBS constructs for each of the senders populations (X1,X2,X3).

,meaning that we should use 1 strong RBS and two medium to low ones.


For the sender's 1 population we decided to use our strongest RBS, BCD2 due to the fact that we want output = 1 whenever the sender’s population is induced. As for the other 2 populations, we decided to use a medium RBS for each one, BCD12, since we would like to have no fluorescence when only each one of these sender’s population is induced (see index 1,2). On the other hand, when both senders 1 and 2 are induced we expect that the strength of two medium RBS combined would lead to an OC6 production similar to that of BCD2 and thus lead to high fluorescence.


It should be mentioned that we didn’t use synthetic 6 RBS-which was the result of the training on the application-for the simulation due to the lack of experimental data for that construct. However, the BCD12 RBS is also a relatively suitable candidate as a medium to low RBS, always in comparison to the strongest one, which is BCD2 (as shown in Table 1).


Simulation Results


We performed a simulation of our system which can be found on the ‘proof_all.ipynb’ code on gitlab. In particular, we used a 1/7 senders/receivers cell culture ratio since we observed that it was the ideal ratio for our pattern recognition task. The culture volume remained 200μL. For more details on that you can take a look at our model analysis on the proof of concept page .


As we’ve seen on the characterization part of our model, the sender’s model takes as input the aTc concentration for induction and gives as output the molecules of external AHL. On the other hand, the receiver’s model uses concentrations of input OC6-instead of molecules-in order to estimate the output mneongreen. So, for that specific task, we added a simple expression to the implementation code that converts molecules of OC6 to Molarity:


Molarity_AHL=Molecules_AHL/(NA*V_bead), where:


NA =6.023*10**(23) is avogadro’s number

V_bead=200μL, the total cell volume


Other than that, we simulated the entire model system and received the following results:


Figure 4. Direct Comparison of the OpLo and PFR receivers constructs.

From the results above, the significance of the steep activation function is displayed; for the high value patterns the AHL produced from BCD2 is enough to give us the wanted output 1. However, the analog nature of OpLo means that even the medium RBS alone can offer a fluorescence output similar to that of fully induced BCD2. On the other hand, such a concern doesn’t exist on the PFR construct, where the receivers are not activated only with the induction of one medium RBS. Moreover, we witness that in the pattern with index 2, 2 medium RBS’s are more than capable to induce the receiver’s population.



Pattern Recognition Code

Diffusion

Overview


So far in the model we have analyzed the change of chemical species in the subpopulations of senders and in the population of receivers separately. In our project this year, however, we did not intend on a simple mixing of the two populations in a common container. Our purpose was that the cells of the senders and the receivers did not mix but only the chemical species, and mainly the lactone molecule, "travelled" from one population to another. For this reason, the team decided to use 6-well plates which had a geometry that allowed us to separate the two populations with a membrane.

More specifically, each well consists of two concentric disks, the inner and the outer. In the outer ring the populations of the senders were placed, while in the inner disk the population of the receivers. As mentioned above, the populations were separated by a membrane that has pores with a size of 0.4 nm, so it only allows chemicals to pass through it and not cells. The one with the suffix geometry was ordered by the company SPL Life Sciences [15].

So while we were able to model the change in the concentrations of the chemical species in each culture separately using differential equations, we had to think of an alternative way to approach the mechanism of diffusion across this membrane. We thought that this particular membrane can constitute and be expressed as "resistance" to the lactone molecule as it passes from the senders to the receivers. And so one way we could schematically show this was by using the COMSOL software, which is suitable for analyzes of mass transfer phenomena.

Creation of Our “Geometry” Using COMSOL


For this specific analysis we used the COMSOL software and more specifically version 5.2. Initially, it was chosen to perform a two-dimensional and time dependent analysis with “Transport of Diluted Species” as a selected Physics. We created the geometry of each well using concentric disks with the corresponding dimensions, one for the outer ring (where the senders populations would be), one for the membrane and one for the inner disk (where the receiver's population would be). We added “water” as the material of the outer ring and the inner disk and added a blank material for the membrane. For this material we tried to determine as many of the basic properties as we could since we did not have much data available for the membrane except those written on the technical data sheet of the product (density, permeability etc). So basically our basic geometry for the analysis was:

Figure 1: Basic geometry

The Analysis


At first we had to import the data concerning the production of the lactone by the senders in the outer ring. We solved our ODEs for the time frame of 1000 min and the data we collected were imported in COMSOL as an interpolation:

Figure 2: Data collected for ODEs for the production of the lactone by the senders.

Τhen we fitted these data using matlab and obtained a time polynomial equation for the change in the concentration of lactone produced by the senders. So in this way, we imported the software how to change the lactone in the outer ring. Our substance of focus was of course the lactone and that’s why she was the only diluted molecule we analyzed and is added as a “c” consentration. As for the inner ring, we again determined through the differential equations of our kinetic model and thus introduced the lactone consumption rate as a “Reaction” to the inner disc. Our “Mesh” had an appearance that looks like this:

Figure 3: The Mesh.

As is evident near the membrane the lactone will find greater resistance to pass than in the aqueous medium of the disks. The time dependent analysis took place for the first 1000 min with a step of 1 min.

Results


The results of this analysis were not as accurate as we expected, and the phenomenon was taking place much faster than we believed, but in general followed the pattern we had in mind. This may be due to the fact that we did not have the necessary data to fully characterize the membrane, or to the fact that no member of our team was fully familiar with the software. Precisely for this reason, we present some indicative qualitative results of our analysis with the aim of showing how we expect the change in concentration and movement of the lactone spatially in our system. The concentration of lactone in each disk at different time points of the analysis is shown schematically and annotated.

Figure 4: Indicative results.

As can be seen from the image, at 6 seconds the production of lactone has started from the centre (we consider moment 0 the moment when the lactone production begins) and has started to diffuse to the outer ring. Afterwards, the senders begin to produce a large amount of lactone, but it has difficulty getting through the membrane, while at the same time the amount that manages to actually go through the membrane is consumed by the cells of the receivers located in the inner disc. As the phenomenon progresses, more and more of the lactone produced by the sender penetrates the membrane for the first 450 minutes, where its production rate reaches its maximum. Afterwards, a decrease in the amount of lactone in the outside is observed. This probably happens because the rate of its production by the senders decreases, and it diffuses to the inside of the membrane as it continues to be consumed there. A reasonable question that arose at this point is why the lactone would still be able to diffuse the interior. After communicating with Dr. Kavousanakis, we were informed that this happens because it is consumed in the internal disc, and we can consider the membrane as a resistance whose value is proportional to the difference in the lactone concentration inside and outside the membrane.

In the future, we would like to import all our ODEs in COMSOL so that we don’t lose any useful information because of the fitting. We would also like to find a way to add a flux at the interface of the senders and the receivers to simulate the membrane as an actual “resistance”. We believe that in this way we could have much more reliable results, but due to the lack of time and knowledge we implemented a simpler system.
1. Team HSiTAIWAN 2016 2. Team Valencia UPV 2018 3. Team Zurich 2014 4. Harvard Bionumbers https://bionumbers.hms.harvard.edu/bionumber.aspx?id=107873&ver=6&trm=rna+polymerase+footprint&org= [Accessed: 24/7/2022] 5. Meyer, A., Shapiro, T., Glassey, E., Zhang ,J., Voigt, C. (2018). Escherichia coli “Marionette” strains with 12 highly optimized small-molecule sensors . Nature Chemical Biology volume 15, 196–204. https://doi.org/10.1038/nmeth.2404 6. Barati, Z. (2011). Parameter Estimation of Nonlinear Muskingum Models Using Nelder-Mead Simplex Algorithm. Journal of Hydrologic Engineering 16(11):946-954 . DOI: 10.1061/(ASCE)HE.1943-5584.0000379 7. A. W. K. Harris, C. L. Kelly, H. Steel and A. Papachristodoulou, "The autorepressor: A case study of the importance of model selection," 2017 IEEE 56th Annual Conference on Decision and Control (CDC), 2017, pp. 1622-1627, doi: 10.1109/CDC.2017.8263882. 8.Meyer AJ, Segall-Shapiro TH, Glassey E, Zhang J, Voigt CA. Escherichia coli "Marionette" strains with 12 highly optimized small-molecule sensors. Nat Chem Biol. 2019 Feb;15(2):196-204. doi: 10.1038/s41589-018-0168-3. Epub 2018 Nov 26. PMID: 30478458. 9. Process Control: Understanding Dynamic Behavior [Accessed 9/10/22] 10.Analysis and design of analog integrated circuits [Accessed 9/10/22] 11.Structural basis of ClpXP recognition and unfolding of ssrA-tagged substrates [Accessed 9/10/22] 12.Structural basis of ClpXP recognition and unfolding of ssrA-tagged substrates [Accessed 9/10/22] 13. Team Fudan-CHINA 2018 14. Team Hong_Kong_HKUS 2016 15. SPL life sciences, 6-well plates. http://www.spllifesciences.com/en/m21.php?cate=1&idx=307

© 2022 - Content on this site is licensed under a Creative Commons Attribution 4.0 International license