Model

Abstract

This year, our team created a mathematical model to represent the series of reactions from Iron(II) ion concentrations in saliva to amilCP concentrations for our team’s system. We first had to decide whether we should model our final cell-free system (in vitro), or our in vivo lab work. An in vivo model would be advantageous, as it would allow us to compare the lab work we did to a model. However, we decided an in vitro model would be more appropriate. This is because the in vivo lab work was to test whether our construct worked, rather than finding a relationship between iron (II) ions and amilCP. We therefore felt a model to represent the relationship in an in vitro environment would better represent the final product our project aimed to create.

It should be noted that in a cell-free system, we would have to provide the FUR (Ferric Uptake Regulator) protein; in an in vivo system, the bacteria provides the FUR. Although it means we can easily increase the FUR concentration, FUR is very expensive (£130/5μg) [11]. We converted this value into something more useful. The molecular mass of FUR is 4047 [25] and the approximate volume of a paper-based system is 1µl, giving a value of 0.01372303¢ per nM.

With this in mind, these are the aims of our model:

(1) Show a relationship between initial Iron (II) ion concentration and amilCP concentration by means of a graph.

(2) Prove that our construct would or would not work in an in vitro environment in a cell-free system.

If the system does work, we then aim to:

(3a) Find an appropriate FUR concentration for an effective test, while minimising its cost.

(3b) Find an appropriate time to take a measurement.

(3c) Find an appropriate DNA concentration.

Model formation

We used mass action kinetics to model our reactions. This works on the assumption that the rate of reaction is proportional to the concentration of the reactants and products. We first needed to establish what our reactions were. They can be seen below:

(1)Dimerisation of the FUR protein.

(2) Transcription of LacI that is inhibited by the activated FUR.

(3) Translation of LacI protein.

(4) Transcription of amiCP that is inhibited by the LacI protein.

(5) Translation of amilCP protein.

If we let initial FUR concentration be a set value between tests, then more iron leads to more activated FUR, resulting in less LacI production and therefore more amilCP. In short, a lower iron concentration in saliva leads to a less intense blue, while a high iron concentration leads to a more intense blue.

We can break down our model into 9 different reactions to consider separately, then put them together once they are all defined.

For the transcription of LacI and amilCP, the transcription factors (FURactivated and LacI respectively) are not constant. We therefore cannot describe their rates of formation as a constant- we must consider the concentration of the transcription factors (TF). To do this, we use the Hill equation, Where:

<img src="

">

--θ is the fraction of transcription factors that are bound to the promoter.

--[TF] is the concentration of the transcription factor.

--Kd is the dissociation constant of the transcription factor to the promoter.

   <img src="">

--n is the Hill coefficient. This can be approximated to the number of TFs that binds to one promoter.

Chandran et al. discussed in their review “Mathematical modelling and synthetic biology” [1] that the rate of transcription, also known as the velocity (v) of a reaction, will be proportional to this fraction:

If the TF is an activator, v = kTX(θ)

If the TF is a repressor, v = kTX(1-θ)

Where KTX is the maximum transcription rate of the mRNA.

   <img src=""> <img src="">

Because both of our transcription reactions are regulated by a repressor TF, both our velocities will be equal to kTX(1-θ).

The Hill equation is actually used to describe the fraction of a macromolecule saturated by ligand as a function of the ligand concentration. It can be used to describe transcription, but we can also describe our FURactivated reaction by using the Hill equation, where [TF] becomes the iron (II) ions and KTX becomes initial FUR concentrations (FURt).

We can now consider each reaction. We will use v as an arbitrary variable, simply there to denote a rate of reaction. The r0x notation will be used later in the programming to remove clutter from the code

(0) Dimerisation of the FUR protein

   <img src="">

For the FURactivated reaction, we make a quasi-steady state assumption. This is used when a part of a system reacts much quicker than the other parts, and thus reaches completion on a much smaller timescale. While FURactivated reaction will reach equilibrium on a scale of seconds and minutes, translation and transcription takes hours, so instead of considering the velocity of the reaction, we will assume it happens instantly.

(1) Transcription of LacI mRNA (r01)

   <img src="">

(2) Degradation of LacI mRNA (r02)

   <img src="" align="middle">

Note that this rate will be negative. This is because it will be decreasing the concentration of the mRNA. We can now put them together, to form an ordinary differential equation (ODE) for the rate of mRNA production.

   <img src="" align="middle">

