IGEM:Imperial/2010/Mass Action

Model based on Law of Mass Action (Weeks 6 - 8)

Michaelis Menten kinetics does not apply

We cannot use Michaelis-Menten kinetics because of its preliminary assumptions, which our system does not fulfil. These assumptions are:

• [math]\displaystyle{ V_{max} }[/math] is proportional to the overall concentration of the enzyme.

But we are producing enzyme, so [math]\displaystyle{ V_{max} }[/math] will change! Therefore, the conservation [math]\displaystyle{ E_0 = E + E_S }[/math] does not hold for our system.

• Substrate >> Enzyme.

Since we are producing both substrate and enzyme, we have roughly the same amount of substrate and enzyme.

• Enzyme affinity to substrate has to be high.

Therefore, the model above is not representative of the enzymatic reaction. As we cannot use the Michaelis-Menten model we will have to solve from first principle (which just means writing down all of the biochemical equations and solving for these in Matlab).

Change of output

Instead of GFP, it is now dioxygenase acting on catechol (activating it into colourful form). Catechol will be added to the bacteria manually (i.e. the bacteria will not produce Catechol). Hence, in our models dioxynase is going to be treated as an output as this enzyme is the only activator of catechol in our system. This means that the change of catechol into its colourful form is dependent on dioxygenase concentration.

Models

Model preA: Simple Production of Dioxygenase

Simple production of Dioxygenase

This model includes transcription and translation of the dioxygenase. It does not involve any amplification steps. It is our control model against which we will be comparing the results of other models.

• Here is the Matlab code for Model pre-A.

Model A: Activation of Dioxygenase by TEV enzyme

1-step amplification

The reaction can be rewritten as:

[math]\displaystyle{ {TEV} + {split Dioxygenase} \rightleftarrows {TsD} \rightarrow {TEV} + {Dioxygenase} }[/math]

This is a simple enzymatic reaction, where TEV is the enzyme, Dioxygenase the product and split Dioxygenase the substrate. Choosing [math]\displaystyle{ k_1, k_2, k_3 }[/math] as reaction constants, the reaction can be rewritten in these four sub-equations:

1. [math]\displaystyle{ \dot{T} = -k_1[T][sD] + (k_2+k_3)[TsD] + s_T - d_T[T] }[/math]
2. [math]\displaystyle{ \dot{sD}= -k_1[T][sD] + k_2[TsD] + s_{sD} - d_{sD}[sD] }[/math]
3. [math]\displaystyle{ \dot{TsD} = k_1[T][sD] - (k_2+k_3)[TsD] - d_{TsD}[TsD] }[/math]
4. [math]\displaystyle{ \dot{D} = k_3[TsD] - d_D[D] }[/math]

These four equations were implemented in Matlab, using a built-in function (ode45) which solves ordinary differential equations.

• Here is the Matlab code for Model A.

Implementation in TinkerCell

Another approach to model the amplification module would be to implement it in a program such as TinkerCell (or CellDesigner). This would be useful to check whether the Matlab model works.

Network implemented in TinkerCell

Model B: Activation of Dioxygenase by TEV or activated split TEV enzyme

2-step amplification

This version includes the following features:

• 2 amplification steps (TEV and split TEV)
• Split TEV is specified to have a and b parts
• TEVa is forbidden to interact with TEVa (though in reality there could be some affinity between the two). Same for the interaction between Tevb and Tevb
• Both TEV and TEVs are allowed to activate dioxygenase
• Dioxygenase is assumed to be active as a monomer
• Activate split TEV (TEVs) is not allowed to activate sTEVa or sTEVb (this kind of interaction is accounted for in the next model version)
• This model does not include any specific terms for time delays
• Here is the Matlab code for Model B.

Model C: Further improvement

This model has not been implemented because of the results from Models A and B.

This version adds the following features:

• activated split TEV (TEVs) is allowed to activate not only sD but sTEVa and sTEVb

Results

The first results that were obtained seemed to be flawed since they indicated negative concentrations would be obtained from the amplification step. In particular, for concentrations smaller than [math]\displaystyle{ 10^{-4} mol/dm^3 }[/math] the results were inconclusive since they were oscillating around zero. We realised that this could be due to the ode-solver that we were using (ode45 in Matlab).

Trying to correct this problem with the ode-solver, the following precuations were implemented:

