Team:Colombia/Modeling/Paramterers

Parameters of the equations

When we want to model a biological process, it is necessary to write the differential equation that model the system which requires a number of constants. If the real value for the constants are unknown, the system can not have any biological sense. These entries are called parameters and they are crucial elements in the mathematical model made in this project.

There are three possible ways to find this parameters:

i) Literature. There are a lot of studies trying to find biological parameters, such as basal rate of protein production, kinetic constants, Monod's constants, etc. Moreover, there are so many biological systems and only a few of them have been characterized. Thus, it is difficult to find the parameters that are needed.

ii) Experimental way. If an experiment is made using the biological system of interest, it is possible to find the parameters for the equations that models the whole system. For this project, it was necessary to model the biological system first. Thus, experiments couldn't be performed to find the constant for the differential equations.

iii) Screening of parameters. Sometimes we don't have the exact number that we need, but we have a rank where it could be or the parameter for a similar biological system, then we can perfom a screening of parameter, where we try to find the value that perfectly fits the reponse of our circuit.

How did we do it?

For Pest-busters we used options i) and iii). In the literature, it was possible to find the parameters for the complete CI activation box. For the rest of values, we found an approximately rank where they could be. Thus, we proceeded to step iii).

We used a lot of resources and methods to approximate the limits for this ranges. Here we show an explanation:

a.Basal levels of the proteins:

Determining a rank which this value could be, we described mathematically a very simple process. Suppose there is a usual differential equation:

Where αo is the basal level, γ is the protein parameter for degradation and P is the protein concentration. Thus, we can assume that in steady state conditions we have that dP/dt= 0:

But γ is approximately 1/T, then:

