# Team:UPF Barcelona/Model

UPF_Barcelona

In this page you will find the development and results of the mathematical models used for our T3 biosensor and the lactone proof of concept.

Modelling

# Introduction: the importance of modelling

Mathematical modelling is a powerful tool that helps us to understand the behaviour and dynamics of living systems. This tool may explain the effects of different components of the system and, more importantly, dynamics that are difficult or impossible to test in the laboratory.

Moreover, it allows us to characterize and quantify the amount of protein that is being expressed, which is something really important for any kind of project, but specially for Hormonic, as we are dealing with feedback restoration.

In this page you can find two mathematical models developed in order to characterize and understand the behaviour of a intein-mediated T3 biosensor and a lactone feedback system. Furthermore, these models were fitted into experimental data and sensitivity and stability analysis were performed.

On top of that, a model to characterize the ASV tag was designed and then, we mathematically demonstrated its ability to increase the protein degradation rate, thus reducing sensing time.

# Model development

## Introduction

The hormonal triiodothyronine (T3) biosensor is a key aspect for Hormonic, and one of the project bases. This biosensor works at the post-translational level and it is mediated by inteins, proteins that carry out protein splicing [1]. When the hormone binds to the receptor it splices, reconstructing the split fluorescent protein that is at both intein extremes (this is widely explained and detailed at the T3 Biosensor section). Therefore, the fluorescent protein is only active when the splicing of the intein occurs, that is, if the hormone is present in the medium. Thanks to this biological agent capable of translating biological information into measurable values we are able to sense the amount of triiodothyronine (T3) in the medium. Consequently, a model was developed to understand the biosensor behaviour and to characterize it.

## T3 Biosensor model

A schematic representation of the interactions found on the T3 sensor cell is shown below (Fig.1). As it can be seen, when T3 binds to the receptor, the intein splices, reconstructing the green fluorescent protein (GFP), thus generating fluorescence that can be measured. We must consider that this reconstruction is not 100% reliable, thus leading also to non-functional fluorescent proteins [2]. Moreover, the construct is preceded by a constitutive promoter, therefore, the intein is produced constitutively in the cell, although ideally there is only fluorescence when T3 is present in the media and the intein domain is activated.

From these processes we can develop 4 reactions (Fig.2) which are the basis for our T3 sensor cell model. In these reactions P refers to the intein mediated biosensor protein as a whole (including the receptor, the intein and the two halves of the green fluorescent protein).

Due to the experimental results, we must consider that there is a basality in the production of fluorescent protein. Moreover, we are taking into account that the variable P always refers to the whole intein construct, including the split fluorescent protein at the end sites. Therefore, when fluorescent protein is generated, this construct is degraded. This makes sense as the residual construct is non-functional, since it has already spliced the green fluorescent protein (GFP). Furthermore, as said before, it is necessary to multiply the production of GFP for a rate, as the reconstruction can lead to non-functional GFP. The final ODE system and its parameters and variables can be seen below (Eq.2, Tab.1-2). The parameters were fitted with the fluorescence normalized from 0 to 1 and non-normalized.

 Parameters Description Units Fitted values Normalizedfitted values αP Production rate of P (a.u.)/time 7.14·103 5.89·103 δP Degradation rate of P 1/time 4.85·10-2 5.03·103 β Production rate of GFP 1/(time·(a.u.)) 8.74·103 9.07·103 αGFP basal production rate of GFP RFU/time 4.00·103 2.63·103 δGFP Degradation rate of GFP 1/time 5.14 6.28·103 β' Production rate of functional GFP 1/(time·(a.u.)) 7.14·103 5.89·103
 Variables Description Units P Intein mediated biosensor protein concentration Arbitrary units (a.u.) T3 Triiodothyronine (T3) concentration Molar concentration (M) GFP Green fluorescent protein (GFP) concentration Relative fluorescence units (RFU)

