Team:Potsdam Bioware/Project/Details Modeling

Modeling

There is no synthetic biology without modeling, of course. In principle there is structure modeling and system modeling. In structure modeling the conformation and structure of proteins is examined and steric consequences for reactions or the whole system can be estimated. We focused on the second sort of modeling: The system modeling in which the reaction kinetics of the whole system is analyzed, outcomes are predicted and parameters correlated to measurements. Thus, a synthetic biology approach can be chosen because a better understanding of the system is achieved and further changes can be planed - just like in engineering.

Model

We built a model of our in vivo selection system to be able to engineer our system effectively. The following schema shows the major reactions taking place in our cell system. The Romanic numbers indicate the system relevant (partially triggered) reactions that were written down as equations and then numerically propagated through time.

Figure 1:Simplified schema of our cell system including labels and markers for the chemical reactions in it. Marked reactions are considered in our system modeling.

The schema shows the triggered expression of mircoviridin (our inhibitor), the protease (that needs to be inhibited for medical reasons) and the β-lactamase (that protects the cell from an antibiotic). The protease cleaves at a specific recognition site in the linker peptide. This will abolish the β-lactamase export and in consequnece the cell will die.
There are three important trigger activation times:

• (t0) - start (microviridin added already)
  
• t1 - start expression of protease
  
• t2 - start expression of β-lactamase
  
• t3 - Ampicillin added into medium
  
• (t4) - end of the experiment: cell cultures survive or die.
  

Keeping these times in mind, the reactions can be written down in seven chemical reaction equations of different sort and order.

Figure 2: Reaction equations between relevant molecules in the microviridin-inhibitor-concept including indications about the triggered time period (t1,t2,t3). The Romanic numbers correspond to the reactions marked in the above schema.

Concentration equations

From the above reaction equations differential equations can be derived that describe the change of substance concentrations. Three concentrations are fixed, however:

Figure 3: Given substance concentrations for differential equations. Those concentrations are added to the system and will not increase over time.

To calculate each concentration at a certain time, several parameters (also refered to as constants or factors) have to be estimated. A literature research was not able to provide us with exact values for our system but a typical range for comparable systems could be assigned to all values. Most of the parameter units had to be transformed. The robustnes of the system towards changes of the most important constants was examined (see results) and some parameters could be estimated by fitting the system to lab measurements (see parameter estimation).

All not-constant concentrations can be represented in form of differential equations. Between time t1 and t2 four parameters are introduced:

• k+1 in (1/s*molecules)- parameter for the association of MdnA (microviridin) and Prot (protease)
  
• k-1 in (1/s) - parameter for the dissociation of the inhibited protease
  
• kexpr.prot in (molecules/s) - parameter for the expression of protease
  
• kdeg1 in (1/s) - parameter for the degradation of protease

Figure 4: Differential equations for substance concentrations from t1 until t2. At t1 the expression of protease begins.

Between time t2 and t3 six additional parameters are introduced:

• k+2 in (1/s*molecules)- factor for association of Prot (protease) and TorABla (β-lactamase)
  
• k-2 in (1/s) - parameter for the dissociation of protease and substrate
  
• kcat in (1/s) – parameter for the catalytic enzyme reaction that cleaves the signale TorA sequence from the β-lactamase inactivating the export
  
• kexpr.Tor in (molecules/s) - parameter for the the expression of β-lactamase
  
• kdeg2 in (1/s) – parameter for the degradation of β-lactamase
  
• ktransTor (1/s) – parameter by which β-lactamase in the cytoplasm is able to pass the membrane and get into the periplasm

Figure 5: Differential equations for substance concentrations from t1 until t2. At t2 the expression of lactamase begins.

After t3 there are only two more factors that we should introduce:

• ktransAmp (1/s) – parameter by which the added Ampicillin in the medium is able to pass the outer membrane and get into the periplasm
  
• kcat2 in (1/s*molecules) – parameter for the catalytic enzyme reaction of ampicillin deactivation

Figure 6: Differential equations for substance concentrations from t1 until t2. At t1 ampicillin is added.

Results

The equations above were solved using MATLAB. WARNING: These equations are moderately to very stiff! A solution can only be obtained using the functions ode23t or ode23s!
Using the Avogadro constant and the volume of a E.coli cell (9*10^-16L) and a periplasm-volume of (3*10^-17L), following constants were calculated and compared with similar literature values:

• k+1 = 8.e-6 (1/molecules*s)
  
• k-1 = 4.e-4 (1/s)
  
• kexpr.Prot = 0.6 (molecules/s)
  
