Team:Valencia UPV/Modeling/Simulations

Valencia UPV iGEM 2015

Do you want to download the Matlab files to test and see our mathematical model? Here you have the data file and the executable file.

After developing the mathematical model and thinking about the equations we developed our model in Matlab, which can be used to perform other experiments with different values or light pulses. In Matlab we performed different experiments in order to see how our biological decoder would react to different light pulses and to try different ideas.


Constitutive level: Zero input response

Before studying saturation when light pulses were given, we must take in advice that first level (producing A, B, C and D), starts expressing from minute 0. Activation of next level is directly proportional to the amount of those proteins that were produced firstly, as they are binding domains of secondly activated constructions. Thus, more expression will be obtained as more binding domains are synthetized.

We determine the moment when constitutive expression reaches its plateau by simulating the zero input response of our system during 10000 minutes.

Light Off sequence:

Figure 1.

Constitutive protein production:

Figure 2.

Protein production behaves like a first order system after a step input. In this case, the step from 0 to 1 is the beginning of the protein production inside the cell. Thus, the moment when it will saturate is approx. 5·τ, where τ is the time constant of this process, given by the moment at which production has reached the 63% of its saturation level.

Figure 3.

Searching again with the data cursor:

Figure 4.

Now we can fix the beginning of light irradiation in 250 min, and the initial concentrations of proteins BD1-PIF6, PhyB-VP16, BD2-LOV2 and ePDZ-VP16 in 250 nM. This is the reason why all our inputs will have the structure:

input_red =[zeros(251,1);… ;zeros(100,1)];

input_blue=[zeros(251,1);… ;zeros(100,1)];

The final vector of zeros is added in order to appreciate the effect of degradation when lights are switched off. This could give a qualitative idea of how much of the production can be lost during the extraction of the product for its ingestion.

Switches level.

We opted to make the light pulse “infinitely” long in order to notice until when does production rise and its saturation point. Firstly, we did only a red pulse of 10000 minutes, so the final product should be alpha.

t_final=10000; % with input_t=(0:Ts:t_final)' (minutes)

input_red =[zeros(251,1);red_intensity*ones(length (input_t)-351,1) ;zeros(100,1)];

input_blue=[zeros(251,1); zeros(length (input_t)-351,1) ;zeros(100,1)];

Following graphics let us determine the moment when second and last levels of production, reach their plateau

Figure 5.

Figure 6.

Inputs are on the top of the graphics. Second line belongs to second level expression and green curves represent the outputs performance. Using the data cursor in both elements whose production was being induced (protein E and alpha):

Searching again with the data cursor:

Figure 7.

With inputs of 250 minutes, the plateau will be reached. If carry on with the irradiation after 250 minutes or not, relies on the efficiency of recombinases.

Figure 8.

With inputs of 2510 minutes, maximum output production can be reached. Then:

This will be the duration considered in next simulations.

Recombinases action

Inputs used in our circuit are determined by the value of two variables: color and time. Distinguish different colors is easy due to their different wavelength, resulting in the stimulation of different photoreceptors that regulate each switch. But how can our genetic constructions assign an order to inputs of the same colors?

Recombinases are the biological elements which grant memory to AladDNA. These enzymes work as excisionases and integrases, cutting the gen which is surrounded by recognition sites specific for those recombinases. In our model this excision has been represented by the change from g+ to g-, referring to the capability of the gen for being translated.

In our device, recombinases used are BxB1 and θC31. They are produced with the first input:

Both second levels have the recognition sites for recombinases produced in the analogous level. Thus, enzymes firstly produced recombinate the option which has not been activated. When the second input is introduced, production in the second level will not be possible, as genetic codes have been flipped. In other words, if the second light is a different color, it will only mean production regulated by switches of the final level, resulting in a circuit that reminds which was the first input.

How do they behave?

As we explain in the biochemical reactions section, when recombinases are produced they instantaneously dimerize, and two dimmers will be needed to recombinate each gene. The duration of this process will determine the moment when it will be save to give the second input.

Firstly, we will plot the results of a single red input in order to appreciate the curve described by g+ in the analogous level. Gene copy number used for this performance is of a single copy.

t_final=3100; % with input_t=(0:Ts:t_final)' (minutes)

