Our Model

What is being modelled and why?

The goal of the model is to inform the wetlab with experimental design and simulate expected behavior for the experiments before or in parallel with the experiments in the wet lab. For this project, our proof-of-concept is a dose-dependent detection system using naringenin responsive promotor and mKATE expression as our read-out. Normally, naringenin will control the expression of NarX, hence affecting the dimerization of NarX, which playes an important role in expression of our fluorescent protein mKate. However, NarX mutant can bind with NarX, hence decreasing the amount of NarX dimerization in the cell. Therefore, the job of our model team is to find a threshold of naringenin concentration at which sufficient concentration of wildtype NarX can be produced to surpass NarX mutant (see SHOOT part for more detail) so that the therapeutic protein can be expressed.

In the first part of our Model, we gave an overview of the dynamic system in our experiment, which served as a base for our first mathematical modelling (the second part). Then, we moved to the third part, a novel method for modeling the system that we learned from our supervisor, Vitor. Finally, we conclude with the usability of our models.

Overview of the dynamic system in our experiment (SHOOT part)

The process started with Naringenin introduced into the system. The concentration of NarX is governed by P_Naringenin which responds to varying dosages of Naringenin. Upon increasing concentrations of NarX, they start dimerizing in the presence of nitrate. The dimerized NarX phosphorylates NarL inducing conformal changes in the molecule. These changes cause the phosphorylated NarL to bind and activate the NarL promoter P_NarL. The P_NarL is the final key component that promotes the release of mKate (Figure 1).

These chains of events are the lynchpin for our dose-dependent kinetic model. However, with the presence of NarX mutant as the latent variables, the dose-dependent model will change. NarX mutant, however, is not affected by P_Naringenin or naringenin but P_tet. The NarX mutant does not trigger any chain of subsequent reactions above but it competes with NarX for dimerization which undermines the dimerization of NarX with itself (Figure 2). Therefore, we require higher concentrations of NarX to overcome the NarX-NarX-mutant competitive dimerization. The specific concentrations of NarX required to commence the subsequent reactions reaching mKate are dictated by the concentration of NarX-mutant present in our system.

The overall system is presented in the figure below.

Base model

Kinetic modeling is a functional tool that provides essential information on the time evolution of the reaction mechanisms. With it, is possible to describe the structured interplay of protein interactions that constitute the dynamics of the biosensor circuit. Initially, the chemical reaction chain was formulated (Figure 4). The reactions are then used to propose the block diagram in MATLAB SimBiology editor (figure 5.).

With the SimBiology app, is possible to automatically pose and solve the system of equations that govern the system. However, the novel nature of the project leads to serious limitations on the parameters we can provide to the app. Many of the reactions have not been numerically characterized, and the experimental data required to do so would require multiple different projects spanning several years. For this reason, a different approach was used to model the system of equations posed with the SimBiology app.

Model Division

One of the main principles of engineering and systems theory states that a complex system can be tackled by dividing it in smaller, more manageable subsystems. Following this line of thought, four subsystems were defined for this system (Figure 6).

Step 1. NarX and Narx mutant expression.

As mentioned previously, the expression of NarX mutant is only dependent on P_tet and remains constant during each test. In contrast, the concentration of NarX depends on the concentration of naringenin, which varies over time. For this reason, only the function describing the NarX concentration is considered of interest.

with:

• a₁ = the basal normalized fluorescent signal (leaky expression, au);
• b₁ = the maximum normalized fluorescent signal (au);
• n = the Hill coefficient (cooperativity, sigmoid character);
• M₁ = the Hill constant (half-maximal naringenin concentration, mg/L)

The relationship between NarX and naringenin has been noted to follow a Hill equation (De Paepe et al., 2018). The Hill equation describes the binding of macromolecules (usually proteins) to ligands as a function of the ligand concentration. For this subsystem, Naringenin behaves like a ligand and NarX like the receptor protein. Using experimental data provided by the Center of Synthetic Biology of the University of Ghent , a Hill function was conducted to estimate the values of a₁, b₁ and M₁.

Step 2. NarX and NarX mutant dimerization

In this subsystem, three reaction are occurring simultaneously (equations (2), (3), (4)). The binding of NarX with other NarX, the binding of NarX with NarX mutant and the binding of NarX mutant with other NarX mutant.

To model these reactions a Mass Action Kinetic (MAK) model was implemented. In MAK, the rate of a chemical reaction is proportional to the concentration of the reactants. However, these reactions are reversible which requires association and dissociation rates (equations (5),( 6), (7)). For these equations, the association and dissociation rates were assumed to be equal.

Due to the speed of the reaction, the three types of dimerization can be assumed to be instantaneous. Similarly, the model assumes perfect dimerization (each NarX and NarX mutant binds to other NarX or NarX mutant). These assumptions require equilibrium to be reached instantaneously, leading to the following equality:

Hence,

With K_d as the ratio of the dissociation rate constant and the association rate constant (or the dissociation constant). Since the variable of interest for this subsystem is the NarX bound to other NarX (the only combination that would produce a fluorescent response) the system can be described in terms of its concentration. Let Ɵ be the ratio of NarX bound to NarX vs the concentration of all three possible bindings.

