Team:ULB-Brussels/Modeling

From 2010.igem.org

Revision as of 11:26, 27 October 2010

•   Home  
•   Team  
•   Project  
•   Parts  
•   Modeling  
•   Notebook  
•   Safety  
•   Sponsors  

Questions from the Wetlab part

Questions asked by the dihydrogen module were mainly general questions : how much the deletion of the whole, or of one part, of the genes proposed would increase the amount of dihydrogen produced ? To answer these questions, the metabolic pathway of mixed acid fermentation had to be modelled before we can, eliminating equations for some of the arms of the pathway, observe the advantages in term of dihydrogen production.

Initial equations and further enhancements

We posited the differential equations expressing the variation in the main intermediates of the metabolic pathway, vs the time. Generally speaking, these equations comprise a creation term (positive) representing the formation of this intermediate from a reactive which is before it in the pathway and from one or several destruction terms (negative) that stand for the degradation of this intermediate in the next product. We distinguish two types of terms (creation and destruction ones). First of all, they are linear terms, characteristic of non-catalyzed reactions that we assume being spontaneous. They are of the type coefficient.reagent concentration. There is, then, the big category of the non-linear terms, for more complex situations, in which enzymes act. In the H2 module, these non-linear terms are Michaëlis-Menten like, which is the typical cinetic of enzyme-catalyzed reactions. We'll then see, in the next section, other non-linear terms, more particulary used when we modelize systems with important genetic regulations. Here are, then, the first equations which represent the metabolic pathway of the mixed acid fermentation.

v0 is here a constant of production of the phosphoenolpyruvate by glycolysis. The two next terms are enzymatic terms of first order, characterizing the catalyzed reaction ADP + PEP <=> ATP + PYR + 2H+ (v1) and the catalyzed reaction PEP + CO2 + H2O <=> OXA + P + H+ (v2)

The first term is the term of production of the pyruvate, from phosphoenolpyruvate The second term (in v3) is a term of Michaëlis-Menten for the catalyzed reaction PYR+ CoA  FOR + AcCO. The third term (in v4) is once again a non-linear term traducing the catalyzation of pyruvate disappearance, which is transformed in lactate, following the reaction NADH + PYR + H+ <=> LAC + NAD+ )

The first term is a term of formation of the formate, from pyruvate. The second term (in v5) traduce the catalysation of the reaction FOR + H+ <=> CO2 + H2. The third term (in v6) represents facilitated diffusion of formate, by specific and non-specific permeases. It's again an enzymatic reaction, with a typical non-linear kinetic.

The first term is a term of formation of the acetyl-CoA, from pyruvate. The second (in v7) is a non-linear term of the second order in OXA and AcCO, traducing the reaction OXA + AcCO+ H2O <=> CIT + CoA + H+

The first term is a term of apparition of the oxaloacetate, from PEP The second term is a non-linear term of the second order in OXA and AcCO, traducing the reaction OXA + AcCO+ H2O <=> CIT + CoA + H+

The first term is a term of apparition of the citrate from acetyl-CoA and oxaloacétate The second term is a term of consumption of the citrate, in following reaction.

The first term is the term of apparition of the lactate from pyruvate. The second term (in v9) is the facilitate diffusion term for lactate.

We then realized the importance of an equilibrium term in the last reaction leading to the formation of hydrogen  (FOR + H+  <=> CO2 + H2.),this product can accumulate, which encourages the reverse reaction.
We also refined the kinetic of the reaction OXA + AcCO+ H2O <=> CIT + CoA + H+, using a non-linear term which considers the concentration and Michaëlis-Menten constant specific of each substrate.
Finally, and to more easily observe the benefit of each deletion, we added an equation modeling the output speed of H2 diffusing out of the bacteria.    
Here are the final equations :

Where the term in v5r is the equilibrium term characterizing the reverse reaction of the H2 production. This term has been added to account for the balance of the reaction, and avoid an irrealist of H2 in the cell.

Where the term in v7 is a modification of the preceding term in v7, modeling a more general kinetic, to better take into account the differences between the two substrates. Where the term in K7b is a term of utilisation of acetyl-CoA in later reactions.

With the same modification than above.

With the same modification than above.

Where the term in v5r is the equilibrium term characterizing the reverse reaction of the H2 production. This term has been added to account for the balance of the reaction, and avoid an irrealist of H2 in the cell.

Results

Such a differential equation system can't be analytically solved. We then had to use a specialized software able to numerically integrate the system. To do that, we used the XPP software, specifically designed to model biological systems. This software allows us to study the evolution of each variable depending on time. It also has an extension, Auto, able to draw stability diagrams of the systems. At first, we only vary the parameters, depending on the expected wetlab results. For example, v2=0 means the PPC gene has been Knock Out successfully (see table 1 for the corresponding between genes deletions and the parameters.). We can see, in the diagram 1, the theoretical effects of such a deletion on the intra-cellular concentration of dihydrogen at the steady-state: it's almost doubled (from 0.057 to 0.118). In the same way, we put each parameters to a nul value, alternately and in a combined way, to observe the effects of the deletions and their optimal combination.