(3) Translation of LacI protein (r03)

   <img src="">

The translation reactions are fairly simple, as they are just affected by the concentration of the mRNA of the corresponding protein and can be modelled as follows:

(4) Degradation of LacI protein (r04)

   <img src="">

These two together form the ODE for the rate of LacI production:

   <img src="">

(5) Transcription of amilCP mRNA (r05)

   <img src="">

(6) Degradation of amilCP mRNA (r06)

   <img src="">

These come together to form an ODE for the rate of amilCP mRNA production:

   <img src="">

(7) Translation of amilCP protein (r07)

   <img src="">

(8) Degradation of amilCP protein (r08)

   <img src="">

These two form the ODE for amilCP concentration:

   <img src="">

Altogether, we have 5 reaction equations:

   <img src="" class="left">




   <img src="" class="left">




   <img src="" class="left">



   <img src="" class="left">




   <img src="" class="left">

Defining our constants

Now that we have our equations, we must now define our constants. This step is very important, but unfortunately we didn't have the resources to find the constants experimentally ourselves. So instead, a lot of effort was put into ensuring we had accurate values based on the literature we could find online. These are the values we found, along with their source (note that if multiple values or a range of values were found, we indicated the value we used in a bracket {x} after the values).

For the transcription and translation rates, we used the table from Eyal Karzbrun et al [7], as it describes the rate in a cell-free system using RNA Polymerase. Transcription is described in terms of rNTP/s (ribonucleoside tri-phosphates per second) (1.0±0.5 -we used 1.5, as it was closer to other literature values [10]). We could then work out kTX by dividing the number of rNTPs in each mRNA by 1.5s ((1128/1.5) = 752s for LacI, (669/1.5) = 446s for amilCP) to get the time take to transcribe one mRNA, then finding its reciprocal ((1/752)~1.33e-3/s for LacI, (1/446) = 2.24e-3/s) to find the transcription rate.

We then did the same process for translation. Eyal Karzbrun et al [7] also describes the rate of protein synthesis in terms of a rate of amino acid production. It describes the rate as >4 rNTP/s, but from other literature [10], their maximum is 21 rNTP/s, so we calculated a range of values and estimated that the rate would be somewhere in the middle. LacI is 360 amino acids long [13], so 0.014≤ kTL ≤ 0.058. We used 0.04 as our value. AmilCP is 223 amino acids long [14], so 0.022≤ kTL ≤ 0.094. We used 0.06 as our value.

Model Constant Table

Constant	symbol	Value	Source
Hill Coefficient (no units)	n1	2	Approximated from Hamed et al [2]
	n2	2.04	Bagg et al [3]
	n3	2	Becskei A et al [4]
Dissociation constant (nM)	Kd1	5.50E+04	Hamed et al [2]
	Kd2	4.6-71E+3 {4.6E+3}	Deng et al [5]
	Kd3	0.1-1E-3 or 800 {800}	Setty Y et al [12] and Basu S et al [6] respectively
Transcription constant (/s)	k_tx1 (LacI mRNA)	1.33E-03	[7] for rate, [16] for length
	k_tx2 (amilCP mRNA)	2.24E-03	[7] for rate, [14] for length
	k_tl1 (LacI)	0.014≤ Ktl ≤ 0.058 {0.04}	[7] for rate, [13] and [15] for length
	k_tl2 (amilCP)	0.022≤ Ktl ≤ 0.094 {0.06}	[7] for rate, [14] for length
Degradation rate (/s)	δ(LacI mRNA)	3.04E-03	Semsey et al [18]
	δ(LacI)	0	[21, 22]
	δ(amilCP mRNA)	2.20E-03	Selinger [9]
	δ(LamilCP)	0	[21, 22]
Range of Iron levels in saliva (nM)	Fe	2.6 - 61.8 e6	[19]

Solving the equations

There are fundamentally two ways to solve a series of ODEs: analytically and numerically, usually using software. We first attempted to solve analytically. This involves integrating the equations so that you get a variable as a function of time. We actually have 4 variables: FURt, Fe2+. time and amilCP. The first 3 are all independent and amilCP is dependent. Trying to solve this with 4 variables is near impossible, so we let FURt and Fe2+ be constant and let time and amilCP be our two variables. In this way, after integrating fully, we would get amilCP concentration as a function of time. An attempt was made to do this. However, it proved to be unsolvable. The final equation was a second order differential of amilCP as a function of time and amilCP mRNA and even Wolfram Alpha was unable to solve it. This being the case, we then opted to try using software, namely MatLab. Note that the k_tx and k_tl values are 0 up until a specified time. This is because in Eyal Karzbrun et al [7], they observed a delay in mRNA and protein synthesis:

