IGEM:Imperial/2010/Modelling

Have a look at this link: Synthetic Biology (Spring2008): Computer Modelling Practicals

Have a look at Cell Designer to easily generate images of the system.

Example on how Valencia 2006 team used SimulLink to simulate their project: Valencia 2006 PowerPoint presentation

Objectives

Week 6

Day	Monday	Tuesday	Wednesday	Thursday	Friday	Weekend
Date	09	10	11	12	13	14-15
Objective				Find constants	Find protein production constants and TEV reaction rate constants
Completion				We didn't manage to complete the task	The orders of magnitude established - ready to run simulations

Week 7

Day	Monday	Tuesday	Wednesday	Thursday	Friday	Weekend
Date	16	17	18	19	20	21-22
Objective	Implementing the constant ranges in the output model. Comparing the results between the models.	Start modelling the protein display signalling to find the concentrations. Explain the oscillations that are occuring in the output amplification model	Research stiff differential equations Research on receptor (especially MAP kinase) Finalise results for amplification model Prepare presentation	During meeting it was decided that we should test amplificiation model for its sensitivity to change of parameter values. This is to determine the key parameters to be measured in the lab.	Discuss the experiments for amp. module further Start on modelling the protein display on cell surface
Completion	The task is accomplished. However, unexplained oscillations are observed for some specific values	We got rid of the oscillations in our model (by using ode15s instead of ode45).	Presentation ready Reading on stiff equations done	Key sensitive parameters determined Experiments initially discussed with Chris	Further consult on experiments for determining parameters needed for modelling with James and Wolf Initial thought for a protein display model implemented in MatLab

Week 8

Day	Monday	Tuesday	Wednesday	Thursday	Friday	Weekend
Date	23	24	25	26	27	28-29
Objective	Test the protein display model and find to what kind of values it's sensitive to	define the control volume around the bacterial cell Start on testing the model	We realised that it would be worth adding the last step (colour production by dioxygenase) into our amplfication models	We need to finalise testing and formulation of certain assumptions regarding the amplfication (cell death in particular) Consider helping RMIT- Australia in modelling wiki	Keep exploring the prospective solutions to our system going unstable Try to implement the last bit of amplification in TinkerCell. Maybe it will deal with it well.
Completion	We couldn't reach the testing as the first model "Display 1" was simulating reactions as if the were taking place inside the bacterial cell while in reality the take place outside. Then we hit an issue of defining the size of volume around the cell. We didn't manage to resolve that problem today.	Control volume defined Testing of model done Protein model pretty much finalised	The colour production almost complete, however complications were encountered regarding the cell death (colour compound slowly kills cells)	Basically we're stuck: scaling didn't help neither did trying idfferent solvers Meeting with Matthieu gave us some prospective solutions that we can explore Alternatively we can have a look at SimBiology or MatCont helping RMIT - their request is for protein engineering simulation - we cannot do that	We know how to use the splinetool in Matlab - it's tool to smoothen noisy graphs allows to manipulate graphs We found that specifying time span (initial and final time) instead of specific time array with defined time steps helps to reduce the simulation time, but it does not improve much on issue we're trying to tackle at the moment

Notebook

Output Amplification Model

Motivation: We have come up with a simple concept of amplification of output done by enzymes. Before the final constructs are assembled within the bacterial ogranism, it is beneficial for us to model the behaviour of our design.

The questions that need to be answered:

1. How beneficial is the use of amplification? (compare speeds of response of transcription based output to amplified outputs)
2. How many amplification steps are beneficial to have? (if too many amplification steps are involved, the associated time delay with expressing even amplfiied output may prove not to be beneficial.)
3. Does mixing of amplfication levels have a negative influence on the output? Is it better to use TEV all the way or HIV1? Modelling should allows us to make a decision on which design is more efficient.

First attempt

At each stage of amplification a distinct protease is being used

At each stage of amplification a distinct protease is being used

---

TEV is used at both stages of amplification

TEV is used at both stages of amplification

Second attempt

Model improved to account for the enzymes (protease action)

Implementation in Matlab

The Matlab code for the different stages of amplification and diagrams can be found here.

Kinetic constants

Quality	GFP	TEV	split TEV	split GFP
Km and Kcat	Doesn't apply	TEV constants (Km and kcat)	40% of whole TEV	Doesn't apply
half-life or degradation rate	Half-life of GFP in Bacillus = 1.5 hours - ref. Chris	?	?	Half-life shorter than GFP
production rate in B.sub	?	?	?	?

Conclusions

We couldn't obtain all the necessary constants. Hence, we decided to make educated guesses about possible relative values between the constants as well as varying them and observing the change in output.