input_red =[zeros(251,1);red_intensity*ones(length (input_t)-351,1) ;zeros(100,1)];

input_blue=[zeros(251,1); zeros(length (input_t)-351,1) ;zeros(100,1)];

Figure 9.

Graphic on the left down corner represents the amount of BxB1 produced. In comparison to the graphic from Picture 10, this gradient is minor, which is due to the usage of monomers produced to form dimmers as soon as proteins are produced. Looking closer to the second line of results, we can determine from the left graphic the moment when G will reach its maximum (and so will G2).

Figure 10.

Figure 11.

Figure 12.

On the other hand, the moment when recombinases switch off gH+, gI+ and gJ+, is approx. the minute 280, as we can see in the graphic below:

Figure 14.

In order to prove recombinases efficiency, we will simulate two followed inputs of different colors. If they work, the expected result is that there is no expression of that genes later, when the other light is given.

t_final=3100; % with input_t=(0:Ts:t_final)' (minutes)

t_recomb=280; %time needed for recombinases work

input_red =[zeros(251,1); red_intensity*ones(t_recomb),1); zeros(length(intput_t)-t_recomb-251,1)];

input_blue=[zeros(251+length(t_recomb),1); blue_intensity*ones(length(input_t)-t_recomb-351,1); zeros(100,1)];

Figure 15.

As we guessed, no expression of the second level is produced after a previous input of 280 minutes, which is the time that ensure us the end of recombinases action.

Influence of gene copy number

Different vectors can be used to introduce our genetic constructions in the cell. The amount of gene’s copies varies in each kind of plasmid, and this is a critical parameter. There will be a trade-off between the amount of genes expressing enzymes, and the number of genes that those will have to recombinate.

We performed simulations with different values of gene copy numbers relying on some vectors typically used to transform E.coli.

Values considered for this parameter were:

gin=[1,5,10,17,400,600]; %Initial amount of genes

Figure 16.

From 1 to 17 copies:

Figure 17.

With more genes, there is more basal expression of recombinases, but also of those genes that BxB1 and θc31 must recombinate. Moreover, levels of basal expression of the analogous recombinases will also increase. This is the reson of the decay in gE+, which is more abrupted as the gene copy numer rises. That basal recombination is undesired, because the amount of binding domain produced also decays, as we can appreciate in the second group of curves (E). However, with 5 copies the difference is positive, but this arise does not persist when the number of copies keeps on increasing.

Figure 18.

Graphic above suggests that recombinases efficiencyis not affected by the numer of copies, because they can get the same result (g+  0) in the same time. Then, if increasing the number of copies means less amount of expression and it makes no difference in recombinases efficiency, these results suggest that low numbers of copies are more suitable to get the desired behavior.

Results with 400 and 600 copies support that higher amounts of copies induce more basal expression of recombinases, which leads in an undesired behaviour of our circuit:

Figure 19.

No protein production is appreciated. Looking closer to the first 100 minutes of simulation:

Figure 20.

Basal expression of both recombinases is so high that even when light is still off, they are produced in enough quantities to totally switch off production in the second level.

Light stimuli

Pause between pulses

We knew from literature, that changes induced by light are instantaneous. However, there are several differences between switching on and off these switches. For example, if we had implemented a red toggle switch, we would be able to ensure the stop of production using far red light, whereas dark periods induce a slower decay.

This means that although light irradiation is stopped, the presence of photosensible elements may will decrease slowly, and production will not occur only when light is given. In other words: production will not be totally over light control.

In order to ensure this coordination between light and presence of active elements, we suggested our wet lab mates tagging PhyB and ePDZ, they may will be degraded fast enough to minimize their presence when light has been eliminated. This is the reason why B and D have higher degradation rates than other proteins.

So, if dark periods mean loss of activators, pauses between first and second inputs should be as short as they can be, in order to maximize production. The shorter pause possible, is a no-pause, i.e continuous irradiation.

Firstly, we will compare production of alpha with different pauses between both red inputs, in order to analyze the influence of the dark period length. We implemented in Matlab a for loop to facilitate the comparison:

pause_dur=[0,2,30,60,80,100]; % vector containing different duration of pauses.

for i=1:length(pause_dur)

light_dur=length(input_t)-351-pause_dur(i);% duration of irradiation