“Under TX/ degredation, mRNA appeared after a delay of τ0 = 15 ± 5 … We attribute the delay to the time required for TX of a single mRNA” [Fig. 1(c)]:

“Assuming TX and TL are concurrent [20], we attribute the additional delay τf ≤ 5 min to the protein folding time, which is similar to previous measurements [21,22]” [Fig. 1(d)]

   <img src="" class="left" width="20%">

(a) Information flow from DNA (D) to mRNA (m) to protein (p) is carried by synthesis enzymes RNA polymerase Rp and ribosome R. Degradation by RNAse Xm and protease Xp.

(b) Exponential degradation of 200 nM mRNA with lifetime m ¼ 12 min.

(c) mRNA dynamics in the presence of 30 nM DNA exhibiting a delay and exponential accumulation with a time scale m.

(d) Kinetics of protein synthesis at various DNA concentrations.

(e) Normalized protein synthesis rate p_ syn=p_ max fitted to an exponential rise with a time scale m.

(f) Protein maximal accumulation rate as a function of DNA fitted with a MM curve using Rp ¼ 29 nM and KTX ¼ 1:5 nM.

These delays are not accounted for in our equations, so instead, we let k_tx and k_tl equal 0 up until the delay has ended- we estimated these to be 15 minutes for k_tx and 17 minutes for k_tl.

Because it is a cell-free system, transcription is 10 to 100 times slower than in vivo [7], which could be an issue, as the protein may reach steady state at a much lower concentration than we need. However, it should also be noted that because it is cell-free, we do not have to account for dilution of the proteins and we do not have to add protease. Since degradation of proteins is so slow in this environment [23, 24], we can essentially ignore the degradation rates of the proteins. In our code, we set d_LacI and d_CP to 0, but we could just as well have removed r04 and r08 from the ODEs.

ODE method

            <code style="overflow:visible">
timespan = [0 24*3600];  %Time span 
init = [0 0 0 0];  %Initial values of the functions       

%ODE settings
odesettings = odeset('AbsTol', 1e-12, 'RelTol', 1e-6);

%simulating the ODE
[t,f] = ode45(@odefunc, timespan, init, odesettings);

%Graphs
figure;
subplot(2,2,1);
plot(t,f(:,1));
xlabel('Time (s)');
ylabel('Concentration (nM)');
title('Concentration of LacI mRNA over time');
grid on;
grid minor;

subplot(2,2,2);
plot(t,f(:,2));
xlabel('Time (s)');
ylabel('Concentration (nM)');
title('Concentration of LacI over time');
grid on;
grid minor;

subplot(2,2,3);
plot(t,f(:,3));
xlabel('Time (s)');
ylabel('Concentration (nM)');
title('Concentration of amilCP mRNA over time');
grid on;
grid minor;


subplot(2,2,4);
plot(t,f(:,4));
xlabel('Time (s)');
ylabel('Concentration (nM)');
title('Concentration of amilCP over time');
grid on;
grid minor;


function dfdt = odefunc(t,f)
    %Constants
    if (0 < t)&&(t<900); k_TX1 = 0;   %Macimum transcription rate of LacI    
        else; k_TX1 = 1.33e-3;
    end   
    
    if (0 < t)&&(t<900); k_TX2 = 0;   %Maximum transcription rate of amilCP  
        else; k_TX2 = 2.24e-3;
    end   

    if (0 < t)&&(t<1020); k_TL1 = 0;  %Trannslation consant of LacI    
        else; k_TL1 = 0.04;
    end        
   
    if (0 < t)&&(t<1020); k_TL2 = 0;  %Translation constant of amilCP    
        else; k_TL2 = 0.06;
    end   
    
    d_mRNALacI = 3.04e-3;    %Degradation rate of LacI mRNA
    d_mRNAcp = 2.2e-3;       %Degradation rate of amilCP mRNA
    d_LacI = 0;              %Degradation rate of LacI
    d_CP = 0;                %Degradation rate of amilCP
    Fe = 0;                  %Fe2+ concentration (5e6 to 65e6 nM)
    Kd1 = 5.50e4;            %FUR activator dissosiation constant
    Kd2 = 4.6e3;             %pAceB repressor dissosiation constant 
    Kd3 = 8e2;               %pLacI repressor dissociation constant
    
    FURt = 1000;             %FURt concentration
    F = FURt*(Fe^2/(Kd1^2+Fe^2)); %FUR activated concentration
    
    %Stop any values going beloew 0
    f(f<0) = 0; 
    
    %Functions
    mRNALacI = f(1);
    LacI = f(2);
    mRNAamilCP = f(3);
    amilCP = f(4);
    
    %Reactions
    r01 = k_TX1/(1+(F/Kd2)^2.04);
    r02 = d_mRNALacI*mRNALacI;
    r03 = k_TL1*mRNALacI;
    r04 = d_LacI*LacI;
    r05 = k_TX2/(1+(LacI/Kd3)^2);
    r06 = d_mRNAcp*mRNAamilCP;
    r07 = k_TL2*mRNAamilCP;
    r08 = d_CP*amilCP;
     
   %ODEs
    dfdt = zeros(4,1);
    dfdt(1) = r01 - r02;
    dfdt(2) = r03 - r04;
    dfdt(3) = r05 - r06;
    dfdt(4) = r07 - r08;