From the ODE system (Eq.1) a transfer function can be derived if P and GFP are on a steady state (dy/dt=0). This transfer function (Eq.2) allows us to connect the concentration of T3 in the media to the fluorescence that the sensor cells are emitting. It results on a Hill function with a low Hill coefficient (n=1), thus meaning that the sensor will not be hypersensitive and it will have a wide dynamic range where it can sense the T3 concentration.

## Biosensor sensitivity

As we are developing a sensor, it is key to find properties that provide us detail on its ability to detect and react to changes in the measured physical quantity. In order to do so, we decided to calculate its sensitivity.

The sensitivity is defined as the minimum value of an input parameter that is capable of creating a detectable output change. In our case, it will be the minimum change in T3 concentration (M) that generates a change in the GFP concetration, expressed in relative fluorescence units (RFU). In order to calculate the sensitivity, we used its mathematical definition (Eq.3.1), which says that the sensitivity is the derivative normalized to the input and output. The resulting sensitivity value can be found below in the modelling results section.

## Stability analysis

In order to know if the steady state calculated on the transfer function (Eq.2) is stable, the Jacobian of the ODE system needs to be computed. The stability of the system is an important metric to know if our system will converge to this steady state (stable) or it will diverge from it (unstable). To obtain the Jacobian (J) (Eq.4.1) we must perform the derivative of the ODE of the system (f()=dP/dt, g()=dGFP/dt). From the Jacobian we can compute the trace [tr(J)] and the determinant [det(J)] so that we can analyze the stability of the T3 sensor cell system.

Since the determinant is positive and the trace is negative (Eq.4) we can affirm that the system has a stable steady state [3]. This means that, even if the system is far from the steady state, for example when we include T3 to the media and GFP has not yet been formed, the system will evolve towards the steady state, thus sensing the T3 concentration.

## Introduction

As a proof of concept for our device, a biological circuit capable of producing and sensing lactone in a controlled way was designed. With this circuit, we can simulate and understand how the thyroid axis works and, consequently, how our device would work. This circuit is formed by two types of cells. On the one hand we have the producer cell, which produces lactone in presence of arabinose and in absence of glucose in the media. On the other hand we have the sensor cell, which will produce super-folder GFP (sfGFP) in presence of lactone in the media.

## Lactone sensor cell model

A schematic representation of the interactions that happen in the sensor cell is shown below (Fig.3). As it can be seen, the sensor cell constitutively produces LuxR. When formed, this protein can bind to lactone (AHL), thus generating a complex that activates the pLux promoter, generating the sfGFP whose fluorescence can be detected.

From these interactions we can develop 5 reactions (Fig.4) which are the basis for our lactone sensor cell model.

If the binding and the unbinding of LuxR and AHL is in equilibrium, and the production of sfGFP has some basality (full demonstration on the Modelling notebook), then the model can be simplified into 2 ordinary differential equations (ODE) by applying the law of mass action (LMA) (Eq.5). The LMA states that the rate of the chemical reactions (dLuxR/dt) is directly proportional to the concentrations of the reactants (LuxR). The final ODE system and its parameters and variables can be seen below (Eq.5, Tab.3-4). The parameters were fitted with the fluorescence normalized from 0 to 1 and non-normalized.

 Parameters Description Units Fitted values Normalizedfitted values αLuxR Production rate of LuxR (a.u.)/time 9.99·103 9.99·103 δLuxR Degradation rate of LuxR 1/time 1.68 1.68 BLuxR Type of regulation (If BLuxR>1 LuxR acts as an activator) no units 6.62 6.62 wLuxR Dissociation constant 1/(M·(a.u.)) 9.99·103 9.99·103 αGFP Production rate of GFP RFU/time 9.12·103 8.63·103 δGFP Degradation rate of GFP 1/time 17.68 1.36·103
 Variables Description Units LuxR LuxR concentration Arbitrary units (a.u.) AHL Lactone (AHL) concentration Molar concentration (M) GFP Green fluorescent protein (GFP) concentration Relative fluorescence units (RFU)