[t,x] = ode23s(@model_switch_v2, time_span, ic, opt,input_t,input_red,input_blue,param);

figure (1)

subplot(221);plot(t,[x(:,19)],'g');xlabel('t');ylabel('alpha[···]');grid on;hold on

Here we plot the results. All time axis are expressed in minutes.

Figure 21.

Fixing our attention in alpha:

Figure 22.

It is clear that production decays with the increase of the pause between inputs. In other words, production is increased with more irradiation. Then, combinations of the same color repeated (red_red and blue_blue), will provide better results if they are one single pulse instead of two.

Inputs combination: Sequences

From conclusions extracted above, we will assume that both pulses must be continuously given, with no pause between them. Furthermore, in combinations red_blue and blue_red, there will be a necessary condition that ensures waiting time enough in order to let recombinases finish their work. The duration of this period was estimated before, concluding that 280 minutes of simulation were enough in order to give the second pulse. In summary:

Another criteria which may be useful, is the moment when production reaches its plateau. This leads in abundance of binding domains required by next level. Output amounts are directly proportional to the presence of those binding domains. From saturation experiments we determined that since production begins, it reached its maximum 250 minutes later. As second level expression arises after an initial pause of 250 minutes, the absolute time will be:

Our first idea was comparing three kinds of structure:

- Both inputs with the same length.

- First input longer than the second one.

- Second input longer than the first one.

Inputs introduced in Matlab in order to test each possibility, and results obtained were:

input_red =[zeros(251,1);red_intensity*ones(2750/2,1) ;zeros(length(input_t)-(2750/2)-251,1)];

input_blue=[zeros(251+(2750/2),1); blue_intensity*ones(2750/2,1) ;zeros(100,1)];

Figure 26.

Figure 27.

Figure 28.

Figure 29.

input_red =[zeros(251,1);red_intensity*ones((2750/2)+500,1) ;zeros(length(input_t)-(2750/2)-500-251,1)];

input_blue=[zeros(251+(2750/2)+500,1); blue_intensity*ones((2750/2)-500,1) ;zeros(100,1)];

Figure 30.

Figure 31.

Figure 32.

Figure 33.

input_red =[zeros(251,1);red_intensity*ones((2750/2)-500,1) ;zeros(length(input_t)-(2750/2)+500-251,1)];

input_blue=[zeros(251+(2750/2)-500,1); blue_intensity*ones((2750/2)+500,1) ;zeros(100,1)];

Figure 34.

Figure 35.

Figure 36.

Figure 37.

In simulations performed above, the expected output was β, and the third scenario is that one which gets higher levels of its expression. Then, the second input must be longer that the first one. If we optimize this condition, our aim in following simulations will be to find the first and second inputs as short and as long as they can be, respectively.

Since necessary condition is achieved in (tmax.BD > 280 min), theoretically, the best option would be giving the second input in minute 500. However, we will simulate different scenarios in order to test which one gets the best result. As the order of inputs is red_blue, the expected output is β.

input_red =[zeros(251,1);red_intensity*ones(281-251,1) ;zeros(length(input_t)-281,1)];

input_blue=[zeros(281,1); blue_intensity*ones(length(input_t)-281-100,1) ;zeros(100,1)];

Figure 38.

In graphics below prove recombinases condition is achieved, but not max BD production:

Figure 39.

Figure 40.

Outputs amounts obtained:

Figure 41.

input_red =[zeros(251,1);red_intensity*ones(250,1) ;zeros(length(input_t)-501,1)];

input_blue=[zeros(501,1); blue_intensity*ones(length(input_t)-501-100,1) ;zeros(100,1)];

Figure 42.

Maximum BD production is achieved, and also is action of recombinases:

Figure 43.

Figure 44.

Outputs produced:

Figure 45.

If we give time enough for recombinases to actuate, giving the second input at minute 280, but we do not stop production of the BD until it reaches its maximum (t1st input=500 minutes), there will be an overlap between pulses, but both conditions are accomplished as well.

input_red =[zeros(251,1);red_intensity*ones(250,1) ;zeros(length(input_t)-501-100,1)];

input_blue=[zeros(281,1);blue_intensity*ones(length(input_t)-281-100,1) ;zeros(100,1)];

Figure 46.