end
</code>

</figure>

Unfortunately, when we ran the code, we noticed that the model had very little sensitivity to changes in FURt and iron (II) ion concentrations, as can be seen below:

Fe=0, FURt=0

   <img src="" class="left" width="48%"> <img src="" width="48%">

Fe=1000, FURt=0

   <img src="" class="left" width="48%"> <img src="" width="48%">

Fe=0, FURt=1000

   <img src="" class="left" width="48%"> <img src="" width="48%">

Fe=1000, FURt=1000

   <img src="" class="left" width="48%"> <img src="" width="48%">

It appears that there is very little difference in changes in FURt and iron (II) ion concentrations. We tested many ranges and values, but nothing gave any sensitivity unless it was unrealistically extreme (and certainly not in the range of values we could use for FURt and we would test for Fe).

However, it did point something out to us that we had not considered before: the RNA polymerase does not transcribe plasmid 1 then plasmid 2 - it does them both at the same time. This is shown by the amilCP mRNA concentration graph. Because there is no LacI present, there is nothing inhibiting the transcription of mRNA so it is transcribed very quickly. Then, as LacI concentration builds up, it inhibits amilCP mRNA transcription more and more until transcription is equal to degradation (turning point). As transcription inhibition continues to increase, there is a net decrease in amilCP mRNA concentration.

Although our model is testing in vitro values, it appears that the lack of sensitivity is shown in our in vivo lab work too, suggesting it is not necessarily a mistake in the model, but rather a problem with our construct. That being said, we contemplated the idea that there was something wrong with the model, so we checked back over the constants to confirm their values. After clarifying them, we decided to try a different code.

Euler Method

We programmed our own ODE function using the the Euler method instead of using MatLab's preset ODE function to see if that made a difference. MatLab uses an iterative method to find each value of the amilCP after a small time period dt, then repeat up until a specified tmax This is the code we made to represent this model:

            <code style="overflow:visible">%Define constants
d_mRNALacI = 3.04e-3; %Degradation rate of LacI mRNA
d_mRNAcp = 2.2e-3;    %Degradation rate of amilCP mRNA
d_LacI = 0;           %Degradation rate of LacI
d_CP = 0;             %Degradation rate of amilCP
Fe = 2.6e6;           %Initial Fe2+ concentration (5e6 to 65e6 nM)
Kd1 = 5.50e4;         %FUR activator dissosiation constant
Kd2 = 4.6E3;          %pAceB repressor dissosiation constant 
Kd3 = 8e2;            %pLacI repressor dissociation constant

FURt = 1; %FURt concentration
F = FURt*(Fe^2/(Kd1^2+Fe^2)); %FUR activated concentration

%Define times
dt = 10;         %Time period of each iteration
tmax = 144*3600; %when to stop the loop

clf;
hold on;
for Fe=2.6e6:59.2e6:61.8e6
    
mRNALacI = 0;
LacI = 0;
mRNAamilCP = 0;
amilCP = 0;