From the ODE system (Eq.5) a transfer function can be derived if LuxR and GFP are on a steady state (dy/dx=0). This transfer function (Eq.6) allows us to relate the concentration of lactone (AHL) in the media with the fluorescence that the sensor cells are emitting. The resulting graph has the shape of a Hill function with a low hill coefficient (n=1), thus meaning that the sensor will not be hypersensitive and it will have a long dynamic range where it can sense the AHL concentration [6][7].

## Biosensor sensitivity

As we are developing a sensor, it is key to find properties that provide us detail on its ability to detect and react to changes in the physical quantity measured. In order to do so, we decided to calculate its sensitivity.

The sensitivity is defined as the minimum value of an input parameter that is capable of creating a detectable output change. In our case it will be the minimum change in AHL concentration (M) that generates a change in the fluorescence (RFP). In order to calculate the sensitivity (Eq.7) we used its mathematical definition, which says that the sensitivity is the derivative normalized to the input and output. The resulting sensitivity value can be found below in the modelling results section.

## Stability analysis

In order to know if the steady state calculated on the transfer function is stable, the Jacobian of the ODE system needs to be computed. The stability of the system is an important metric, it is used to know if our system will converge to this steady state (stable) or if it will diverge from it (unstable). To obtain the Jacobian (J) we must perform the derivative of the ODE of the system (f()=dLuxR/dt, g()=dGFP/dt). From this matrix we can compute the trace [tr(J)] and the determinant [det(J)] so that we can analyze the stability of the AHL sensor cell system.

Since for positive δLuxR and δGFP the determinant is positive and the trace is negative (Eq.8) we can affirm that the system has a stable steady state [3]. This means that even if the system is far from the steady state, for example when we include AHL to the media and GFP has not yet been formed, the system will evolve towards the steady state, thus sensing the AHL concentration.

## Introduction

As a proof of concept for our device, a biological circuit capable of producing and sensing lactone in a controlled way was designed. With this circuit, we can fully simulate and understand how the thyroid axis works and, consequently, how our device would work. This circuit is formed by two types of cells. On the one hand we have the producer cell, which produces lactone in presence of arabinose and in absence of glucose in the media. On the other hand we have the sensor cell, which will produce super-folder GFP (sfGFP) in presence of lactone in the media.

## Lactone producer cell model

A schematic representation of the interactions that take place in the producer cell is shown below (Fig.5). As it can be seen, this cell produces LuxI in response to high concentration of arabinose (ARAB) and low concentration of glucose (GLU). When LuxI is formed, this enzyme produces lactone (AHL).

Due to the fact that the functioning of the pBad promoter, which controls the expression of LuxI, is complex and has several intermediate agents, a detailed description of these intermediate reactions falls out of the scope of this project. A model’s complexity should always be dependent on the results that need to be obtained, and for us, a detailed characterization of the pBAD promoter is not necessary. So, if we take into account the overall mentioned process we can find that the system has 5 reactions, which can be found below (Fig.6).

If the reactions involving the Pbad promoter are in equilibrium (full demonstration on the Modelling notebook), then the model can be simplified into 2 ordinary differential equations (Eq.9) by applying the law of mass action (LMA). The LMA states that the rate of the chemical reactions (dLuxI/dt) is directly proportional to the concentrations of the reactants (glucose or arabinose). The final ODE system and its parameters and variables can be seen below (Eq.9, Tab.5-6).

 Parameters Description Units Fitted values αLuxI Production rate of LuxI (a.u.)/time 26.72 δLuxI Degradation rate of LuxI 1/time 1.15 wLuxI Dissociation constant 1/M 1.28·104 αAHL Production rate of AHL 1/time 26.72 δAHL Degradation rate of AHL 1/time 1.15
 Variables Description Units LuxI LuxI concentration Arbitrary units (a.u.) AHL Lactone (AHL) concentration Arbitrary units (a.u.) ARAB Arabinose (ARAB) concentration Molar concentration (M) GLU Glucose (GLU) concentration Molar concentration (M)