Establishing the rank, we assumed T as φ and we looked for a normal rank of concentration of a protein in a cell [http://bionumbers.hms.harvard.edu/bionumber.aspx?&id=104520&ver=10&trm=protein]. Another important aspect is that there were constitutive proteins and inducible proteins. Thus, we divided the basal levels in two class: i) one for basal proteins and ii) one for constitutive proteins. Finally, following ranges were defined:

Basal proteins range (LuxR, LuxI, HipB and Salicylic Acid):

Constitutive proteins range (HipA7,Sensor, Chitinase, Chitoporin and CBP):

By the other hand, the basal production of HipA has to be constitutive and inducible because the cell has to have three stages. The first stage is always sleep (constitutive) due to HipA effect. Once an impulse of chitin or 3-OH-PAME is presented, it is going to awake the cell due to concentration decrease of HipA . This is the second stage of the cell. The third stage is when cell is slept again which it is achieve through a positive control (inducible). The basal concentration used in the equations was the constitutive one because it is going to have more effect in the response.

b.Protein degradation:

The chosen unit of time was the bacteria life time because it represents when the proteins are diluted and decrease in the cell. We can establish easily the protein degradation using values related to life time of E.Coli. We defined the protein degradation for almost all values as follows:

Again, we found a special case. Degradation of HB is greater, then we assigned it a value of 4/φ.

c. Reaction constant:

This value can vary depending of the described process. Fortunately, we found in literature how to approximate this value [http://www.biokin.com/dynafit/scripting/html/node23.html#SECTION00431000000000000000]. However, we took the approximation and turned it into units described above.

d. Hill coefficients k:

We looked for some reference in different data sources and we found the coefficient of CI ([http://www.sciencemag.org/content/307/5717/1962.short Nitzal et al.2005]). Then, we decided to approximate a rank taking this value as a reference.

e. Hill coefficients n:

We based on the value "n" found for CI( [http://www.sciencemag.org/content/307/5717/1962.short NItzal et al.2005]). In this case, there was an analysis that took in account the number of molecules that interfere with gene activation.

f. Maximum cell production (β):

Considering the concentration of Rubisco, which is about 20 times the average concentration and rate of production of a protein, [http://bionumbers.hms.harvard.edu/bionumber.aspx?&id=107431&ver=1&trm=rubisco], we assign the limits for this range as follows.

CI PARAMETERS

All the parameters for CI activation system were found in the literature. ([http://www.sciencemag.org/content/307/5717/1962.short Nitzal et al.2005]). We only made a change of units.

Once the ranks of the parameters were set, we proceed to do the screening. We tried a lot of techniques but we failed. Here we present the steps for the successful one:

Step 1.Sensitivity Analysis

We asked ourselves for the weighted importance of each of these parameters to our design’s main outputs, namely Salicylic Acid, the toxin HipA7 and the antitoxin HipB. Hence, we came up with a sensitivity test in order to understand how every parameter adds up to our model. This test does not know about the optimal or most useful parameter values, but instead informs us of how and how much will our final response vary in function of each parameter.

The test consisted of the following: Knowing the ranges of the parameters (see above) we determined appropriate range divisions to obtain enough data for analysis (to be known as step size from here onwards). Then we iterated each parameter values while leaving the others fixed in our MATLAB code.

Note: the exact value of each parameter is not important. What matters is their relevance and how they change the response.

Here is an example:

Parameter: Chitinase leakage production

Range: 0.1-0.9 μM

Step size: 0.01

The diagram above shows that the chitinase leakage production does not have an important effect in salicylic acid, toxin, or antitoxin behaviour.

We make this procedure for all the parameters in both systems (Ralstonia and Rust), here are the results:

Rust:

The following table show us the different parameters involved in our coffee rust detection system and how they affect the final behaviour.

The following graphs represent some more relevant examples. It seems our coffee rust detection system is not directly sensible to some parameters. The blue line is the Salicylic Acid, the green is HipA7 and the red one is HipB:

Ralstonia

We expose the results obtained by the sensibility test for Ralstonia. The next table shows the information obtained by the test so that we can extract two important aspects we can stand out from our analysis; the first one is to distinguish between those parameters that, when varying, make a relevant change in the model, from those that don’t, and the second aspect is, having relevant parameters, the classification of them taking into account those ones that are directly or inversely proportional to the salicylic acid response denoting them respectively as “+” and “-“.

The following image shows the results obtained for almost all the relevant parameters graphing the response of the salicylic acid in function of the value of a determined parameter:

Now we have an idea of which parameters affect the response the most. Although we still don't know what is the exact value of the parameter we know how the system is supose to response.With this two things we can proceed to step 2.

Step 2. Optimization of parameters

Within all the ranges of the parameters there is a number that combined with the other parameters makes the system behave just like is expected. There are many possible combinations where this can happen, and is represented as an area where the system woks. The figure below shows an example of the acceptable area of two parameters.

In two dimensions this sees easy to identify but when we are dealing with more than 10 parameters, doing this by hand or with iterations in a computer takes a lot of time (weeks or even months).

The goal in this step is to find a point within this area and from here the screening can begin. And one way for doing this is by a process of optimization

A process of optimization takes a function and maximizes o minimizes changing some variable of it. The changing variables are named optimization variables and function is called objective function. If we want to have some limitations or behaves with the variable we can add restrictions to the system, and the function will be optimized without breaking this limitations.

In our case the objective function is a minimum difference of squares between the points of our expected behavior and the response of the equations with a set of optimization variable. When the distance between these points is minimum the parameters are found!!!

This optimization process was done discretizing the differential equations and putting them as restrictions, the same was done for the stable state concentration where the differential equation must be equal to zero when the chitin is zero. All was programed in the software Gamside.

Unfortunately the firsts results were not satisfactory, so we are now developing a new code to run.

3.Screening

If we have a set of parameters why do we do a screening?

In the last step we found a set of parameters that give the expected response of the system. But this is only a point... What does happen with the rest of points of the area? Does the optimization method “know” if these parameters are real or have biological sense? As long as we specified one single response for the parameters, we don’t know if there are other possible responses of valid solutions. Do they exist?

To answer these questions we performed a screening of the parameters. Beginning at the resulting point of the optimization, we started to move in all directions in a fixed range and then we examined when the acceptable area ends.

How much do we move, depends on the sensitivity analysis – If the response is not very sensible to the parameter, we take longer steps than the case in which the response is largely sensible.

Now we need to know if we are in the acceptable area or not. To solve this problem we defined some ranges where we want the Salicylic Acid to be depending on the impulse of chitin or 3-OH-PAME. The figure below shows this:

 -Note: Due to the fact that this takes a lot of computational time we are parallelizing the code 
 and we expect to run the tests shortly