for t=0:dt:tmax
    if (0 < t)&&(t<900); k_TX1 = 0;    
        else; k_TX1 = 1.33e-3;
    end   
    
    if (0 < t)&&(t<900); k_TX2 = 0;  
        else; k_TX2 = 2.24e-3;
    end   

    if (0 < t)&&(t<1020); k_TL1 = 0;  
        else; k_TL1 = 0.04;
    end        
   
    if (0 < t)&&(t<1020); k_TL2 = 0;    
        else; k_TL2 = 0.06;
    end   
    
    %Define the iterations
    dmRNALacI= k_TX1/(1+(F/Kd2)^2.04) - d_mRNALacI*mRNALacI;
   
    dLacI = k_TL1*mRNALacI - d_LacI*LacI;
    
    dmRNAamilCP = k_TX2/(1+(LacI/Kd3)^2) - d_mRNAcp*mRNAamilCP;
    
    damilCP = k_TL2*mRNAamilCP - d_CP*amilCP;
    
    %Define the species
    mRNALacI = mRNALacI + dmRNALacI*dt;
    LacI = LacI + dLacI*dt;
    mRNAamilCP = mRNAamilCP + dmRNAamilCP*dt;
    amilCP = amilCP + damilCP*dt;
   
   %Plot for iron=2.6e6, then plot iron=61.8e6
   if mod(t,60)==0
       if Fe==0
           plot(t/3600,amilCP,'b.')
           xlabel('Time (s)');
           ylabel('Concentration (nM)');
           title('Concentration of amilCP over time');
           grid on;
           grid minor
       else
           plot(t/3600,amilCP,'r.')
           xlabel('Time (s)');
           ylabel('Concentration (nM)');
           title('Concentration of amilCP over time');
           grid on;
           grid minor
       end
   end
end
end
</code>

</figure>

This time, we decided to plot both the minimum iron (II) ion concentration (2.6e+6) and the maximum on the same graph (61.8e+6) so the difference is more clear. The blue line represents the lower bound, red represents the upper bound. This is the result:

FURt=1, dt=10, t_max=144+3600 (6 days)

   <img src="" class="left" width="100%">

There is no distinguishable difference and the blue line is not even visible. We checked that this was not just an error by pausing the code midway through the run (MatLab plots everything it has covered so far). As you can see, the red line is effectively in the process of covering up the blue line. This is because MatLab plots whatever is chronologically first in the code. It then plots whatever is next and importantly it appears on top of any previous plots. For us, the blue line is the lower bound, red is the upper. MatLab first plots the blue, then the red, so the red appears to cover up the blue.

   <img src="" class="left" width="100%">

We then tried to vary our FURt value. This is the result:

FURt=1000, dt=10, t_max=144*3600 (6 days)

There is an ever so slight increase in final amilCP concentration, but again, there is no distinguishable difference between the two bounds for iron levels. Finally, we put the iron (II) ion concentration values to the extremes to see what would happen. This is the results:

FURt=1000nM,Fe = 0 (blue), Fe = 1e10nM (red), dt=10, t_max=144+3600 (6 days)

   <img src="" class="left" width="100%">

As you can see, there is still no difference, despite the range of values being so extreme. For clarity, we also coded to make MatLab draw the concentration of amilCP against iron concentration after a given time t (4 hours), to see if there was any change we could observe. This is the code we used:

<code style="overflow:visible">%Define initial concentrations
mRNALacI = 0;
LacI = 0;
mRNAamilCP = 0;
amilCP = 0;

%Define constants
d_mRNALacI = 3.04e-3; %Degradation rate of LacI mRNA
d_mRNAcp = 2.2e-3;    %Degradation rate of amilCP mRNA
d_LacI = 0;           %Degradation rate of LacI
d_CP = 0;             %Degradation rate of amilCP
Fe = 2.6e6;           %Initial Fe2+ concentration (5e6 to 65e6 nM)
Kd1 = 5.50e4;         %FUR activator dissosiation constant
Kd2 = 4.6E3;          %pAceB repressor dissosiation constant 
Kd3 = 8e2;            %pLacI repressor dissociation constant

FURt = 1;             %FURt concentration
F = FURt*(Fe^2/(Kd1^2+Fe^2));  %FUR activated concentration

%Define times
dt = 10;         %Time period of each iteration
tmax = 24*3600; %when to stop the loop
    