From the ODE system (Eq.9) we can derive a transfer function if we take into account that LuxI and AHL are on a steady state (dy/dx=0). This transfer function (Eq.10) allows us to relate the concentration of arabinose (ARAB) and glucose (GLU) with the concentration of lactone (AHL) that the cells will produce.

## Stability analysis

In order to know if the steady state calculated on the transfer function is stable, the Jacobian of the ODE system needs to be computed. The stability of the system is an important metric, it is used to know if our system will converge to this steady state (stable) or if it will diverge from it (unstable). To obtain the Jacobian (J) we must perform the derivative of the ODE of the system (f()=dLuxI/dt, g()=dAHL/dt). From this matrix we can compute the trace [tr(J)] and the determinant [det(J)] so that we can analyze the stability of the AHL producer cell system.

Since for positive δLuxI and δAHL the determinant is positive and the trace is negative (Eq.11) we can affirm that the system has a stable steady state [3]. This means that even if the system is far from the steady state, for example when we include arabinose to the media and AHL has not yet been formed or when we add glucose and the AHL levels have not yet been reduced, the system will always evolve towards the steady states found on the transfer function (Eq.10).

## Introduction

Since there is a need for faster response in the system we aim to build, a faster recovery to the initial state of the system is also needed. Faster responses are a key aspect when it comes to our system, as the feedback of the PID controller needs discrete data in a short period of time to be able to regulate. To achieve this we must increase the protein degradation rate but we also have to consider that a strong increase would imply a loss of accuracy on the measurements, as the range will become narrower. Therefore, the addition of the ASV tag, which is supposed to be a weak one, becomes ideal [8]. If you want a more complete information on why we introduced the ASV tag please check the Implementation page.

Two models have been developed to characterize the implications of ASV tag in the degradation rate and to demonstrate that a tag is needed for a faster recovery to the initial state.

## ASV tag model

A schematic representation of the interactions on the two systems is shown below (Fig.7). Both systems are preceded by the PJ23100 constitutive promoter, meaning that the production rate of sfGFP will be the same for both systems. However, the introduction of a tag will have an impact on the sfGFP degradation rate, thus changing this parameter in the second system. All in all, both systems consist of a simple constitutive expression plus a degradation rate.

From these interactions we can develop both system reactions (Fig.8) which are the basis for our ASV tag characterization models.

From the above system reactions (Fig 8.) we are able to develop an ordinary differential equation (ODE) for each system. Both systems have a constitutive production rate and a degradation rate that depends on the amount of protein in the medium. The ASV tag will imply an increase in the degradation rate, therefore, the rate will be different for both systems. The final ODE system and its parameters and variables can be seen below (Eq.12, Tab.7-8).

 Parameters Description Units αP Production rate of sfGFP RFU/time δP Degradation rate of sfGFP 1/time δPASV Degradation rate of sfGFP with the impact of the ASV tag 1/time
 Variables Description Units sfGFP Super-folder green fluorescent protein (sfGFP) concentration Relative fluorescence units (RFU)

Once we know the ODEs from both systems, we integrate them to know how the concentration of both proteins will change over time (Eq.13). This function will allow us to calculate the steady-states and to demonstrate that an increase of degradation implies a reduction of time. Full integration and development is available in the Modelling notebook below.

## Steady-states & stability analysis

In order to obtain the steady-state values of both systems, we perform the limit , when t tends to infinity, of the functions obtained before. This is due to the fact that the steady-state values are the ones that the system tends to end. The steady-states of both systems can be seen below (Eq.14).

Once the steady-states are known, we must compute the stability. The stability of the system is an important metric to know if our system will converge to the steady state (stable) or it will diverge from it (unstable).

As the derivative of both ODEs , with respect to the protein, is negative we can conclude that the systems will converge to a stable equilibrium point [3].