Figure 47.

Figure 48.

Comparing the production of all products in each of the three last simulations:

Figure 49.

Figure 50.

Figure 51.

Figure 52.

The third option presents the worst results of all. As irradiation with both lights is more prolonged, more undesired outputs are produced.

Between the first and the second option, the main difference is between the production of alpha and gamma in the first case, and the production of Omega in the second one. At the end, ratios achieved in both cases:

Figure 53.

Figure 53.

Ratio beta versus alpha is higher if the second input begins just after recombinases work, i.e minute 280.

Although the logical structure of inputs would be one after another, we found interesting the idea of simulate some unusual conditions such as two inputs partial or totally simultaneous, or sequences with three structure color1_color2_color1.

Three light inputs

input_red=[zeros(276,1);red_intensity*ones(900,1) ;zeros(900,1);red_intensity*ones(900,1);zeros(125,1)];

input_blue=[zeros(1176,1);blue_intensity*ones(900,1) ;zeros(1025,1)];

Figure 54.

Figure 55.

input_red=[zeros(1176,1);red_intensity*ones(900,1) ;zeros(1025,1)];

input_blue=[zeros(276,1);blue_intensity*ones(900,1) ;zeros(900,1);blue_intensity*ones(900,1);zeros(125,1)];

Figure 56.

Figure 57.

Three pulses are not valid to obtain just one output and consequently, this combination should be avoided if only one output is wanted, because the presence of all which is produced and its possible intake, may have negative consequences.

However, if multiple production is desired, this kind of sequence permits the production of two outputs in different proportions, which could be useful in situations such as obtention of products which need a “helper” compound.

Simultaneous inputs

input_red =[zeros(251,1);red_intensity*ones(length(input_t)-351,1);zeros(100,1)];

input_blue=[zeros(251,1);blue_intensity*ones(length(input_t)-351,1) ;zeros(100,1)];

Figure 58.

Figure 59.

input_red =[zeros(251,1);red_intensity*ones(2750,1); zeros(100,1)];

input_blue=[zeros(251,1);blue_intensity*ones(1375,1) ;zeros(1475,1)];

Figure 60.

Figure 61.

input_red =[zeros(251,1);red_intensity*ones(2750,1); zeros(100,1)];

input_blue=[zeros(1626,1);blue_intensity*ones(1375,1) ;zeros(100,1)];

Figure 62.

Figure 63.

input_red =[zeros(251,1);red_intensity*ones(1375,1); zeros(1475,1)];

input_blue=[zeros(251,1);blue_intensity*ones(2750,1); zeros(100,1)];

Figure 64.

Figure 65.

input_red =[zeros(1626,1);red_intensity*ones(1375,1); zeros(100,1)];

input_blue=[zeros(251,1);blue_intensity*ones(2750,1); zeros(100,1)];

Figure 66.

Figure 67.

With this sequences we can produce majorly 2 outputs, or three with one of those three majorly produced. This sequence offers the possibility of multiple production, combining the effect of different compounds, which could be useful in scenarios as the production of multivitamin complex.

Inputs choice for each output

Using information recovered from simulations described, we were able to determine light sequences that obtained ratios of nearly 100, being the desired output the one majorly produced with a great difference. In summary, these are the relations between inputs and outputs that describe AladDNA’s performance:

Output alpha:

input_red =[zeros(251,1);red_intensity*ones(length(input_t)-351,1) ;zeros(100,1)];


Figure 68.

Figure 69.

Output beta:

input_red =[zeros(251,1);red_intensity*ones(281-251,1) ;zeros(length(input_t)-281,1)];

input_blue=[zeros(281,1); blue_intensity*ones(length(input_t)-281-100,1) ;zeros(100,1)];

Figure 70.

Figure 71.

Output gamma:

input_red =[zeros(281,1) ;red_intensity*ones(length(input_t)-281-100,1);zeros(100,1)];

input_blue=[zeros(251,1); blue_intensity*ones(281-251,1) ;zeros(length(input_t)-281,1)];

Figure 72.

Figure 73.

Output Omega:

input_red =[zeros(length(input_t),1)];

input_blue=[zeros(251,1);blue_intensity*ones(length(input_t)-351,1) ;zeros(100,1)];

Figure 74.

Figure 75.