clf;
hold on;
for Fe=0:0.1e6:61.8e6
    mRNALacI = 0;
    LacI = 0;
    mRNAamilCP = 0;
    amilCP = 0;
    F = FURt*(Fe^2/(Kd1^2+Fe^2));
    for t=0:dt:tmax
    
        
        if (0 < t)&&(t<900); k_TX1 = 0;    
            else; k_TX1 = 1.33e-3;
        end   

        if (0 < t)&&(t<900); k_TX2 = 0;  
            else; k_TX2 = 2.24e-3;
        end   

        if (0 < t)&&(t<1020); k_TL1 = 0;  
            else; k_TL1 = 0.014;
        end        

        if (0 < t)&&(t<1020); k_TL2 = 0;    
            else; k_TL2 = 0.022;
        end   
         %Define the iterations
        dmRNALacI= k_TX1/(1+(F/Kd2)^2.04) - d_mRNALacI*mRNALacI;

        dLacI = k_TL1*mRNALacI - d_LacI*LacI;

        dmRNAamilCP = k_TX2/(1+(LacI/Kd3)^2) - d_mRNAcp*mRNAamilCP;

        damilCP = k_TL2*mRNAamilCP - d_CP*amilCP;

        %Define the species
        mRNALacI = mRNALacI + dmRNALacI*dt;
        LacI = LacI + dLacI*dt;
        mRNAamilCP = mRNAamilCP + dmRNAamilCP*dt;
        amilCP = amilCP + damilCP*dt;

        %Plot the value for iron and amilCP concentration after tmax
       if t==tmax
                plot(Fe,amilCP,'r.')
                 xlabel('Iron Concentration (nM)');
           ylabel('AmilCP Concentration (nM)');
           title('Concentration of amilCP against Iron levels');
           grid on;
           grid minor;
       end
    end
end
</code>

</figure>

These are the results

FURt=0nM

   <img src="" class="left" width="100%">

As you can see, when FURt=0, iron has no effect on amilCP concentrations, which is what we would expect.

FURt =1,000nM

   <img src="" class="left" width="100%">

There is some variation, but it is over a very small region and it is very steep. It is almost like a threshold iron level: below this value, there will be 777nM of amilCP, but above it, there will be ~778.46nM. Additionally, the range of amilCP values is very small, which may prove to be problematic.

FURt=1e9nM

   <img src="" class="left" width="100%">

Interestingly, more FURt actually made the amilCP concentration increase, but increase over a smaller region (steeper increase). This tackles the issue of a very small increase, but worsens the issue of having a small range. We opted to try smaller FURt values to see what would happen. We tried many values to see if anything would work. Although it widened the increasing region slightly, the increase in amilCP is very small. Here is the graph showing what it looked like when FURt=1nM

FURt=1nM

   <img src="" class="" width="100%">

So it seems as though increasing FURt increases maximum amilCP, but shortens the region of sensitivity, while decreasing FURt decreases maximum amilCP and widens the region of sensitivity. It doesn't look as though it would be possible to find a balance between the two. Even when the change in amilCP is just ~0.005nM, the region still doesn't come close to covering our range of values for iron.

We wanted to zoom in on the increasing region, so instead of testing up to 61.8e6nM, we kept testing smaller and smaller regions until we isolated just the increasing section. This graph represents iron levels varying from 0 to 0.6e6 (note that our lowest saliva iron level is 2.6e6)

FURt=1000nM

   <img src="" class="left" width="100%">

Unfortunately, we were unable to find a source showing an image of amilCP at a set concentration, and we didn't have the equipment to find this ourselves. This would give us an indication of when it became visible (eg,777nM could be invisible or very intense). We took this into account and tested to see if there would be any steady increase in amilCP levels on smaller timescales. We used t=1050 (30 seconds after translation begins) and t=1200 (arbitrary increase). These are the results:

time = 1050s

   <img src="" class="left" width="100%">

time = 1200s

   <img src="" class="" width="100%">

As you can see, there is still the notable threshold iron level, below which amilCP is lower, above which amilCP is higher. There is no steady increase over the range of values we would need to test. We simulated for many different times, and even varied FURt to see if that made a difference, but the pattern continued.

Finally, we simulated a larger time scale of 24 hours and varied FURt again, but the same pattern emerged:

FURt=1000nM

<img src="

" width="100%">

Finally, we tried varying DNA concentration. We did this by quoting a different value for k_tx from our literature. All previous values had been based on a DNA concentration of 60nM, but we can also try simulating 10nM, but changing k_tx. This is the graph it produced:

FUR=1000nM, T=4 hours, DNA=10nM

<img src="

" width="100%">

As you can see, the same issues arise, just at a lower amilCP range.

Analysis

Summary