## Relationship between the degradation rates and its implications

As both systems are constitutively regulated by the same promoter, the PJ23100 (BBa_J23100), the production rate will be the same in the two systems. Due to this fact, we are able to compare the two degradation rates, thanks to the steady-states relations computed above (We take the last point of the experimental data results as the steady-states).

In the relation between the degradation rates (Eq.15) we can clearly see that the ASV tag has an impact on the degradation rate. Mathematically it has been proved that the degradation rate is increased , approximately, 25% when the ASV tag is added on the sfGFP gene sequence. Then, a faster response will be achieved without the implication of a loss of accuracy on the measurements, as seen in the Proof of Concept.

## Demonstration that an increase of degradation rate implies a reduction of time

In this section, we will mathematically demonstrate that an increase of degradation rate, as the one performed by the ASV tag, will imply a reduction of the system recovery time. This reduction is needed for having faster responses in our feedback system as we will have more data in the same period of time, meaning that we will be able to better regulate the system.

To demonstrate this we calculated the protein half-life, the time when the protein decays until half of its initial concentration. In our context it will be when the steady-state is reached and there is no more input, for example in our sensor cell. For this, we will model a cell that starts at the steady-state concentration and only degrades. Due to this model, we will be able to demonstrate that the tag is needed to reduce the recovery time. All the calculations and demonstrations to obtain the following expression are in the Modelling notebook.

In conclusion, as seen above (Eq.16), if the degradation rate of a protein is increased it will always imply a decrease in the protein half-life.

# Model results and discussion

## Characterization of the intein-mediated T3 biosensor

The system transfer function (Eq.16) was fitted into the experimental data in order to see if our model could follow it. This data was extracted from a Plate-Reader analysis on GFP fluorescence emission with concentrations of T3 going from 1 pM to 100 µM. From these results, we can conclude that there is a high correlation between the experimental and the modelling results, thus indicating that the developed T3 biosensor model reflects the real phenomena that happen in the T3 sensor cell. Consequently, the fitted transfer function could be used to predict specific T3 concentrations measuring the fluorescence (the inverse of the fitted transfer), which will be useful for our close-loop system. All the calculations of the functions seen in the model development section can be seen at the modelling notebook.

From the graph above (Fig.8) we can obtain the biosensor dynamic range, that goes from 10^-7 to 10^-3 M, which is the range of T3 concentrations at which the biosensor can produce a change in GFP expression, that is, the sensor operative range.

Moving on, using the sensitivity formula calculated before (Eq.3), and the fitted values, we can obtain the biosensor sensitivity. A sensitivity around 0 and 1 (ξ=0.2107) means that the sensor is neither robust (which would mean that the biosensor cannot react to different T3 concentrations), nor hypersensitive (which would mean that the biosensor overreacts to minimal changes in the T3 concentration, thus shrinking the dynamic range), respectively.

## Phase planes and portraits

To get important information about how the biosensor performs and works, phase planes and phase portraits were computed (Fig 9.). Phase planes indicate the nullclines and the vector field, with phase portraits adding some trajectories with different initial conditions into the plot.

On the one hand, the nullclines are the curves where an ODE is zero. Therefore, since our system has two ODE (LuxR and GFP), we will have two nullclines in all of our phase planes and portraits. As nullclines are where an ODE is zero, the point where the two nullclines cross will correspond to the steady state of the system (all ODE are zero).

On the other hand, the vector field indicates the evolving direction of the system for all kinds of initial conditions. This allows us to see the biosensor behaviour for different initial conditions without having to do experiments in the lab.

As we can see above (Fig.9), if the concentration of T3 increases we find the steady state at a higher GFP value, this is due to the fact that GFP nullcline will have a steeper slope. Moreover, it can be seen that the steady states, found in the phase planes and portraits, are the same ones from the transfer function, which is coherent. In the phase planes, we can also see that the vector field indicates that the steady state is stable, since all the vector field arrows finally converge at it. This result coincides with the stability analysis done in the T3 biosensor modelling section [4][5]. Furthermore, we can see that in the phase portraits, although the trajectories have different initial conditions of GFP and P (Intein mediated biosensor protein), they all converge at the steady state.