We noticed that for double deletions, the results cluster in three increased levels of production (cf. Diagram 2). A deletion with particularly strong associated effect: FocA. LdhA and PPC are middle effect deletions, whereas the deletion of both hydrogenases seems to have a lesser effect on productivity. Indeed, the combinations including FocA are always in the two higher levels of production, while those including the hydrogenases are always in the two lowest levels.

In triple deletions, the above results are confirmed and show even more the importance of FocA in dihydrogen production, while the weakness of the hydrogenases deletion seems less important. Indeed, we notice here two very marked dihydrogen production levels, one being FocA and two other genes (no matter what ever they are) KO bacteria and the other one with a functional FocA.  

There is an important clarification to make here. The intracellular concentration in dihydrogen presented here are in no case absolute values, they are given in arbitrary units. Generally speaking, all of the numbers presented in this section are only for a comparative purpose. This is easily understood if we remember the largely empirical way to establish the parameters of the differential equations.

Questions from the Wetlab part

Here, after a general modelisation to theoretically validate the system chosen, we had to answer more specific questions, and with an immediate application : it was to choose the relative and absolute strength for the Ribosome Binding Sites (RBS) upstream of the genes used in this module, and to model the possible influence of the sequence the genes are built into. To obtain a functioning module, we wanted the system created by the interaction of both plasmids not to compromise the long-term population density; and that after a dilution of a ten to a hundred times, the density of bacteria would immediately fall and no new bacteria growth would be observed in the river

Initial equations and further enhancements.

When we have, like in this model, a relative complexity in the genetic regulations and we do modelize them, we generally use in addition to Michaëlis-Menten terms, to Hill terms.Hill terms are non-linear, directly derived from Michaëlis-Menten like laws. They are proportional to C^n/C^n+K^n with n>=1 and where C is the concentration of the reactive and K in a constant.
To determine which order of the Hill terms to choose, we used Hill term order=number of effective regulations levels. This is experimentally because we see, at each regulation level, proteins characterized by a Michaëlis-Menten type kinetic. If there are two regulations levels, we can reasonably assume that these saturations effects combine.
In a simple model, like the one of the dihydrogen, we can admit that approximate amount of transcripted genes= the amount of produced protein, with a lack of regulations systems modeled. This order is thus equal to one.
It's interesting to note that the order of the Hill's term is representative of the complexity degree of the dynamic of the system, and that terms which an order equal to two or more often have complex stability diagrams, including bistability phenomenons.
This model is richer in regulations level, since we here consider the detail of the genetic expression and his regulations. We cannot posit anymore the approximation that the amount of ARN is constant, and to take into account this variability of the transcription, and the complexity of this ARN stage, with the cooperative effects of the ARN-polymerase complex, we attribute the second order to the Hill terms put for the creation and consumption steps of the proteins which information are on the genes found in the plasmids we are talking about. Indeed, we have two proteins involved in protein synthesis. The first equation is more specific, since it models the bacterial density. This density classically increase, with an exponential phase followed by saturation of the system, hence the non-linear term of the logistic growth. The population, in addition to that density-limited growth, is also killed by the poison, proportionally to the bacterial density.

The first term is the natural growth tendency of a bacterial population.
The second term is the cell death term due to the poison.