We set out to prove whether our simple iron test would work in an in vitro environment in a cell-free system. We first successfully modelled our reactions in terms of ODEs. We then attempted to solve the equations in with three different methods: analytically by integration, using the MatLab ODE function and using the Euler method. Unfortunately, it was not possible to solve analytically. but both the ODE function and Euler method worked and gave similar results. The most conclusive simulations being the iron vs amilCP graphs, as they show exactly how big the difference is in amilCP and iron varies, and with different FURt and time values.

Conclusion

Our cell-free system is only sensitive over a very small region of iron levels (0 to 0.6e6nM) which is below our minimum iron level of 2.6e6nM. Although decreasing FURt made it sensitive over a larger region, it still didn't cover any significant amount of our ranges and the change in amilCP was very minute when FURt was small. Likewise, increasing FURt made amilCP increase by more, but made the test sensitive over a smaller region. The model therefore quite conclusively states that the simple iron test we originally intended to make will not work. Although our wet-lab work was in vivo, it should be noted that a similar issue was observed - the test had no sensitivity to iron levels.

However, this certainly doesn't mean the end of our project. We believe this global issue is very important, so a solution would be incredible useful. This being the case, we have a few things we could adjust to our project, based on the model. Perhaps our most promising was diluting the saliva sample. If we diluted the sample in a ratio of 1:99 with distilled water (which we could provide with the test kit), then the concentration of iron would be 100 times smaller, meaning it would fit almost exactly onto our 'zoomed in' region. Suddenly, our system is sensitive to the iron levels we need it to be. This being the case, our goals of finding appropriate values for time, FURt and DNA become feasible. As an example, here is a graph of varying FURt values.

FUR= 1nM (blue), 1000nM (red), 1M (green); T=4 hours; DNA=60nM

However, we still have the big issue of whether we can see a visible change in colour over such a short amilCP range. We ideally wanted zero iron concentration to correlate to zero amilCP concentration (y-intercept), but it doesn't. We would need a plate reader and a means to find amilCP concentration experimentally, but we fear there may not be a visible change in amilCP no matter what we do, but we didn't have the time nor equipment to explore this avenue in more detail.

Alternatively, we could vary the system itself. We hypothesised that our promoters were quite leaky, which was causing a lot of sensitivity issues. Unfortunately, we could only find 1 iron sensitive promoter in the iGEM registry, so we didn't have much choice. Perhaps once a new one is characterised (or an alternative system entirely is created by a future team), our project could be revitalised.

Evaluation

As with any good model, we must evaluate the limitations of our model, highlighting any factors we have not accounted for and any errors that may have occurred, as well as how they could affect our conclusion.

We assumed protein degradation rate was zero: Protein degradation is very slow [23, 24], so we felt it was valid to assume it was zero, but nonetheless we could include it.

Did not account for reaction conditions: The reaction conditions of the system were ignored- temperature, pressure and pH.

Other limitations on growth were not taken into account: Although we considered degradation, there are also other factors that restrict the velocity of our reaction. One example would be space. We essentially assumed there was an excess of space. However, this may not be the case, especially as we are dealing with concentrations on a μM scale. Another would be the concentrations of the enzymes were not considered.

Did not account for diffusion rates: We assumed that everything was perfectly and evenly spread among the reactants and ignored the idea that it may take time for the various parts to reach each other before they react. However, because we were considering very long periods of time and a very small surface area, this is a valid assumption and would have little effect on the outcome of our model.

Possible overuse of sources: We used Eyal Karzbrun et al [7] quite a lot in our model. This was primarily because it was a study of protein synthesis in a cell-free system, which is exactly what our project is. However, we also considered other sources for each of our values and took these into account when we decided on our final values for our constants.

Did not check extensively the reliability of sources: Little time was spent evaluating the reliability of sources. Although we checked that the values we found were the values we were trying to find, we did not necessarily check whether the source was particularly reliable.

Only 1 member of the team was working on the model: By only having one person working on the project, we were unable to have a second opinion at every step. In this way, mistakes are significantly more likely to occur and having multiple team members working on the model would make it much more reliable.

We started quite late, meaning there was a lot of time pressure: Because our team had no experience in modelling, it took the team a while to get started on the model. All the work was therefore done in the course of about one month, by a highly inexperienced student. Although this could not have been avoided, it should still be stated.

With these limitations on our model, it would be improper to suggest our model is 100% accurate and that our conclusion is correct. However, we feel as though we have tried to account for as much as possible, and the parts we did not account for would not have a big effect on the model. The fact that our lab work supports this conclusion makes us even more confident that the model is accurate.