## Discussion

Being able to understand the behaviour of the intein-mediated T3 biosensor, and to characterize it through mathematical modelling, is important to ensure that the dynamics of the system are known. If the dynamics of the T3 biosensor are understood, we could be able to predict the concentration of T3 when fluorescence is detected, which is the first key aspect to start restoring the thyroid feedback loop. This is due to the fact that, if the hormone concentration is known, we could be aware of the thyroid state, and more importantly, , if there is a lack or overexpression of T3. The stability analysis performed could be extrapolated to any type of system, providing a powerful tool to better understand the system dynamics. Moreover, this T3 biosensor model sets a precedent for any kind of future intein-mediated biosensor.

In conclusion, this model represents the first feasible step towards the restoration of the thyroid feedback loop. As the dynamics of the T3 biosensor have been understood, we are now able to predict the concentration of T3 when fluorescence is detected.

## Sensor cell results

The transfer function (Eq.22) was fitted with experimental data to see if our model could follow real data. This data was extracted from a Plate-Reader analysis on GFP fluorescence emission with concentrations of lactone going from 100 pM to 100 µM. From these results we can conclude that there is a high correlation between the experimental and the modelling results, thus indicating that the sensor cell model we have developed reflects the real phenomena that happen in the cell.

Using the formula calculated on the lactone sensor modelling section we can calculate the sensitivity of our sensor. A sensitivity around 0 and 1 (ξ=0.3688) means that the sensor is neither robust (which would mean that the sensor can’t react to different AHL concentrations), nor hypersensitive (which would mean that the sensor overreacts to minimal changes in the AHL concentration, thus shrinking the dynamic range). Another important parameter that can be extracted from the graph above (Fig.10) is the dynamic range, which is the range of AHL concentrations at which the sensor can produce a change in GFP expression. This was found to be from 1nM to 1uM, giving 3 orders of magnitude where AHL can be sensed.

To get important information about how the biosensor performs and works, phase planes and phase portraits were computed (Fig.11). Phase planes indicate the nullclines and the vector field, with phase portraits adding some trajectories with different initial conditions into the plot. The nullclines are the curves where an ordinary differential equation (ODE) is zero. Therefore, since our system has 2 ODE (dLuxr/dt and dGFP/dt), we will have 2 nullclines in all of our phase planes and portraits. As nullclines are where an ODE is zero, the point where the 2 nullclines cross will correspond to the steady state of the system (all ODE are 0). Moreover, the vector field indicates the evolving direction of the system for all kinds of initial conditions [4][5]. This allows us to see what is the behaviour of the sensor cell at different initial conditions without having to do experiments in the lab.

As we can see above (Fig.11), if the concentration of AHL increases we find the steady state (intersection between the 2 nullclines) at a higher GFP value, this is because the GFP nullcline (dGFP/dt=0) goes from an horizontal line to a Hill function. In the phase planes we can also see that the vector field indicates that the steady state is stable, since all the vector field arrows finally converge at it, which coincides with the stability analysis done in the proof of concept modelling section. Furthermore, we can see that all trajectories, which correspond to different GFP and LuxR initial conditions, will converge to the steady state of the system.

## Minimal hormonal regulation model results

For the whole genetic circuit, we characterized its behaviour under different concentrations of both, arabinose and glucose. This data was extracted from a Plate-Reader analysis on GFP fluorescence emission with concentrations of glucose ranging from 100mM to 10uM and concentrations of arabinose ranging from 1mM to 10nM. This data was fitted using the transfer functions of the producer and sensor cells (Eq.23).

As we can see above (Fig.12), lactone will only be produced with low glucose concentrations and in high arabinose concentrations. Furthermore, we can see that our model follows the experimental data, thus showing that it reflects the real phenomena happening in the producer cells even though some assumptions were made.