The first term express the fact that the signal-molecule HSL (cf. the Cell death module in Wetlab part, for a more detailed introduction to the role of signal molecules and a characterization of HSL) induce his own expression (it's here a typical induction, which obeys, too, to the saturability of substrates law), and that this is modulated by the bacterial density: the more bacteria we have, the more probable will be the entry of HSL in the cell. The second term is a degradation term of HSL by the specific enzyme; the third one is a non-specific degradation or diffusion term

The first term express the inhibition by HSL of the promoter which controls the gene coding for the HSLase.
The second term is a non-specific degradation term.

The first term is identical to the one of the above equation, although with a different coefficient, because of the possibility of a different strength RBS. The second term expresses the mutual neutralization of the antidote and poison, by complex formation. The third term is a non-specific degradation term.

The first term express the activation by HSL of the promoter, just as in HSL equation but with a different coefficient due to RBS.

In a second time, to refine the system, we also chose to better take into account the subtleties of the culture medium, in which we add HSL in growth conditions, by clearly separating the two co-activators (or inhibitors) of the promoter. It was also to add equations, permitting to better take into account the formation of complexes, important in this system, since these complexes activate or inhibit the whole system. We again explicitly made the difference between the concentration in ARN and in protein.
Finally, we have distinguished the intra and extra-cellular HSL, to better take into account the HSL massive outflow phenomenon, by diffusion, during the dilution. This will allow our model prediction permitted by our model to predict a possible peak of HSL , causing a certain stability (and even a slight initial growth) in the population, immediately after the dilution.
Here are the final equations :

LUXR is at the same time the transcript of the gene and the corresponding protein: the coactivator (or inhibitor) with HSL of the promoters. R2 is a coefficient put to model the relative strength of the RBS. LUXR is under control of the promoter activated by the LUXR-HSL complex, hence the creation term, and we notice the second order of the Hill term. The second term is a non-specific degradation term and the last one comes from the dissociation of the complex.

LUXI designates the gene coding for HSL. We notice again the coefficient of the RBS and the activation by the complex, together with the non-specific degradation term.

HSL is translated using the LUXI transcript, hence the first term. The second one is a degradation term by the HSLase, proportional to this enzyme and to HSL. The third term is a loss term due to diffusion outside of the cell and the last one symbolizes the reverse phenomenon: the entry of extra-cellular HSL. Fke is a function of N, which allows us to model the fact that, the more diluted the population is, the more HSL will leave the cells, by passive diffusion. This will explain the massive outflow of HSL observed at the very first stages following an important dilution.

HSLe (extra-cellular) is modulated by the diffusion in and out of bacteria. The last term is a loss term: when HSLe is too far away, it has no influence on the system anymore.

LUXRHSL represents the complex formed by these two signal molecules. Its concentration is proportional to the one of its constitutive parts, minus the quantity of complex that dissociates and a non-specific loss term.

Once again, we can see a coefficient put for the RBS. Here, the complex formed by the signal molecules inhibits the activity of the HSLase, which is under control of the poison-promoter. There is also a non-specific degradation term. .

PARD is the antidote, thus under control of the promoter activated by the signal-complex, hence the first term. The next one is due to the complex formation of PARD with the poison, and the third one to the dissociation of the complex. The last term is of non-specific loss.

The antidote-poison complex is formed from these two molecules and degraded by their dissociation. There is also a non-specific degradation term.

Results.

Before we tried to answer the questions asked by the wetlab part, we tried to establish the different values of the parameters (and the order of the Hill's terms) allowing a satisfying modelisation of the module, and thus mathematically validing  the chosen system. These results will be introduced at the end of this section, because they are a summary of the modelized behaviour of the bacterial population, but we have to keep in mind that these stage is the very first one.
The first question answered a practical interest, because the main goal was to model the ideal relative value of RBS for each genes carried by plasmids. So we first simulated various configurations in this way.
We managed to vary the values of coefficients used in the above equations, which symbolized coded proteins on genes.
We studied the behavior of our population before and after a ten-order dilution as we were varying the RBS relative values (ones compare to each other). All parameters that gave an unstable population in the reactor or, to the contrary, stable or increasing population after dilution were immediately rejected.

Here are the main results obtained during our simulations:

Note that:
- According to the equations r1, r2, r3 represent the strength of "antidote" RBS and r4, r5 represent the strength of "poison" RBS.
- r1 represents LUXI RBS, r2 represents LUXR RBS, r3 represents PARD RBS, r4 represents HSLase RBS, r5 represents PARE RBS.
- We estimated that 6 units of time on the graphics represent 1 hour. We obtained this estimation by looking at the time a population needs to stabilize itself.

Step 1:
Vary R=r1=r2=r3=r4=r5

In the reactor:
For small values of R (<1 "relatively") the population collapses drastically after 5 or 6 hours spent in the reactor.
For higher values of R (>1) the population stabilizes itself. The bigger the value, the higher the population stabilizes itself.

In the river:
For small values of R (<1) the population is stable, which is a situation to absolutely avoid.
For higher values of R, we first observe a "baby-boom" quickly followed by a fall of population. (The importance of both of these phenomenons is proportional to R.)

A first result is obvious: we will have to work with powerful RBS.

Step 2:
Vary the poison (r4 = r5) and the antidote (r1 = r2 = r3) separately.
According to previous results, we will only use high values of RBS. (>1 "relatively")


Step 2.1 Improve the poison/antidote ratio

Results in the reactor:
- For a ratio of 1.05 (poison 5% stronger than antidote): during the 24 first hours, we only loose 5% of the population (compare to a ratio of 1), but later, the population does not stabilize itself.
- For a ratio of 1.1 (poison 10% stronger than antidote): during the 24 first hours, we loose 15% of the population (compare to a ratio of 1), and after 24h the population collapses.
- For a bigger ratio, the population collapses too quickly.
cf picture 1

Results in the river:
The population falls down more quickly when the ratio increases, but the difference is not that interesting.

Step 2.2 Improve the antidote/poison ratio

Results in the reactor:
- Increasing the ratio until 1.5 enables the population to reach a stable state at 92% of its initial situation. Compare to 85% for a ratio of 1. The bigger the ratio, the higher the population will stabilize itself.

Results in the river:
- If the ratio is over 1.5 the population will not decrease any more.
- It is better that we do not go over a ratio of 1.2 in order to avoid a huge baby-boom (inevitable if the ratio is bigger than 1).

Conclusion for step 2:

In the reactor:
If we increase the poison by 10% (max), it does not have significant effect while t<24h.
If we increase the antidote by 20%, the population stabilizes itself at 90% of its initial situation. Compare to 85% for a ratio of 1 . As shown on picture 2.

In the river:
If we stand in a bracket of +10% poison/+20% antidote, our population collapses in a significant way.

Step 3 Vary antidote ones compare to each other, and similarly for poisons.

If we keep in mind the last results (rules) we obtained (mainly, r>1 and the fact that ratio antidote/poison stands between 0.9 and 1.2),
the effects on the population are not that different from the previous cases, whatever the configuration is. So we can keep r1=r2=r3 and r4=r5
Our recommendations to the wetlab part were thus simple, design the genes with RBS of strength 1 (practically speaking, it's hard to find strong RBS which are only a little bit different). We then could start a more analytic part of our work.

Here, because of the order of Hill terms strictly upper than one in many equations of the module, we were hoping for special stability situations to occur and we tried to establish stability diagrams of the systems.

V1 is the key parameter of the antidote plasmid, since it's the coefficient standing for the transcription speed of the system. We generally have, in the absence of regulations systems or of special RBS value, arbitrarily attribute the value 1 to each parameter of that type. With this value v1=1, the population diluted 10 times collapses very fast. We then increased step-by-step that value of v1 to observe which value will lead to the stability of the population and if they exist intermediate population densities between stability and collapse of the population.

The result is a very strong step effect. Indeed, either the population remains up to 80% of its intial density, or it collapses completely. This is very clearly shown on graph1, where the curve in red and blue respectively represent the evolution of the population for the values of v1=1 and v1=1.19, when the green curve represents the evolution of the population over time if v1=1.2 .

In the same way, we fixed the value of v1 to 1 and studied the behavior of the population depending of its initial density. We can see (picture3) that for Ninitial=39, the population grows and stabilizes around 80%, when for Ninitial=38, this population density collapses. We still notice that the bacterial density, few instants after the beginning of the simulation, grows in a similar way to what happens for Ninitial=39. This is explained by the non-nul initial concentration of HSL which maintain the population even if the quorum sensing isn't sufficient to produce it during some time (see graph for the parallel between concentration in bacteria and HSL). We chose to give a non-null initial value to HSL because some is added in initials growth conditions and, if some bacteria escape in the river- which is similar to a dilution-, their intracellular rate of HSL won't be null. However, we can emphasize an important step effect and a bistability phenomenon.

This step effect is studied once more in the next graph (graph3) which shows the evolution of the minimal relative transcription speed of the antidote plasmid (v1) such that the population is stable, depending on the initial population density. We logically see the decrease of the v1 value when the initial density increases, but more unexpected, we notice a limite value of v1 (=0.985) such that, no matter how big the intial population is, this population finally collapses. This emphasizes the great importance not to use a too small value of RBS for this plasmid, to avoid the collapse of the bacterial population in the reactor.

Then, using the datas computed by the auto software, we draw the stability graph for N (population density, in % of the maximal value in reactor; this symbol will be used in all this section) for different values of v1 (graph4). We see here two disjoint stable curves; the first one, which is a constant equal to zero is the only one for v1<0.985, as expected, and continues for any value of v1. This is explained by the fact that, whatever the value of v1, it's always possible to find a Ninitial such that the population collapses. The second curve only starts at v1=0.985, and the equilibrium value of N increases with v1. This means that, whatever the initial value of the population density, if it's bigger that the step value, the equilibrium value of N will be the same, for a given v1.

Finally we did a little simulation of a real situation, with a population in the reactor initially equal to 100% and a dilution factor of ten, at t=10. The first graph (picture4) is quite clear: during the dilution, the density falls below the critical value, and the population then self-destructs. Zooming on different variables for the population (picture5) we notice the decrease in "growth" proteins (HSL, ParD, LuxR-HSL) and the increase in « destructive » proteins (HSLase, ParE), from the dilution time through a compensation point around 17, in which values are inverted which is the point of no return for the population and until there are no more growth proteins. The system chosen for the programmed cell death seems thus to work perfectly well.

Bacteria planned death system

Extra-Scientific issues