We wish to thank George Smith and Adam Jones from the City of London School for Boys iGEM team and Chun Ngai Au from the Oxford iGEM team for guidance, assistance and discussion for our model.

References

[1] Chandran, D. et al. “Mathematical Modeling and Synthetic Biology.” Drug discovery today. Disease models 5.4 (2008): 299–309

[2] Hamed, Mazen Y., J. B. Neilands, and V. Huynh. "Binding of the ferric uptake regulation repressor protein (Fur) to Mn (II), Fe (II), Co (II), and Cu (II) ions as co-repressors: Electronic absorption, equilibrium, and 57 Fe Mössbauer studies." Journal of inorganic biochemistry 50.3 (1993): 193-210

[3] Anne Bagg and J. B. Neilands “Molecular Mechanism of Regulation of Siderophore-Mediated Iron Assimilation.” Microbiological Reviews. Dec. 1987, p. 509-518.

[4] Becskei A and Serrano L "Genetic network driven control of PHBV copolymer composition", J Biotechnol 122(1):99-121, 2006

[5] Deng, Zengqin, et al. "Mechanistic insights into metal ion activation and operator recognition by the ferric uptake regulator." Nature communications 6 (2015).

[6] Basu S et al "A synthetic multicellular system for programmed pattern formation", Nature 434:1130-1134, 2005

[7] Eyal Karzbrun et al “Coarse-Grained Dynamics of Protein Synthesis in a Cell-Free System” Physical Review Letters, Jan 28 2011, PACS numbers: 87.18.Vf

[8] Semsey et al “The effect of LacI autoregulation on the performance of the lactose utilization system in Escherichia coli” 2013 Jul; 41(13): 6381–6390

[9]Selinger, D. W. "Global RNA Half-Life Analysis In Escherichia Coli Reveals Positional Patterns Of Transcript Degradation". Genome Research 13.2 (2003): 216-223. Web. 18 Oct. 2016.

[10] H. Bremer and P. P. Dennis, in Escherichia Coli and Salmonella: Cellular and Molecular Biology, edited by F. C. Neidhardt et al. (ASM Press, Washington, DC, 1996), pp. 1553–1569, 2nd ed.

[11]https://www.abbexa.com/ecoli-ferric-uptake-regulator-protein-recombinantgclid=Cj0KCQjw95vPBRDVARIsAKvPd3J1bjFD-mag6B0xR1MofivOa6KgOigQoMOsmw8TbRZZDVBj-_C8WjAaApo5EALw_wcB

[12] Setty Y et al "Detailed map of a cis-regulatory input function", P Natl Acad Sci USA 100(13):7702-7707, 2003

[13] http://www.nature.com/nature/journal/v274/n5673/abs/274765a0.html

[14] http://parts.igem.org/Part:BBa_K592009

[15] http://www.uniprot.org/uniprot/P03023

[16] http://parts.igem.org/Part:BBa_K1163003

[18] https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3711431/#gkt351-B31

[19] https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3435124/table/t1-mjhid-4-1-e2012051/

[20] O.L. Miller, Jr., B. A. Hamkalo, and C. A. Thomas, Jr. “Visualization of Bacterial Genes in Action” Science 169, 392, (1970)

[21] J. Shin abd V. Noireaux, J. Biol. “Study of messenger RNA inactivation and protein degradation in an Escherichia coli cell-free expression system” Eng. 4, 9 (2010)

[22] J. A. Megerle, G. Fritz, U. Gerland, K. Jung, and J. O. Ra¨dler, Biophys. J. 95, 2103 (2008)

[23] https://slack-redir.net/link?url=http%3A%2F%2Fpubs.acs.org%2Fdoi%2Fabs%2F10.1021%2Fja954077c

[24] http://pubs.acs.org/doi/abs/10.1021/ja00230a041

[25] Ali Tiss et al "Characterization of the DNA-binding site in the ferric uptake regulator protein from Escherichia coli by UV crosslinking and mass spectrometry" FEBS Letters 579 (2005) 5454–5460

[26] Pardee et al., Paper-Based Synthetic Gene Networks, Cell (2014), http://dx.doi.org/10.1016/ j.cell.2014.10.004

Team:Judd UK/Pages/Model

Contents

Model

Abstract

Model formation

Defining our constants

Model Constant Table

Solving the equations

ODE method

Euler Method

Analysis

Summary

Conclusion

References