## Discussion

Understanding how the concentration of a molecule can be controlled by an external device is key for the treatment of several diseases such as thyroid disease and diabetes. With the proposed models of the proof of concept we have demonstrated that we can produce and sense lactone in a controlled way via the action of effector molecules (glucose and arabinose). These results show us that a lactone closed circuit that incorporates a PID controller is a viable option for the demonstration of Hormonic’s functionality. This means that the next step of this project would be to combine the sensor and producer cells of the proof of concept with the turbidostat circuit.

# Experiments and simulations

For the characterization of all the cell types described in this page Plate-Reader analysis were made. Although GFP emission was measured to check the sensor activity in all the experiments, RFP emission was only used when the experiment required 2 cells per well. Measuring the RFP emission allowed us to know the optical density (OD) of the 2 cells separately. All the information on the experimental conditions and parameters used are described on the table below (Tab.9).

 Experimental parameters and conditions Description Plate-Reader model Synergy HTX Plate type Thermo Fisher 96-well microplates black-walled clear bottom Cell medium LB Time 24 hours Shake Liner: Continuous Frequency: 567 rpm (3 mm) Temperature 37 ºC Gain 50 Optical Density (OD) measurement (absorbance) 660nm GFP excitation wavelength 485nm GFP emission wavelength 528nm RFP excitation wavelength 575nm RFP emission wavelength 625nm

For the data fitting algorithms, used in the characterization of the T3 biosensor and lactone proof of concept, and the generation of graphs, Python 3 was the selected programming language, using Jupyter 6.0.3 notebooks as our interface. In addition, a notebook by Doulcier G. was adapted to produce the phase planes and portraits seen in the results section [9]. All the code used and the modelling notebook can be downloaded from the table below (Tab.10).

 Files T3 biosensor characterization and phase planes Minimal hormonal regulation model characterization and phase planes Modelling notebook

## References

[1] Skretas G, Wood DW. Regulation of protein activity with small-molecule-controlled inteins. Protein Sci. 2005;14(2):523-532. doi:10.1110/ps.04996905

[2] Gangopadhyay JP, Jiang SQ, Paulus H. An in vitro screening system for protein splicing inhibitors based on green fluorescent protein as an indicator. Anal Chem. 2003;75(10):2456-2462. doi:10.1021/ac020756b

[3] Layek, G.C. “Stability Theory” in An Introduction to Dynamical Systems and Chaos 1rst ed. Springer, 2015, pp. 129-158.

[4] Strogatz, S.H. “Phase Plane” in Nonlinear Dynamics and Chaos, 1rst ed. Perseus Books Group, 1994, pp. 145–181.

[5] Layek, G.C. “Phase Plane Analysis” in An Introduction to Dynamical Systems and Chaos, 1rst ed. Springer, 2015, pp. 83-127.

[6] Ballestero, M. C., Duran-Nebreda, S., Monta, R., Solé, R., Macía, J., Rodríguez-Caso, C. A bottom-up characterization of transfer functions for synthetic biology designs: lessons from enzymology. Nucleic Acids Research, 2014, Vol. 42, No. 22.

[7] Garcia-Ojalvo, J., Elowitz, M.B., Strogatz, S.H. Modeling a synthetic multicellular clock: repressilators coupled by quorum sensing. Proc. Natl Acad. Sci. U.S.A., 101, 10955–10960.

[8] Andersen, J.,Sternberg, C., Poulsen, L. K., Bjørn, S. P., Givskov, M. New Unstable Variants of Green Fluorescent Protein for Studies of Transient Gene Expression in Bacteria. Applied and environmental microbiology, June 1998, p. 2240–2246.

[9] Fitzhugh-Nagumo model: an excitable system: Lecture Notes, retrieved August 31th, 2020, from https://www.normalesup.org/~doulcier/teaching/modeling/excitable_systems.html