• NonNegative function in Matlab preventing solver from reaching negative values - still some marginally negative values show
• Scaling - all the values were scaled up by a factor of [math]\displaystyle{ 10^6 }[/math] as working on small numbers could be problematic for Matlab. Once the result is generated by the solver the resulting matrix is scaled back down by [math]\displaystyle{ 10^6 }[/math].

When we entered the real production and degradation rates into our model, we once again obtained nagetive values. This was due to our set of differential equations being stiff. Since ode45 cannot solve stiff differential equations, we had to switch to using ode15s - an ode-solver designed to handle stiff equations.

Model pre-A

This is the result of the simulation of simple production of Dioxygenase. It can be seen that the concentration will tend towards a final value of approximately [math]\displaystyle{ 8 }[/math]×[math]\displaystyle{ 10^{-6} mol/dm^3 }[/math]. This final value is dependent on the production rate (which has been estimated for all of the models).

Results of the Matlab simulation of Model preA

Model A

• Initial Concentration

The initial concentration of split Dioxygenase, c0, determines whether the system is amplifying. The minimum concentration for any amplification to happen is [math]\displaystyle{ 10^{-5} mol/dm^3 }[/math]. If the initial concentration of split Dioxygenase is higher, then the final concentration of Dioxygenase will be higher as well (see graphs below). Note that the obtained threshold value is above the maximum value that can be generated in the cell according to Model pre-A!!!

Comparison between Model pre-A and Model A. Initial concentration of split Dioxygenase: [math]\displaystyle{ 10^{-5} mol/dm^3 }[/math]

• Changing Km:

Km is indirectly proportional to the "final concentration" (which is the concentration at the end of the simulation), i.e. the bigger the evalue of Km, the smaller the "final concentration" will be. Different Km values determine how quickly the amplification will take place.