As the result, we concluded that the amplification happens at each amplification level proposed. It's magnitude varies depending on the constants. There doesn’t seem to be much difference in substitution of TEV with HIV1.

Output Amplification Model - Modified version

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:

• Vmax is proportional to the overall concentration of the enzyme.

But we are producing enzyme, so Vmax will change! Therefore, the conservation E0 = E + ES 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).

Changes in the system

GFP is not used any more as an output. It is dioxygenase acting on the catechol (activating it into colourful form). Catechol will be added to bacteria, it won't be produced by them. Hence, basically in our models dioxynase is going to be treated as an output as this enzyme is recognised as the only activator of catechol in our system. This means that change of catechol into colourful form is dependent on dioxygenase concentration.

Models:

Model preA: 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.

Model A: Activation of Dioxygenase by TEV enzyme

The reaction can be rewritten as: TEV + split Dioxygenase <-> TEV-split Dioxygenase -> TEV + Dioxygenase. This is a simple enzymatic reaction, where TEV is the enzyme, Dioxygenase the product and split Dioxygenase the substrate. Choosing k1, k2, k3 as reaction constants, the reaction can be rewritten in these four sub-equations:

1. [T'] = -k1[T][sD] + (k2+k3)[TsD] + sT - dT[T]
2. [sD']= -k1[T][sD] + k2[TsD] + ssD - dsD[sD]
3. [TsD'] = k1[T][sD] - (k2+k3)[TsD] - dTsD[TsD]
4. [D'] = k3[TsD] - dD[D]

These four equations were implemented in Matlab, using a built-in function (ode45) which solves ordinary differential equations. The Matlab code for this module can be found here.

Results of the Matlab simulation, setting all constants to 1

Implementation in TinkerCell

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

LHS: Network implemented in TinkerCell, RHS: constants and results

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

1. T + sTa + sTb + sD <-> T-sD --> T + D + sTa + sTb
2. T + D + sTa + sTb <-> T-sTa --> T + Tsa + D + sTb
3. T + D + Tsa + sTb <-> T-sTb --> T + Tsb + Tsa + D
4. T + Tsa + Tsb + D <-> Tsa-Tsb --> Ts + T + D
5. T + sTa + sTb + sD -- T-sTa --> T + Tsa + sTb + sD
6. T + Tsa + sTb + sD -- T-sD --> T + Tsa + sTb + D
7. T + Tsa + sTb + D -- T-sTb --> T + Tsa + Tsb + D
8. T + Tsa + Tsb + D -- Tsa-Tsb --> T + Ts + D
9. T + sTa + sTb + sD -- T-sTb --> T + Tsb + sTa + sD
10. T + Tsb + sTa + sD -- T-sD --> T + Tsb + sTa + D
11. T + Tsb + sTa + D -- T-sTa --> T + Tsb + Tsa + D
12. T + Tsb + Tsa + D -- Tsa-Tsb --> T + Ts + D
13. T + sTa + sTb + sD -- T-sTa --> T + Tsa + sTb + sD
14. Ts + T + sD -- Ts-sD --> Ts + T + D

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 TEVa (though in reality there could be some affinity between the two). Same for interaction between Tevb and Tevb
• Both TEV and TEVs are allowed to activate dioxugenase molecule
• Dioxugenase 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)
• There is no specific terms for time delays included

The MatLab code can be found here. Note that no final conclusions should be drawn before realistic estimates for kinetic constants are included. It wasn’t done so far.

All chemical species appearing in the model	Network of the improved model
Resulting graphs part 1. Compare the production graphs of TEV (transcribed and translated from scratch and the Dioxugenase which is the final species in the whole cascade	Resulting graphs part 2.

Model C: Further improvement

This model is not implemented yet.

This version adds the following features:

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

Network of the further improved model

Network of the further improved model (continued)

Output Amplification Model - Results

The major concern of the results that we got (in particular for small concentrations < 10^(-4) mol*dm^(-3)) was that the solver was oscillating about positive or zero values but marking concentration values below zero too. It was recognised as faulty and probably leading the solver to false solutions.

Trying to prevent the ode solver going crazy, the following precuations were implemented:

• function preventing solver from going to negative values (does it really work) - still some marginally negative values show
• Scaling - all the values were scaled up by a factor of 10^6 as working on small numbers could be problematic for matlab (especially at the beggining of the cascade). Once the result is generated by the solver the resulting matrix is scaled back down by 10^6

Model pre-A

This is the result for the simulation of simple production of Dioxygenase. It can be seen that the concentration will tend towards a final value of around 8*10^-6 mol dm^(-3). This final value is dependent on the production rate (which we have just estimated!).

Results of the Matlab simulation of Model preA

Model A

When we entered production and degradation rates into our model, it did not seem to work properly (e.g. we got negative concentrations as our output). We found out that this is due to our set of differential equations being stiff. Since ode45 cannot solve stiff differential equations, we had to switch to using ode15s.

• Initial Concentration

The initial concentration of split Dioxygenase, c0, determines whether the system is amplifying. The minimum concentration for any amplification to happen is 10^-5 mol dm^(-3). 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: 10^-5 mol dm^(-3)

• 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 k1 and k2 entered into Matlab does not change the outcome as long as the ratio between them (Km~k2/k1) is kept constant. This is important when simulating (in case entering very high values for k1 and k2 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 kcat=k3. That indicates that k3 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, c0, gives change of oen order of magnitude in the output concentration (range: 1>c0>10^-5). Loses sensitivity for extremely high or low values.
Km	Change of one order of magnitude results in change of output concentration by one order of magnitude (0.01<Km<100). 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.
kcat	kcat proportional to dioxygenase production (1-to-1 sensitivity for all values) for an initial concentratio of 0.01 mol/dm^3. For very high initial concentrations, the system is very sensitive to changes in kcat.
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 10^-6;. 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, c0, 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 (c0=10^-5 mol dm^(-3)). 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.

Constants for Modelling

Type of constant	Derivation of value
TEV Enzyme dynamics	Enzymatic Reaction: E + S <-> ES -> E + P Let k1 = rate constant for E + S -> ES k2 = rate constant for E + S <- ES kcat = rate constant for ES -> E + P We know that Km = (kcat + k2)/k1 Assuming that kcat << k2 << k1, we can rewrite Km ~ k2/k1 From this paper constants for TEV can be found: e.g. wildtype TEV Km = 0.061 +/- 0.010 mM kcat = 0.16 +/- 0.01 s^-1 These values correspond with our assumption that kcat ~ 0.1 s^-1 and Km ~ 0.01 mM. Hence, we can estimate the following orders of magnitude for the rate constants: k1 ~ 10^8 M^-1 s^-1 k2 ~ 10^3 s^-1 Using these values should be a good approximation for our model.
Degradation rate (common for all)	Assumption: To be approximated by cell division (dilution of media) as none of the proteins are involved in any active degradation pathways Growth rate (divisions/h): 0.53<G.r.<2.18 Hence on average, g.r.=1.5 divisions per hour => 1 division every 1/1.5=0.6667 of an hour (40mins) To deduce degradation rate use the following formula: τ_(1⁄2)=ln2/k Where τ_(1⁄2)=0.667 hour, k – is the degradation rate k=ln2/τ_(1⁄2) = 0.000289 s^(-1)
Production rate (TEV and dioxygenase)	We had real trouble finding the production rate values in the literature and we hope to be able to perform experiments to obtain those values for (TEV protease and catechol 2,3-dioxygenase). The experiments will not be possible to be carried out soon, so for the time being we suggest very simplistic approach for estimating production rates. We have found production rates for two arbitrary proteins in E.Coli. We want to get estimates of production rates by comparing the lengths of the proteins (number of amino-acids). As this approach is very vague, it is important to realise its limitations and inconsistencies: Found values are taken from E.Coli not B.sub. The two found rates are of the same value for quite different amino-acid number which indicates that protein folding is limiting the production rates (we use the chosen approach as the only way of getting the estimate of order of reaction) LacY production = 100 molecules/min (417 Amino Acids) LacZ production = 100 molecules/min (1024 AA) Average production ≈ 100molecules/min 720 AA That gives us: TEV production ≈ 24 molecules/min = 0.40 molecules/s (3054 AA) As production rate needs to be expressed in concentration units per unit volume, the above number is converted to mols/s and divided by the cell volume → 2.3808*10^(-10) mol*dm^(-3)*s^(-1) C23D production ≈ 252 molecules/min = 4.2 molecules/s (285 AA) → 2.4998*10^(-9) mol*dm^(-3)*s^(-1) We will treat these numbers as guiding us in terms of range of orders of magnitudes. We will try to run our models for variety of values and determine system’s limitations.
Kinetic parameters of dioxygenase	Initial velocity of the enzymatic reaction was investigated at pH 7.5 and 30 °C. Wild type (we use that one) Km = 10μM kcat = 52 s^(−1) Mutated type Km = 40μM kcat = 192 s^(−1) Consequently, the kcat/Km 4.8 of the mutant was slightly lower than kcat/Km 5.2 of the wild type, indicating that the mutation has little effect on catalytic efficiency. reference
Dimensions of Bacillus subtillis cell	Dimensions of Basillus subtilis (cylinder/rod shape) in reach media (ref. bionumbers): diameter: d= 0.87um length: l=4.7 um in rich media This gives: Volume= π∙(d/2)^2∙l=2.793999 μm^3≈ 2.79∙10^(-15) dm^3
Split TEV production rates	*Assume the both parts of split TEV are half of size of the whole TEV (3054/2=1527 AA) The length of the coil is 90 AA The whole construct is then: 1617 AA → split TEV production rate ≈ 1.2606*10^(-10) mol*dm^(-3)*s^(-1)
Relevant concentrations of catechol	We have catechol in the lab in powder form so we're limited only by catechol's solubility. For concentration of 0.1 M with built up levels of dioxygenase the colour change happens within seconds! We will run our models for 0.1M ± several orders of magnitude to determine the smallest catechol concentration that will give appreciable difference between simple production response and the amplified response.

Receptor and Surface protein model

The aim of this model is to determine the concentration of Schistosoma elastase or TEV protease that should be added to bacteria to trigger the response. It is also attempted to model how long does it take for protease or elastase to cleave enough peptides.

Cleavage of protein is an enzymatic reaction, which can be written as:

• S + E <--> E-S --> P + E
• Substrate (S) = Protein
• Enzyme (E) = TEV (Protease)
• Product (P) = Peptide

This can be modelled in a very similar way to the 1-step amplification model, however, all the constants and initial concentrations will be different.

[TEV](t=0) - initial concentration of TEV will be arbitrary for us to choose. However, ultimately we would need to measure the concentration of elastase the schistosoma releases.

Threshold concentration of peptide (20/08/2010)

The peptide concentration required to activate ComD is 10 ng/ml (see this paper). This is the threshold value for ComD. We want to know how long it takes until the threshold will be reached.

• The mass of a peptide is 2.24kDa = 3.7184*10^-21g.
• The number of molecules in a ml is 10ng/3.7184*10^-21g = 2.6893*10^12. In a litre, the number of molecules is 2.6893*10^15.
• Dividing this value by Avogadro's constant gives the threshold concentration of c_th=4.4658*10^-9 mol/L.

Protein production in B.sub (23/08/2010)

Reference for protein production (Page 4, Paragraph 1)

• The paper mentions that each cell displays 2.4*10^5 peptides.
• 2.4*10^5 molecules = 2.4*10^5/6.02*10^23 mol = 0.398671*10^-18 mol
• Volume of B.sub: 2.79*10^-15 dm^3

• Concentration = [mol/L]
• c = 1.4289*10^-4 mol*dm^(-3)
• This is the concentration of protein that will be produced in a single cell of B.sub.

Hence, we can deduce the concentration that the protein expression will tend to (c = 1.4289*10^-4 mol/dm^3 = c_final). Therefore, we can model the protein production by transcription and translation and adjust the production constant so the concentration value will converge to c_final.

Using a similar model to the simple production of Dioxygenase for the Output Amplification Model (Model preA), we obtain the following graph:

Production of protein by transcription and translation

The degradation rate was kept constant, and the production rate was changed according to the final concentration.

Control volume 1 - initial choice (23/08/2010)

All enzymatic reactions that we have modelled so far were happening in the fixed volume (whithin bacterial cell). This case is different as the molecules are not confined by the bacterial membrane and can diffuse freely out of the cell.

The control volume: In order to avoid trying to account for transient diffusion, we decided to determine combination of imaginary and real boundaries for the ssytem. The inner boundary is determined the bacterial cell (proteins after being displayed and cleaved cannot diffuse back into bacterium). The outer boundary is more time scale dependent. We have decided that after mass cleaving of display-proteins by TEV many of them will bind to receptors quite quickly (eg.8 seconds). Our volume is determined by the distance that AIPs could travel outwards by diffusion within that short time. In that way, we are sure that the concentration of AIPs outside our control volume within given time is approximately 0.

This approach is not very accurate and can lead us to false negative conclusions (as in reality there will be concentration gradient with highest conentration at the call wall, not evenly spread concentration).

Control volume (volume of B.sub. to be excluded. x indicates the distance travelled by AIPs from the surface in time t)

Difussion distance was calculated using Fick's 1st Law: x=sqrt(2*D*t), where: x - diffusion distance, D - diffusion constant, t - time of diffusion

D_average = 10.7*10^(-11) m^2/s - for a protein in agarose gel for pH=5.6 according to paper

t = 8s (for now chosen arbitrarily by us - we hope this is long enough time for AIP to bind to ComD)

This gives: x = 4.14*10^(-5) m

The control volume can be calculated by adding 2x to the length and the diamter of the cell. This gives a control volume (CV) = 4.81*10^(-7) m^3

Production of protein by transcription and translation in control volume

Protein production in Control Volume (23/08/2010)

The previously determined constants of protein production in B.sub to obtain the concentration of proteins per cell indicated by paper is not valid in Control Volume. It has to be adjusted (multiplied) by the following factor:

factor= V_bacillus/V_control_volume = 5.7974*10^(-6) (for particular numbers presented above)

It is enough to adjust the production rate by the 'factor'. The resulting concentration will follow.

Control volume final choice (23/08/2010)

We have realised that initial choice of control volume was not reflecting the reality well. It was treating single bacterium as if it was by itself in the medium. However, in reality bacteria live in colonies very close to each other. They are closer to each other that the diffusion distance (1.9596*10^(-5) m) derived above even if placed in water solution.

Using CFU to estimate the spacing between cells (24/08/2010)

CFU stands for Colony-forming unit. It is a measure of bacterial numbers. For liquids, CFU is measured per ml. Since we already have data of CFU/ml from the Imperial iGEM 2008, this is an easy way to estimate the number of cells in a given volume using spectrometer at 600nm wavelength. The graph below is taken from the Imperial iGEM 2008 Wiki page.

CFU/ml vs. Optical Density at 600nm (OD600)

The graph shows values of CFU/ml for different optical densities. The range of CFU/ml is therefore between 0.5*10^8 - 5*10^8.

In this calculation, we will assume that only one cell will grow and become one colony (i.e. no more than one cell will form no more than one colony). Therefore, the maximum number of cells in 1ml of solution is 5*10^8. Taking the volume of 1 ml = 10^-3 dm^3 and dividing by the (maximum) number of cells in 1ml gives the average control volume (CV) around each cell: 2*10^(-12) dm^3/cell. As this volume is 3 magnitudes bigger than the cell itself, its shape is arbitrary. For simplicity we choose it to be cubic. Taking the 3rd root of this value gives the length of one side of the control volume. Side length of CV = y = 1.26*10^(-4) dm = 1.26*10^(-5) m.

Choice of Control Volume allows simplifications (24/08/2010)

• Firstly, let us assume that on average the cells will be placed in the centre of the cube. This gives that protein after cleavage will have distance y/2 (assume the bacterium is dimensionless for a moment) to travel to get out of control volume. This is calculated to happen within 0.18s. Even if the bacterium was not in the centre, the protein will travel from one end of the cube to the other in less than a second (~0.74s). Hence, it will take between 0.18 and 0.74s for concentration of AIPs around the cell to be uniform. Noting that those time values are really small, we can approximate our model to be having uniform concentration across the volume. In that way we are understimating the value of AIP concentration right next to the cell's surface. If our model will happen not to be able to reach the threshold we will consider this approximation as one of the reasons.
• We can neglect the diffusive fluxes across the CV border (see figure below). Assuming that adjacent cells are producing the peptide at the same rate and the concentration of TEV is the same around them, then the fluxes should be of the same value giving net 0. Hence, we can neglect diffusion and have our model limited to 1 bacterium!

Figure showing 2 cells with their control volumes. It shows the reciprocal fluxes in the opposite directions. The situation is true for all 3 dimensions. In our model we concentrate on 1 cell that is surrounded by others

Matlab Simulation (24/08/2010)

• Here is the Matlab code for the Matlab simulation.

Graphs showing the simulation using [TEV]0 = 4*10^-4 mol/dm^3. The graph on the right hand-side below shows that the AIP threshold (red line) is reached after 22 s.

Sensitivity of our model (24/08/2010)

• Changing initial concentration of TEV

Whether the threshold concentration of AIP is reached is highly dependent on the initial concentration of TEV. The smallest initial concentration of TEV, [TEV]0, for which the threshold is reached is 6.0*10^-6 mol/dm^3. On the graoh below, it can be seen that the optimal [TEV]0 is a concentration higher than 10^-4 mol/dm^3, which corresponds to a threshold being reached within 1.5 minutes.

Graph showing when threshold AIP concentration is reached (for different initial TEV concentrations). Notice log-log scale.

• Changing the production rate

1 order of magnitude change in production rate results in at least 50s (50 seconds is the smallest step for 1 order change - for others is way beigger) delay of the AIP concentration riching the threshold concentration.

• Changing production rate

It influences a lot the time duration of the AIP concentration being above the threshold level. The higher it is the shorter the receptor is activated (at extreme values, AIP concentration never get to the threshold). However, it has not much influence on how fast the threshold is being reached.