From Equation 8:

Then,

Additionally, since most of the reactions of the step 1 are quasi-instantaneous, by the time step 2 occurs the concentrations of NarX and NarX mutant are expected to be constant, resulting in:

With NarX ₀ the concentration of NarX after step 1.

Step 3. Phosphorylation of dimerized NarX and transphosphorylation of NarL

The phosphorylation NarX occurs as a response to the presence of nitrates in the medium and can be described by the following reversible reaction:

With ka and kd corresponding to the association and dissociation rates respectively.
Then the transphosphorylation process (the passing of the phosphoryl group from one molecule to other) occur at a rate of k₁. Simultaneously, the phosphorylated NarL degrades at a rate of d₁. Both processes are described by equations 11 and 12 accordingly.

As with previous steps, the reactions of the previous subsystem are considered instantaneous, making all the dimerized molecules and ligand concentrations constants. Since the reaction corresponds to the interaction of a macromolecule with its ligand, the subsystem can be modeled with a Hill equation. For the equation 10:

with K_d = k_d/k_a

Then the concentration of NarLNarX_dim is:

Hence,

As it can be observed, the resulting equation 13 is equivalent to the Michaelis-Menten Equation. To account for the degradation of NarLP, an additional term needs to be included, resulting in the equation:

Step 4. Phosphorylated NarL activates the promoter PNarL, prompting the expression of mKATE

For this step, the ratio between the concentration of phosphorylated NarL and the promoter PNarL is assumed to be 1:1. As a consequence of this assumption, the concentration of phosphorylated NarL can be directly linked to mKATE without the intermediate reaction with PNarL. The relation between the concentration of phosphorylated NarL and the expression of mKATE is described by k2 and the degradation of mKATE is given by d2.

As with other subsystems, the reactions previous to the start of this step are considered to have reached equilibrium. Thus, the equation 12 becomes:

Including the degradation of mKATE results in the following modified equations:

After reaching saturation, the expression of mKATE will be maximal. Since the concentration of mKATE will not change at this point, d[mKate]/dt = 0. The concentration of mKATE then becomes:

Integrated model

According to the postulates of system theory: “the whole is greater than the sum of its parts”. In other words, a system cannot be fully described only by its components. Following this line of reasoning, the integration of the parts is required to achieve a complete model. MATLAB Simulink was selected for this task due to its versatility and the possibility to observe the behavior of the model at different points of the chain. The integrated model is presented in the figure 7. Each of the scopes at the top mark the output of each of the subsystems (1, 2 and 3-4 respectively). For the constants, the same values as the ones used in each of the smaller models were used.

Results

Step 1

With the Hill function, the values of the unknown parameters were found using Scipy (Python): a1 = 153 (au), b1 = 17,088 (au), n = 2 (au) and M1 = 30 (mg/L). Moreso, the curve presented a goodness of fit R2 of 99.35%, indicating that 99.35% of the variance in the data can be explained by the model.

To establish the fluorescence response, the assumption that each unit of phosphorylated NarX will be equivalent to one unit of fluorescence was made. The range of NarX values was chosen based on the estimate of 1,000-100,000 copies of NarX in one cell (Milo, 2013) and further supported by the results obtained with the integrated model. As presented in the graph (Figure 8), the fluorescence intensity is expected to be saturated at about 16,000 arbitrary units in 1 cell and at a naringenin concentration of over 80mg/L.

Step 2

For this subsystem, the effect of the concentration of NarX mutant over the dimerization of NarX was of great interest. To test it, the model of the subsystem was evaluated in the range 0 to 20000 a.u for both the concentration of NarX and NarX mutant. The initial concentration of NarX (NarX₀) was assumed to be 50,000 (the median point of the expected range of NarX copies in one cell). Additionally, Kd was assumed to be 2 (iGEM 2020 Measurement Webinars). A plot of the function is presented in figure 9.

As observed on the figure 9, for every concentration of NarX, increasing the concentration of NarX mutant will decrease the concentration of dimerized NarX. Conversely, for each concentration of NarX mutant, the amount of NarX dimerized will depend on the concentration of NarX. However, the relationship will correspond to a sigmoid function that reaches saturation at 25,000 units.

Figure 10 provides more insight about the effect of the concentration of NarX and NarX mutant in the concentration dimerized NarX. For each level of dimerized NarX, the optimal concentration of NarX can be established. This concentration represents the minimum concentration required to have a given dimerized NarX concentration, disregarding concentration of NarX mutant. For example, if a concentration of dimerized NarX of at least 15,000 is required, a NarX concentration of 15,000 will suffice. Even if the concentration of NarX mutant were to be over 20,000, a concentration of dimerized NarX of at least 15,000 would be obtained.

Step 3/4