• kdegr1 = 17.e-4 (1/s)
  
• k+2 = 4.e-5 (1/molecules*s)
  
• k-2 = 2.e-4 (1/s)
  
• kcat = 8 (1/s)
  
• kexpr.Tor = 1 (1/s)
  
• kdegr2 = 2.e-3 (1/s)
  
• ktransTor = 1.e-3 (1/s)
  
• ktransAmp = 5.e-3 (1/s)
  
• Kcat2 = 0.001 (1/molecules*s)
  
• mdnA(0) = 3000 (molecules)
  
• Amp(0 Med) = 8.e+2 (molecules)
• alpha = 0.3163 (forumlar parameter)

Now the concentration equations could be visualized in all three time segments:

Figure 7a: MATLAB numerical solution to our system of differential equations for substance concentrations (above) over all time segments.

It can be seen that there is only a negligible amount of Ampicillin in the periplasm if the expression works as indicated. This amount of Ampicillin has only a smal effect on the growth rate of the cells. There is also a wide tolerance left for suboptimal conditions. The trigger time t1, t2 and t3 shall be about one hour after one another.
In the first section the effect of the inhibition can clearly be seen: Even though the continuous expression of protease, only a very small number remains active in the cytoplasm. In the second section β-lactamase is expressed. Most of the β-lactamase in the cytoplasm is destroyed by the protease, however the number is by far large enough to steadily release molecules into the periplasm where it cannot be affected by the protease any more. Because the volume of the periplasm is very small, the concentration there is even higher than the absolute amount of molecules in the graph suggests compared to the amount of molecules in the cytoplasm. The third section shows that the concentrations remain very steady and the cells are Ampicillin resistant.

In case microviridin is not expressed the simulation indicates following substance concentrations that were correlated to our lab data of the same system (see parameter estimation). Under those conditions cells die quickly.

Figure 7b: MATLAB numerical solution to our system of differential equations for substance concentrations under absense of Microviridin.

In case the protease is also not expressed the simulation indicates following substance concentrations. They were correlated to measurements as well. Protected by lactamase cells survive well on Apmicillin plates.

Figure 7c: MATLAB numerical solution to our system of differential equations for substance concentrations under absense of Protease.

Cells in which all parts of our system work are very Ampicillin resistant and grow fast: They double their cell volume about every 20 minutes. If an error would appear and the Ampicillin concentration inside the periplasma would increase over 2µg/ml, the growth rate would slow down drastically and the cells might die. This way defect cell cultures can easily be separated from the good ones.

We constructed a formula for the ratio of the growth rate (cells double during this timespan) to the Ampicillin concentration in n µg/ml out of data that we gained out of the lab. The cell doubling time can approximatelly be calculated by: T(growth) = 20min+10*2^(2*n-1). Analogous to this the cell surviving percentage can be estimated by: surviving cells in % = 25/(20min+10*2^(alpha*n-1)).

Here we display one example of how the change of a simulation constant impacts the growth of the cells and that our system is robust: The parameter K_D (Dissociation of microviridin and protease divided by the association) is central to our model and the lactamase concentration in the periplasm an indicator for the cell survial. To see the dependence of lactamase on this constant, following plot (figure 8) was created:

Figure 8: Starting at t3: Lactamase concentration inside the periplasm. A significant change of the simulation factor K_D (Dissociation of microviridin and protease divided by the association) can result in a lower lactamase concentration and thus less fitness of the cell (because the lactamase protects the cell from the antibiotic ampicillin).

Figure 9 shows the relating cell doubling time that can be measured in the lab.

Figure 9: Starting at t3: Growth (double) rate of cells dependent on the ampicillin concentration in the periplasm. A significant change of the simulation factor K_D (Dissociation of microviridin and protease divided by the association) can result in less fitness of the cell and a reduction of its growth rate.

The lactamase concentration in the periplasm can be very small. Thus for this concentration a semi-stochastig model was produced by handling this value as integer. Such graphs are less smooth (see figures 10 and 11).

Figure 10: Starting at t3: Semi-stochastic lactamase concentration inside the periplasm.

Figure 11: Starting at t3: Semi-stochastic growth (double) rate of cells dependent on the ampicillin concentration in the periplasm..

We were able to see that our system works very well in theory. We learned about correct time-scales for our triggering and we were able to identify expected cell-division rates as reference for lab work.

Parameter estimation

To be filled...

MATLAB code

(update uploaded soon)
Media:conzt1.m
Media:conzt2.m
Media:conzt3.m
Media:UP_plottfuntionendreiaufeinmal2.m
Media:UP_celldivision2.m