(Also, it was found that the absolute value of [math]\displaystyle{ k_1 }[/math] and [math]\displaystyle{ k_2 }[/math] entered into Matlab does not change the outcome as long as the ratio between them ([math]\displaystyle{ K_m }[/math]≈[math]\displaystyle{ k_2/k_1 }[/math]) is kept constant. This is important when simulating (in case entering very high values for [math]\displaystyle{ k_1 }[/math] and [math]\displaystyle{ k_2 }[/math] takes too long to simulate).

• Changing kcat

Model predicts the concentration values for dioxygenase to raise quicker and to higher values for increasing values of [math]\displaystyle{ k_{cat} = k_3 }[/math]. That indicates that [math]\displaystyle{ k_3 }[/math] is actually the slowest step in enzymatic reaction and allows us to appreciate how our system is dependent on kinetic properties of enzyme.

• Changing production rate

At the moment, our biggest source of error could be the production rate, which we weren't able to obtain from literature. So, we had to estimate (see below) the value of the production rate. We hope to be able to take rough measurements of that value in the lab as it has big effect on models' behaviour.

Sensitivity of Model A (20/08/2010)

We want to determine how our system reacts if different parameters are changed. This is to find out which parameters our system is very sensitive to.

Parameter	Sensitivity
Initial concentration of split Dioxygenase	Change of one order of magnitude in the initial concentration, [math]\displaystyle{ c_0 }[/math], gives change of oen order of magnitude in the output concentration (range: [math]\displaystyle{ 1\gt c_0\gt 10^{-5} }[/math]). Loses sensitivity for extremely high or low values.
[math]\displaystyle{ K_m }[/math]	Change of one order of magnitude results in change of output concentration by one order of magnitude ([math]\displaystyle{ 0.01\lt K_m\lt 100 }[/math]). At values smaller than 0.01, the sensitivity is lost. For higher values than 100 the sensitivity is at least the same as the change of order of magnitude.
[math]\displaystyle{ k_{cat} }[/math]	[math]\displaystyle{ k_{cat} }[/math] proportional to dioxygenase production (1-to-1 sensitivity for all values) for an initial concentratio of [math]\displaystyle{ 0.01 mol/dm^3 }[/math]. For very high initial concentrations, the system is very sensitive to changes in [math]\displaystyle{ k_{cat} }[/math].
Production rate of TEV	1-1 sensitivity for most values. At some point the system’s response is limited by the initial concentration of sD, so for very high TEV production rates not much change is observed.
Production rate of split Dioxygenase	Not much influence on 1-step amplification. However, the value seems to be crucial for simple production of Dioxygebase (1-1 order of magnitude sensitivity).
Degradation rates	Sensitive within the relevant range. Not very sensitive for values smaller than [math]\displaystyle{ 10^{-6} }[/math]. For high degradation rates (1>degradation rate>0.01): unexplainable behaviours.

Hence, the system is sensitive to most of the constants (given a particular range of values). The most crucial one, however, seems to be the initial concentration of split Dioxygenase.

Model B

• Initial Concentration

The initial concentration of split Dioxygenase, [math]\displaystyle{ c_0 }[/math], determines whether the system is amplifying. .

The behaviour in varying the initial concentration obeys similar relationship as the one of Model A: If the initial concentration of split Dioxygenase is higher, then the final concentration of Dioxygenase will be higher as well (see graphs below).

• Model A vs. B

Having run both models with the same initial conditions ([math]\displaystyle{ c_0=10^{-5} mol/dm^3 }[/math]). It has been noted that Model B does not generate very siginificant amplfication over the Model A. Hence, it would be more sensible to integrate a one step amplification module.

Comparison between Models pre-A, A and B

Colour response model (25/08/2010)

Initially, dioxygenase was being treated as "output" just because we decided to model what is different between the models in order to determine which one is better. We have found that 2 step amplification (Model B) added to a system presented little improvement over the 1 step amplification (Model A). This conclusion was worrying as action of dioxygenase on catechol to produce colour compound is enzymatic, hence a amplification step by itself. Adding it to all models meant that Model A becomes a 2 step amplifier and Model B becomes a 3 step amplifier. This was worrying having previously drawn the conclusion that 2 step amplifier is not much better than 1 step one. We decided to model that in order to check whether that deduction was true. If it was true that would mean that our construct is not innovative at all.

The important information about colour production is that the coloured compound kills the cells. It is not catechol killing them but the coloured compound.

It is suggested that product of Catechol destroys the cell membrane by inhibiting lipid peroxidation. It causes significant changes in the structure and functioning of membrane components (e.g. disruption of membrane potential, removal of lipids and proteins, loss of magnesium and calcium ions). These effects cause the loss of membrane functions, leading to cell death.

Since the product of Catechol acts on the cell membrane, it will not affect our enzymatic reaction immediately. In our simulation, we will try to model immediate cell death as well as neglect the effect that Catechol has on the cell and compare the two if there is a significant difference in results. If the differences will be appreicable we will try to model the slow cell death in more detail.

• Toxicity of three phenolic compounds and their mixtures on the gram-positive bacteria Bacillus subtilis in the aquatic environment
• Toxic Effects of Catechol and 4-Chlorobenzoate Stresseson Bacterial Cells

Initial conclusions(26/08/2010)

Despite our model not working entirely correctly (the simulator goes a bit funny), we deduced several points.

• The images presented below show cathecol being added at 3 different points in time. Cross section refers to a point in time at which concentration of dioxygenase in amplified systems crosses to be above the concentration of the non-amplified system. From those graphs it is clear that amplification in output is visible only after the cross-section has been reached. Note that those simulations were run for 1M solutions of catechol (which is quite high) which allows to see the differences between various amplification models easily.

Concentration of coloured compound for catechol being added before the cross section is reached

Concentration of coloured compound for catechol being added when the cross section is reached

Concentration of coloured compound for catechol being added after the cross section has been reached

• We noticed that amplification has point only for quite high catechol initial concentrations (>0.01M). For smaller concentrations of catechol the dioxygenase conetrations in different system do not seem to be decisive about the speed of response (no difference between all 3 models). That means that basically, in systems with small catechol concentration added the amplified systems end up being redundant (dioxygenase is overproduced) as concnetration of dioxygenase from simple production seems to high enough to convert catechol almost instantenously. Amplification models become only meaningful when they have a lot of substrate to act on (ie. high concentration of catechol). This leads us back to determination of colour compound concentration threshold for visibility. It will a crucial factor in deciding whether amplifiers designed by our team obtain the fast response or not.

Problems

Once we have implmented the colour change into models, we have noticed that there is a problem with it. After adding catechol, some concentrations were jumping to negative values, which was worrying. We checked our equations and constants but could not find the mistake. Hence, we concluded that there is something wrong with the way that MatLab evaluates the equation or treats numbers. The trouble seemed to have origined from the very rapid concentration change of catechol which disrupts the whole system Prospective solutions to problem:

• Implementation in TinkerCell (31/08/2010)

We hoped TinkerCell will have non-negative conditions imposed on its solutions and will do better job than raw solver manipulating the equations. Hence, we implemented the whole amplification model (including coloured output) in TinkerCell. However, we realized that TinkerCell does not deal very well with very high or low numbers (For example, values higher than [math]\displaystyle{ 10^5 }[/math] are not acceptable - this is important since our rate constants ([math]\displaystyle{ k_1 }[/math]) are usually bigger than [math]\displaystyle{ 10^5 }[/math]. Also, the low degradation rates ([math]\displaystyle{ 10^{-9} }[/math]) result in a zero output line). However, TinkerCell can still be used for testing that our Matlab programs behave the way we anticipated (by using default values of 1), as well as producing diagrams of our system.

• Varying ODE solver options in MatLab (31/08/2010)

We had a close look at the ODE solver options in MatLab. We those that could have helped us. In the end we found that decreasing Relative and Aboslute tolerances (to values as small as [math]\displaystyle{ 10^{-15} }[/math]) was significantly improving the situation. However, this is not an ultimate solution as in simulations still negative numbers appear of order of [math]\displaystyle{ 10^{-15} }[/math]. As no better solution than that was found, we dicded to accept that approach. We decided as well that we are just interested in first 100 to 150 seconds after adding catechol, while concentrations hit the negatve values way after that. The species that were concerned with hitting negative values were approaching to 0.

THe images below show the influence of the relative and absolute tolerance values on the model. Note it was important to allow ODE solver to adjust the time step automatically, as big time steps (1 second) were generating false answers for the catechol model. Adjustign manually time steps to really small values was not efficient (the whole simulation does not require very high definition simluation). Function of automatic choice of step showed to be able to cope well with that issue.

Graphs representing conctrations of catechol, dioxygenase, catechol-dioxygenase complex and colur compound when adding catechol at t=2000s. Note the really negative values that concentration of catechol hits.

Here the Relative and the Absolute Tolerances are lowers to 1E-18 and 1E-12, respectively. From graphs it is not clear that any of the concentration parameters hits negative values. However, after having a look at raw data it is clear that catechol and catechol-dioxygenase complex get slightly negative values like -2E-120

• Using SimBiology to model (31/08/2010)

We again hoped that SimBiology package could be more suited to mdelling our situation then raw ODE solvers, so we implemented our some models in it. This package offers interactive user interface similar to Tinker Cell, but uses MatLab to simulate. Initially, we have confirmed that our simple production model (Model PreA) and 1 step amplification model (Model A) implemented in visual manner in MatLab generated exactly the same results as our manually derived ODE equation based models. The interface allowed us to have clearer control over paramters. As well it allowed modelling special events, like adding catechol at certain point in time. Previously we had to split simulation into to parts.

Final Conclusions (31/08/2010)

The conclusions obtained after fixing the error of negative concentrations in the model, did not change our initial conclusions. However, they are stated here clearly:

• Changing time when catechol is added

If Catechol is added before t= 1000s, then the coloured output will reach its threshold value faster by simple production. If Catechol is added when t>1000s, then the coloured output will increase (marginally) faster through the amplification step in Model A. This does not seem to be appreciable difference between the two models in given time scale. These observations are true for intial concentration of dioxygenase equal to [math]\displaystyle{ 10^{-5}mol/dm^3 }[/math]. However, we noticed that if the initial concentration is raised to [math]\displaystyle{ 10^{-4}mol/dm^3 }[/math], then Model A can be advategous over Model preA within as little as 100 seconds.

• Changing concentration of catechol added

There seem to be 3 regions of catechol concentration that influence the system in different ways. Those regions are: c>1M, 1M>c>0.01M, 0.01M>c. The boundatries of these regions tend to vary depending on the choice of other initial conditions. The values given above apply to boundary conditions that currently are considered physiologically relevant. Varying initial concentration of catechol within the highest region doesn't result in any change in colour output response (probably all enzymes are occupied, the solution is over saturated with catechol). In the middle region the catechol ocncentration has influence on the amplficaytion. Amplification becomes less and less towards c=0.01M. When 3 region is entered, there is no difference in output production by the two models

• Cell death

The coloured product of catechol kills cells by destroying the cell membrane. However, we don't know how quickly the cells will die. Therefore, we examined two different cases: immediate cell death and negligible cell death (i.e. cells death is negligible because it takes too long)

Running the simluation in Matlab (not Simbiology!), our conclusions are:

1. Immediate cell death slows down production of coloured output. Depending on the threshold concentration this can delay the detectable response by a few minutes.
2. If Catechol is added before t=1000s, then cell death slows down the response considerably.
3. In case of cells being modelled as alife, there difference between the amplified and the simple production model is smaller than it is in case of cells dying.