To observe the effect of the concentration of dimerized NarX in the expression of mKate, the function is evaluated in the dimerized NarX concentration range 0 to 1,000. The values for the constants correspond to k2*k1/d1 = 720 mole/min, d2 = 0.02 min^(-1), Kd = 2 mole (all are assumptions, based on iGEM 2020 Measurement Webinars). Additionally, the concentration of dimerized NarX was evaluated from 0 to 10,000 a.u/cell as with previous steps.

Figure 11 presents the concentration of mKate as a function of the concentration of dimerized NarX. At above 100 units of dimerized NarX, the concentration of NarX starts to saturate and stabilize at about 35,000 units.

Integrated Model

The integrated model constitutes one of the first attempts at modeling a biosensor spanning all the known reaction chain, from biomarker concentration to output while also including the induced thresholding interactions of a mutant. With this model is possible to directly relate the concentration of naringenin with the fluorescent output over time. In the figure 12 ,it can be discerned the response (fluorescent output) of the system to a infinitely increasing concentration of naringenin.

As expected, the system exhibits the same behavior as most kinematic models, saturating after achieving equilibrium. It also presents similar fluorescence values to the obtained in the model of step 4.
With the integrated model it is possible to dynamically adjust any of the variables to observe its effect on the response of the system over time. This allows the model to provide insights that none of the individual subsystems can provide by themselves. For example, an interaction of interest is the effect of the concentration of NarX mutant in the saturation of the system and its fluorescent response. By testing different concentrations of NarX mutant with a steadily increasing concentration of naringenin, it is possible to calculate the timeframe required to achieve saturation.

Table 1. Concentration of mKate and time until saturation for each level of NarX mutant

Concentration of Narx Mutant	mKATE concentration	Time for Saturation
0	3.598x10^4	Achieved instantly
2500	3.598x10^4	1.7s
5000	3.598x10^4	2s
10000	3.598x10^4	3.3s
12500	3.598x10^4	3.7s
15000	3.598x10^4	4s
20000	3.598x10^4	5s

As it can be noted in the table 1, the time required for the system to achieve saturation is directly proportional to the concentration of NarX mutant. In contrast, the mKATE concentration is inversely proportional to it. In both cases however, the relationship appears to be linear.

Discussion

The models of the 4 subsystems and the final model allowed a window to characterize a complex process. A drawback of this type of model is that we have used too many assumptions and a lack of reference for the parameter. In this discussion section, the model in each step was evaluated independently and some of the assumptions were backed up by the integrated model.

Step1

The results suggest that modeling this reaction using a hill equation is adequate. However, experimental data obtained from a sensor with a mutant is required to properly assess if the similarity holds

Step2

The existence of an optimal concentration of NarX has interesting implications for the potential uses of this mechanism. It is possible to ensure saturation (i.e., drug delivery or desired response) regardless of the amount of target biomarker.

Integrated Model

Many things were assumed constant and derived from external data of system different to ours (the Ghent paper doesn't even have mKATE). With proper data it would be possible to better tailor the system and make it more robust.
Additionally, the model allows to dynamically test different presets and observe the effect on the system of variables hard to control experimentally. More so, the model allows to deduce data that cannot be directly measured.
The time capabilities allow us to predict behaviors ahead of time and evaluate things like the system stability which could be crucial in biosafety.
Due to the tight schedule of iGEM, more test with other variations weren’t possible.

Conclusion and Future Work

As with every biological model, the system necessitates assumptions and other limitations to lead to something, with further experimental data it would be possible to revise and improve the model. Future experiments can assess a more accurate approximation of the behavior of the reactions involved in this model.
The list of the assumptions is shown in following table:

Table 2. List of assumptions used for the models

Step	Assumption	Discussion
step1	NarX mutant disregarded	The hill function for NarX and NarX mutant may be a convenient approximation, could be not suitable, requires further testing
Step2	All NarXs bind to another NarX (dimerization) or NarXmutant instantaneously.	This doesn’t happen in nature, perfect dimerization never occurs. With further experimental data it could be established the rate of dimerization
Step2	The equilibrium is reached instantaneously, which means that [NarX>sub>dim], [NarXNarXmut], [NarXmutdim] will become constant instantaneously or their derivative over time will be equal.	The integrated model showed that the saturation doesn’t happen instantly, it can take up to 5 seconds for saturation depending on the concentration
Step2	All 3 reactions happen with the same association and dissociation rate which means a3 = a4 = a5 = a and b3 = b4 = b5 = b.	This is obviously not the case in nature.
Step3	the amount of Nitrate is much larger than the amount of NarL.	Wetlab team saturated the system with nitrates so NarX only depends on Naringenin, however, in vivo nitrates will not be saturated. This means that for medical implementation sensor proteins that are only depend on the target biomarker should be preferred
Step3	All NarXdims will bind to NarL.	Not the case in nature, requires experimental data to establish the actual binding rate
Step3	The equilibrium is reached instantaneously.	The integrated model showed that is not the case

The setup of the integrated model allows for future teams to focus in specific parts of the system while having access to the whole system, which can assist with guiding experimental efforts while also improving the robustness of the model. Due to lack of time, more tests were not possible. With more time, further work would have been possible.