EndreGarberg (Talk | contribs) |
|||
Line 161: | Line 161: | ||
limit set by our Human Practices research?</li> | limit set by our Human Practices research?</li> | ||
<li>2) Where are the bottlenecks in the pathway, and can we fix them?</li> | <li>2) Where are the bottlenecks in the pathway, and can we fix them?</li> | ||
− | <li>3) Which signal transduction steps are robust | + | <li>3) Which signal transduction steps are robust, and which are sensitive?</li> |
</ol> | </ol> | ||
<br> | <br> | ||
− | + | To answer these questions, we made use of <b>ordinary differential equations (ODEs)</b>: | |
Equations that illustrate how the state of a variable changes with time. | Equations that illustrate how the state of a variable changes with time. | ||
<br> | <br> | ||
<br> | <br> | ||
− | After working with ODEs for a while, we realized that | + | After working with ODEs for a while, we realized that the limits of using deterministic mathematics to model something as unpredictable |
− | as a cell. | + | as a cell. Therefore, we chose to utilize stochastic |
differential equations (SDEs). | differential equations (SDEs). | ||
<br> | <br> | ||
Line 181: | Line 181: | ||
<br> | <br> | ||
These are the kinds of biologically relevant questions that could be answered | These are the kinds of biologically relevant questions that could be answered | ||
− | with stochastic modeling, partially because SDEs | + | with stochastic modeling, partially because SDEs can represent inherent |
− | randomness present in all existing processes. | + | randomness present in all existing processes. Additionally, the SDEs allowed us to take our lack of knowledge of fine |
− | + | ||
bioprocesses into account. | bioprocesses into account. | ||
<h2 class="titles">How We Model</h2> | <h2 class="titles">How We Model</h2> | ||
Line 192: | Line 191: | ||
<br> | <br> | ||
<br> | <br> | ||
− | We needed to identify values of parameters for the models to | + | We needed to identify the values of our parameters of interest for the models to allow for realistic simulations to run. |
− | + | These values were found by researching literature, other models, and past iGEM wiki-pages. | |
The parameters were selected to be consistent among our different designs. | The parameters were selected to be consistent among our different designs. | ||
− | The parameters found were taken as rough estimates (see | + | The parameters found were taken as rough estimates (see figure 1), because |
− | performance of biological systems can vary for reasons | + | the performance of biological systems can vary for several reasons that have yet to be elucidated. |
However, rough estimates are enough to get great insights into the mechanics of | However, rough estimates are enough to get great insights into the mechanics of | ||
− | cellular pathways. | + | cellular pathways. Optimally, producing more precise laboratory data, and fitting the models hereto, would be the next step to obtain more realistic modeling results. |
− | + | ||
</div> | </div> | ||
<img src="https://static.igem.org/mediawiki/2020/2/2d/T--UCopenhagen--parameters.png"> | <img src="https://static.igem.org/mediawiki/2020/2/2d/T--UCopenhagen--parameters.png"> | ||
− | <figcaption>List of parameters used in our models. | + | <figcaption> |
− | </figcaption> | + | <b>Figure 1: List of parameters used in our models.</b> |
− | + | </figcaption> | |
<div class="txt-btwn"> | <div class="txt-btwn"> | ||
<h2>Ordinary Differential Equations (ODEs)</h2> | <h2>Ordinary Differential Equations (ODEs)</h2> | ||
− | To build our ODEs, | + | To build our ODEs, the below-stated three steps were followed: |
<ol> | <ol> | ||
<li>1) Formulate a diagram that specifies the key players (state variables) | <li>1) Formulate a diagram that specifies the key players (state variables) | ||
Line 216: | Line 214: | ||
</ol> | </ol> | ||
<h5>Step 1 - Making a Diagram</h5> | <h5>Step 1 - Making a Diagram</h5> | ||
− | In order to make a diagram, we needed to have | + | In order to make a diagram, we needed to have an idea of what our design should look like. For this, one could assume that: |
− | + | ||
<ul> | <ul> | ||
<li>a) The presence of a biomolecule in proximity to the cellular membrane | <li>a) The presence of a biomolecule in proximity to the cellular membrane | ||
− | lead to the binding of our receptor to | + | lead to the binding of our receptor of interest to its ligand/biomolecule, and then the |
binding of a coreceptor to this complex.</li> | binding of a coreceptor to this complex.</li> | ||
<li>b) Our receptor and coreceptor were both transmembrane proteins, | <li>b) Our receptor and coreceptor were both transmembrane proteins, | ||
Line 226: | Line 223: | ||
the extracellular association.</li> | the extracellular association.</li> | ||
<li>c) Their intracellular association was facilitated by the natural | <li>c) Their intracellular association was facilitated by the natural | ||
− | affinity of the two halves of a split-ubiquitin molecule ( | + | affinity of the two halves of a split-ubiquitin molecule (for reference, see design page). This association and reconstitution of a full ubiquitin |
− | + | ||
prompted its cleavage by naturally occurring deubiquitinases.</li> | prompted its cleavage by naturally occurring deubiquitinases.</li> | ||
<li>d) Such a cleaving event released a transcription factor that was | <li>d) Such a cleaving event released a transcription factor that was | ||
Line 236: | Line 232: | ||
From the points above, our key players (state variables) were identified as: | From the points above, our key players (state variables) were identified as: | ||
biomolecule, transmembrane receptor protein, transmembrane coreceptor protein, | biomolecule, transmembrane receptor protein, transmembrane coreceptor protein, | ||
− | free transcription factor, gene, mRNA, reporter protein. All other factors | + | free transcription factor, gene, mRNA, and reporter protein. All other factors |
were treated as parameters with a constant value. The interactions between | were treated as parameters with a constant value. The interactions between | ||
− | our variables | + | our variables are described in point a-f. Based on this information, |
− | we made a | + | we made a diagram to make the pathway more intuitively visible, as can be seen below. |
</div> | </div> | ||
− | |||
<img src="https://static.igem.org/mediawiki/2020/8/8d/T--UCopenhagen--pathway-1.png"> | <img src="https://static.igem.org/mediawiki/2020/8/8d/T--UCopenhagen--pathway-1.png"> | ||
− | <figcaption>Diagram of a signaling pathway of the design #1. The design #2 looks virtually identical. | + | <figcaption> |
− | </figcaption> | + | <b>Figure 2: Diagram of a signaling pathway of the design #1. The design #2 looks virtually identical.</b> |
− | + | </figcaption> | |
<div class="txt-btwn"> | <div class="txt-btwn"> | ||
<h5>Step 2 - Writing Interactions in a Quantitative Form</h5> | <h5>Step 2 - Writing Interactions in a Quantitative Form</h5> | ||
− | The next step was to find out | + | The next step was to find out how much product was produced from a given amount of reactant. As we were dealing with signaling pathways, our |
− | + | ||
concentrations were changing in accordance with the association and dissociation | concentrations were changing in accordance with the association and dissociation | ||
− | of the variables. | + | of the variables. The equation derived hereof can be found below: |
+ | |||
<h6>k1[A][B]=>[C]</h6> | <h6>k1[A][B]=>[C]</h6> | ||
− | where k1 | + | |
− | [C] | + | where k1 is a kinetic parameter, [A] and [B] are reactant concentrations, and |
+ | [C] is a product concentration. Note that such a reaction proceeds only in | ||
one direction. However, the reality of many reactions was better represented | one direction. However, the reality of many reactions was better represented | ||
by breakdown into a forward reaction and a backward reaction. This was because | by breakdown into a forward reaction and a backward reaction. This was because | ||
it was empirically observed that chemical reactions tend to reach a state | it was empirically observed that chemical reactions tend to reach a state | ||
of equilibrium with nonzero concentrations of reactants and products. This | of equilibrium with nonzero concentrations of reactants and products. This | ||
− | was understood as a situation described by | + | was understood as a situation described by: |
+ | |||
<h6>k_on[A][B]<=>k_off[C]</h6> | <h6>k_on[A][B]<=>k_off[C]</h6> | ||
− | where k_on and k_off | + | |
− | respectively) | + | where k_on and k_off are kinetic parameters (binding and unbinding constant, |
− | + | respectively). Note that the reaction | |
+ | proceeds in both directions. After some time, equilibrium is reached, which is described by: | ||
+ | |||
<h6>k=k_on/k_off=[C]/([A][B]).</h6> | <h6>k=k_on/k_off=[C]/([A][B]).</h6> | ||
+ | |||
<br> | <br> | ||
− | + | In practice, it was much easier to deduce | |
− | k | + | k compared to its constituents k_on and k_off, because k could be calculated in a |
− | straightforward manner from observed concentrations. | + | straightforward manner from observed concentrations. This was first addressed later. The equilibrium was now written as: |
− | + | ||
<h6>k_on * [biomolecule] * [receptor] * [coreceptor] <=> k_off * [complex]</h6> | <h6>k_on * [biomolecule] * [receptor] * [coreceptor] <=> k_off * [complex]</h6> | ||
+ | |||
However, the next reaction where deubiquitinase cleaved the reassociated | However, the next reaction where deubiquitinase cleaved the reassociated | ||
− | split-ubiquitin proceeded only in a forward direction – once the transcription | + | split-ubiquitin proceeded only in a forward direction, as it is irreversible – once the transcription |
− | factor | + | factor becomes released, it is never again part of the transmembrane protein. |
<h6>k_on * [complex] * [deubiquitinase] => [transcription factor]</h6> | <h6>k_on * [complex] * [deubiquitinase] => [transcription factor]</h6> | ||
The next step was transcription, which was a slightly more complicated | The next step was transcription, which was a slightly more complicated | ||
− | reaction | + | reaction, only proceeding in the forward direction. The kinetics of transcription |
− | + | follows a simple form of Hill equation that was written as <b>[A]/(Kd+[A])=>[B]</b> | |
− | where [A] and [B] | + | where [A] and [B] are concentrations and Kd is a dissociation constant |
(a measure of how tightly the transcription factor bound a promoter). This | (a measure of how tightly the transcription factor bound a promoter). This | ||
equation further needed to be multiplied by a specific rate of transcription | equation further needed to be multiplied by a specific rate of transcription | ||
− | <b>r=g*R/L</b> where g | + | <b>r=g*R/L</b> where g is a number of copies of a given gene, R is a general |
− | rate of transcription and L | + | rate of transcription and L is a gene length. |
<h6>rate * [transcription factor] / (Kd + [transcription factor]) => [mRNA]</h6> | <h6>rate * [transcription factor] / (Kd + [transcription factor]) => [mRNA]</h6> | ||
− | Translation rate was a simpler case | + | Translation rate was a simpler case, as knowing the concentration |
− | of mRNA and the specific translation rate was enough. The specific translation | + | of mRNA and the specific translation rate was enough to determine the translation rate. The specific translation |
− | rate | + | rate is r=R/L where R is a general rate of translation and L is the length of mRNA. |
<h6>rate * [mRNA] => [protein]</h6> | <h6>rate * [mRNA] => [protein]</h6> | ||
<h5>Step 3 - Converting quantified interactions into a mathematical framework.</h5> | <h5>Step 3 - Converting quantified interactions into a mathematical framework.</h5> | ||
− | + | Concerting quantified interactions into a mathematical framework is the first step that was model dependent. Here, we wanted to create a | |
model based on ordinary differential equations. Such equations described how | model based on ordinary differential equations. Such equations described how | ||
values of given variables change over time, or rather how they changed during | values of given variables change over time, or rather how they changed during | ||
Line 305: | Line 305: | ||
equation needed to capture all changes that happen to a single variable. | equation needed to capture all changes that happen to a single variable. | ||
This included both additions and subtractions – if a certain molecule was | This included both additions and subtractions – if a certain molecule was | ||
− | a product of one reaction but a reactant in other reaction, | + | a product of one reaction but a reactant in other reaction, the concentration of the molecule could both rise or drop. |
− | + | ||
<br> | <br> | ||
<br> | <br> | ||
Line 314: | Line 313: | ||
constituent parts. The third term captured spending of the complex by catalytic | constituent parts. The third term captured spending of the complex by catalytic | ||
activity of deubiquitinase.<br> | activity of deubiquitinase.<br> | ||
+ | |||
<h6>D(complex) = k_on * [biomolecule] * [receptor] * [coreceptor] – k_off * [complex] – k_cat * [complex] * [deubiquitinase]</h6> | <h6>D(complex) = k_on * [biomolecule] * [receptor] * [coreceptor] – k_off * [complex] – k_cat * [complex] * [deubiquitinase]</h6> | ||
− | + | We had now made a system of ordinary differential equations that described | |
− | kinetics of our signaling pathway. | + | kinetics of our signaling pathway. Next, we needed to translate our mathematical |
− | model into a programming language of choice, e.g. Matlab. | + | model into a programming language of choice, e.g. Matlab. Each variable was named ‘c’ and an arbitrary |
− | + | number. Parameters and whole equations were handled in a similar fashion: | |
− | number | + | |
<h6>d4 = k1 * c1 * c2 * c3 – k2 * c4 – k3* c4 * k4;</h6> | <h6>d4 = k1 * c1 * c2 * c3 – k2 * c4 – k3* c4 * k4;</h6> | ||
Except from providing initial conditions, we did not need to worry about all | Except from providing initial conditions, we did not need to worry about all | ||
− | the c’s and d’s any longer. The k’s however | + | the c’s and d’s any longer. The k’s, however, informed about |
− | the speed and order by which our system worked | + | the speed and order by which our system worked. |
− | + | ||
− | + | ||
<br><br> | <br><br> | ||
− | + | As scientific work must be translatable at all times, | |
we paid special attention to satisfying standards like MIRIAM. Knowing | we paid special attention to satisfying standards like MIRIAM. Knowing | ||
that many various software tools were used for pathway modeling, we | that many various software tools were used for pathway modeling, we | ||
utilized Moccasin to turn our Matlab code into SBML files that could be | utilized Moccasin to turn our Matlab code into SBML files that could be | ||
− | opened with a wide variety of programs. | + | opened with a wide variety of programs. Below are snippets of code from our models: |
− | + | ||
− | + | ||
</div> | </div> | ||
− | |||
<img src="https://static.igem.org/mediawiki/2020/9/99/T--UCopenhagen--ode-1.png"> | <img src="https://static.igem.org/mediawiki/2020/9/99/T--UCopenhagen--ode-1.png"> | ||
− | <figcaption>Example of the Matlab code from one of our models. | + | <figcaption> |
− | </figcaption> | + | <b>Figure 3: Example of the Matlab code from one of our models.</b> |
− | + | </figcaption> | |
<div clas="txt-btwn"> | <div clas="txt-btwn"> | ||
<br> | <br> | ||
Line 350: | Line 344: | ||
<br> | <br> | ||
<h2 id="SDE">Stochastic Differential Equations (SDEs)</h2> | <h2 id="SDE">Stochastic Differential Equations (SDEs)</h2> | ||
− | This part built upon the previous | + | This part built upon the previous section concerning ODEs. In fact, the ODEs described in the previous section were |
− | so similar that we took our implemented system of differential equations | + | so similar to the SDEs that we took our implemented system of differential equations |
− | and copied it in the right place of our code for the SDE-based model. | + | and copied it in the right place of our code for the SDE-based model. Before that, we simplified the equation of our mathematical model |
− | + | to the form where M denoted all the deterministic parts. In other words, M was all the | |
− | + | ||
− | + | ||
contribution coming from the signaling pathway of our interest. | contribution coming from the signaling pathway of our interest. | ||
<br> | <br> | ||
+ | |||
<h6>D(complex) = k_on * [biomolecule] * [receptor] * [coreceptor] – k_off * [complex] – k_cat * [complex] * [deubiquitinase] = M</h6><br> | <h6>D(complex) = k_on * [biomolecule] * [receptor] * [coreceptor] – k_off * [complex] – k_cat * [complex] * [deubiquitinase] = M</h6><br> | ||
− | Now, we needed a term that was driven entirely by randomness | + | |
+ | Now, we needed a term that was driven entirely by randomness. We called it N. | ||
If we decided that N stood for a stochastic part of the equation, we got an | If we decided that N stood for a stochastic part of the equation, we got an | ||
equation like this:<br> | equation like this:<br> | ||
+ | |||
<h6>D(complex) = M + N</h6><br> | <h6>D(complex) = M + N</h6><br> | ||
− | + | ||
− | + | In the context of SDEs, N represents a Wiener process. N is seen as a variable, whose value changes randomly but without any discontinuity occuring in its graph. | |
− | + | ||
<br> | <br> | ||
</div> | </div> | ||
<img src="https://static.igem.org/mediawiki/2020/2/2b/T--UCopenhagen--wiener-process.png"> | <img src="https://static.igem.org/mediawiki/2020/2/2b/T--UCopenhagen--wiener-process.png"> | ||
− | <figcaption>Illustrative example of continuous noise generated by | + | <figcaption> |
− | </figcaption> | + | <b>Figure 4: Illustrative example of continuous noise generated by Wiener process.</b> |
+ | </figcaption> | ||
<div class="txt-btwn"> | <div class="txt-btwn"> | ||
<br> | <br> | ||
− | This allowed us to model the effects of randomness and probabilities. | + | This allowed us to model the effects of randomness and probabilities. Going back to the mathematical formula, this was expanded further, so |
− | + | that we obtained a stochastic differential equation that we later implemented in Matlab. | |
− | + | ||
− | that we obtained a stochastic differential equation that we later implemented. | + | |
<h6>D(complex) = M + N = k1 * [biomolecule] * [receptor] * [coreceptor] - k2 * [complex] - k3 * [complex] * [deubiquitinase] + N</h6> | <h6>D(complex) = M + N = k1 * [biomolecule] * [receptor] * [coreceptor] - k2 * [complex] - k3 * [complex] * [deubiquitinase] + N</h6> | ||
+ | |||
We implemented our SDE-based model in Matlab. Yet, there was a limited built-in | We implemented our SDE-based model in Matlab. Yet, there was a limited built-in | ||
− | support for SDEs in Matlab, so we used a package SDETools that was available | + | support for SDEs in Matlab, so we used a package called SDETools that was available |
online. It separated the deterministic and stochastic parts of each equation | online. It separated the deterministic and stochastic parts of each equation | ||
− | so that the deterministic parts stayed grouped together, which was why we | + | so that the deterministic parts stayed grouped together, which was why we could |
− | easily | + | easily copy and paste our system of ODEs. The stochastic part was then defined |
separately. | separately. | ||
<br> | <br> | ||
<br> | <br> | ||
− | This was no issue because Wiener processes | + | This was no issue because Wiener processes can be defined in any number of |
dimensions, i.e. if we had ten differential equations, we simply defined N | dimensions, i.e. if we had ten differential equations, we simply defined N | ||
as a 10-dimensional process – one dimension per equation.<br> | as a 10-dimensional process – one dimension per equation.<br> | ||
<br> | <br> | ||
− | As | + | As mentioned earlier, our intention was to make the code readable and |
− | translatable, so we followed the same sets of standards for our SDE models. | + | translatable, so we followed the same sets of standards for our SDE models as for our ODE models. |
− | One exception was that no SBML | + | One exception was that no SBML file was generated as there was no software |
− | support for SDETools | + | support for SDETools due to it being an unofficial package. For the structure of our Matlab code, click here for <a href="https://static.igem.org/mediawiki/2020/8/82/T--UCopenhagen--sde-1.png">contact information</a>, here for <a href="https://static.igem.org/mediawiki/2020/c/c8/T--UCopenhagen--sde-2.png">parameters</a>, here for <a href="https://static.igem.org/mediawiki/2020/1/16/T--UCopenhagen--sde-3.png">deterministic equations</a>, <a href="https://static.igem.org/mediawiki/2020/8/88/T--UCopenhagen--sde-4.png">Wiener process</a>, <a href="https://static.igem.org/mediawiki/2020/4/4f/T--UCopenhagen--sde-5.png">vectors and options</a>, and finally, <a href="https://static.igem.org/mediawiki/2020/3/33/T--UCopenhagen--sde-6.png">running the code</a>. |
</div> | </div> | ||
<h5>To Future iGEMers!</h5> | <h5>To Future iGEMers!</h5> | ||
<img src="https://static.igem.org/mediawiki/2020/6/62/T--UCopenhagen--SDE_quote.png"> | <img src="https://static.igem.org/mediawiki/2020/6/62/T--UCopenhagen--SDE_quote.png"> | ||
<div class="txt-btwn"> | <div class="txt-btwn"> | ||
− | + | There are some limitations to SDETools. Using a programming | |
− | language called Julia together with available packages could give a | + | language called Julia together with available packages could give a better SDE |
− | modeling experience | + | modeling experience (our recommendation for future teams). |
<h4>Note On Alternative Modeling Paradigms</h4> | <h4>Note On Alternative Modeling Paradigms</h4> | ||
− | Apart from ODEs and SDEs, there were many other modeling approaches useful | + | Apart from ODEs and SDEs, there were many other modeling approaches that are useful |
− | for modeling of signaling pathways. The first one to mention | + | for modeling of signaling pathways. The first one to mention is |
− | partial differential equations ( | + | partial differential equations (PDEs). They are very similar to ODEs, but |
− | they | + | they also allow a variable change to be dependent on other variables, and not only |
on time. This extension was essential for modeling physical processes that | on time. This extension was essential for modeling physical processes that | ||
are located in space (e.g. diffusion). | are located in space (e.g. diffusion). | ||
<br> | <br> | ||
<br> | <br> | ||
− | There | + | There are more ways to model that are quite different from differential |
− | equations | + | equations. Some of them are rules-based or stochastic, such as: |
Gillespie algorithm, Markov chains, Petri nets. | Gillespie algorithm, Markov chains, Petri nets. | ||
<br> | <br> | ||
<br> | <br> | ||
− | Different algorithms | + | Different algorithms come with different limitations or assumptions (that |
− | + | may not be immediately visible) so it could be a good idea to model the | |
− | same system in several ways to make sure | + | same system in several ways to make sure that one's predictions are meaningful. |
− | + | ||
− | + | ||
<h4>Pheromone Cascade – Attractive but Complicated</h4> | <h4>Pheromone Cascade – Attractive but Complicated</h4> | ||
<div style="display: flex; flex-direction:row;"> | <div style="display: flex; flex-direction:row;"> | ||
<div style="display: flex; font-family: Avenir, Arial; font-size: 16px; flex: 2.5; align-items: center;"> | <div style="display: flex; font-family: Avenir, Arial; font-size: 16px; flex: 2.5; align-items: center;"> | ||
− | Taking inspiration from last year’s UCopenhagen iGEM team, we modeled the | + | Taking inspiration from last year’s UCopenhagen iGEM team, Ovulaid, we modeled the |
− | pheromone pathway in a similar way to | + | yeast pheromone pathway in a similar way to what the did (using only two equations to |
abstract the whole cascade) and provided it with a parameter fixed for | abstract the whole cascade) and provided it with a parameter fixed for | ||
experimental data found in literature. After some experimentation with our | experimental data found in literature. After some experimentation with our | ||
− | models | + | models, we became convinced that the reactions, in order to lead to empirically observed hypersensitivity, were necessarily |
− | + | more complex compared to what some papers claimed. | |
<br> | <br> | ||
<br> | <br> | ||
Notably, one-to-one reactions could not provide significant amplification | Notably, one-to-one reactions could not provide significant amplification | ||
− | from a theorical standpoint. After | + | from a theorical standpoint. After reading through literature, we found more detailed |
− | descriptions of the pheromone cascade | + | descriptions of the pheromone cascade which included a series of complexes |
and one-to-many interactions. In other words, a lot of amplification. As | and one-to-many interactions. In other words, a lot of amplification. As | ||
− | such, we | + | such, we became convinced that our model now more accurately represented the |
complexity found in the pheromone pathway. | complexity found in the pheromone pathway. | ||
</div> | </div> | ||
<div style="flex: 2; margin:2%;"> | <div style="flex: 2; margin:2%;"> | ||
− | + | <img style="width:100%;" src="https://static.igem.org/mediawiki/2020/5/5a/T--UCopenhagen--one-to-many.png"> | |
− | <figcaption>Comparison of amplification provided by a one-to-one relationship versus one-to-many relationship. | + | <figcaption> |
− | </figcaption> | + | <b>Figure 5: Comparison of amplification provided by a one-to-one relationship versus one-to-many relationship.</b> |
+ | </figcaption> | ||
+ | |||
</div> | </div> | ||
</div> | </div> | ||
</div> | </div> | ||
<div class="txt-btwn"> | <div class="txt-btwn"> | ||
− | We found a brilliant curated model on biomodels.com. We extracted the important | + | We found a brilliant curated model on <a href="https://www.ebi.ac.uk/biomodels/BIOMD0000000032">biomodels.com</a>. We extracted the important |
− | part of the cascade and added it in our model. | + | part of the cascade and added it in our model. This complicated our model, however, making the SDE implementation run into RAM memory issues. |
− | + | ||
− | + | ||
</div> | </div> | ||
<img src="https://static.igem.org/mediawiki/2020/5/53/T--UCopenhagen--pathway-3.png"> | <img src="https://static.igem.org/mediawiki/2020/5/53/T--UCopenhagen--pathway-3.png"> | ||
− | <figcaption>Diagram of a signaling pathway of the design #3. | + | <figcaption> |
− | </figcaption> | + | <b>Figure 6: Diagram of a signaling pathway of the design #3.</b> |
+ | </figcaption> | ||
+ | |||
<div class="txt-btwn"> | <div class="txt-btwn"> | ||
<h3>Impact on Wet Lab - Experiments</h3> | <h3>Impact on Wet Lab - Experiments</h3> | ||
<h4 style="margin-top: 0%;">How our modeling was tied to wet lab experiments</h4> | <h4 style="margin-top: 0%;">How our modeling was tied to wet lab experiments</h4> | ||
− | Since the beginning of our journey, wet lab and dry lab worked closely together. | + | Since the beginning of our iGEM journey, wet lab and dry lab worked closely together. |
Thanks to our supervisors, we always intended to use our dry lab activities | Thanks to our supervisors, we always intended to use our dry lab activities | ||
to guide our wet lab decisions. As our project ended up focusing a great deal | to guide our wet lab decisions. As our project ended up focusing a great deal | ||
Line 466: | Line 461: | ||
<br><br> | <br><br> | ||
The dry lab guided the wet lab in finding out which of the several biosensor | The dry lab guided the wet lab in finding out which of the several biosensor | ||
− | designs was the most viable, in shaping the design of the assays for our | + | designs was the most viable, in shaping the design of the assays for our <em>Saccharomyces cerevisiae</em> |
cells, and also in engineering the proteins with the desired properties. | cells, and also in engineering the proteins with the desired properties. | ||
<br> | <br> | ||
− | + | In all the following experiments, all values were held equal except from the one studied at a time. | |
− | + | ||
− | + | ||
− | + | ||
− | In all the following experiments, all values were held equal except from the | + | |
<h4>Our ODE Models</h4> | <h4>Our ODE Models</h4> | ||
<h5 style="margin-top: -1%;">Dynamic Ranges</h5> | <h5 style="margin-top: -1%;">Dynamic Ranges</h5> | ||
Line 479: | Line 470: | ||
<div style="flex: 2; margin:2%;"> | <div style="flex: 2; margin:2%;"> | ||
<img style="width:100%; margin:2%;" src="https://static.igem.org/mediawiki/2020/3/3f/T--UCopenhagen--dynamic-ranges.png"> | <img style="width:100%; margin:2%;" src="https://static.igem.org/mediawiki/2020/3/3f/T--UCopenhagen--dynamic-ranges.png"> | ||
− | <figcaption>Comparison of dynamic ranges of our designs. Note the logarithmic scale on the x-axis. | + | <figcaption> |
− | </figcaption> | + | <b>Figure 7: Comparison of dynamic ranges of our designs. Note the logarithmic scale on the x-axis.</b> |
+ | </figcaption> | ||
+ | |||
</div> | </div> | ||
<div style="display: flex; font-family: Avenir, Arial; font-size: 16px; flex: 2.5; align-items: center;"> | <div style="display: flex; font-family: Avenir, Arial; font-size: 16px; flex: 2.5; align-items: center;"> | ||
Line 504: | Line 497: | ||
<div style="flex: 2; margin:2%;"> | <div style="flex: 2; margin:2%;"> | ||
<img style="width:100%; margin:2%;" src="https://static.igem.org/mediawiki/2020/c/ce/T--UCopenhagen--sensitivity-1.png"> | <img style="width:100%; margin:2%;" src="https://static.igem.org/mediawiki/2020/c/ce/T--UCopenhagen--sensitivity-1.png"> | ||
− | <figcaption>Sensitivity analysis of the design #1 performed in SimBiology. The reporter was most sensitive to: affinity of a transcription factor to a promoter (k15), interleukin concentration (c_1). | + | <figcaption> |
− | </figcaption> | + | <b>Figure 8: Sensitivity analysis of the design #1 performed in SimBiology. The reporter was most sensitive to: affinity of a transcription factor to a promoter (k15), interleukin concentration (c_1).</b> |
+ | </figcaption> | ||
+ | |||
</div> | </div> | ||
</div> | </div> | ||
Line 519: | Line 514: | ||
us in pathway optimization. | us in pathway optimization. | ||
<img style="width:100%;" src="https://static.igem.org/mediawiki/2020/d/df/T--UCopenhagen--correlation-3.png"> | <img style="width:100%;" src="https://static.igem.org/mediawiki/2020/d/df/T--UCopenhagen--correlation-3.png"> | ||
− | <figcaption>Correlation matrix of all proteins in the signaling pathway of the design #3 and also time. The blue | + | <figcaption> |
− | </figcaption> | + | <b>Figure 9: Correlation matrix of all proteins in the signaling pathway of the design #3 and also time. The blue is associated with a negative correlation (the values closer to -1 signified a higher negative correlation - if one concentration went up, the other went down), the red was associated with a positive correlation (the values closer to +1 signified a higher positive correlation - if one concentration went up, so did the other). Values closer to zero signified less correlation between two concentrations.</b> |
+ | </figcaption> | ||
+ | |||
</div> | </div> | ||
<div class="txt-btwn"> | <div class="txt-btwn"> | ||
Line 527: | Line 524: | ||
<div style="display: flex; font-family: Avenir, Arial; font-size: 16px; flex: 2.5; align-items: center;"> | <div style="display: flex; font-family: Avenir, Arial; font-size: 16px; flex: 2.5; align-items: center;"> | ||
Affinity of a transcription factor (TF) to a promoter was varied in | Affinity of a transcription factor (TF) to a promoter was varied in | ||
− | this case study. High-affinity TF-promoter pair lead to an increased | + | this case study. High-affinity TF-promoter pair (Kd = 1 nM) lead to an increased |
sensitivity accompanied by significantly higher reporter concentrations, | sensitivity accompanied by significantly higher reporter concentrations, | ||
− | compared with a low-affinity TF-promoter pair. | + | compared with a low-affinity TF-promoter pair (Kd = 100 nM). The biosensor’s dynamic range |
− | as well as reporter concentrations can be adjusted by | + | as well as reporter concentrations can be adjusted by the right choice of |
a TF-promoter pair. | a TF-promoter pair. | ||
</div> | </div> | ||
<div style="flex: 2; margin:2%;"> | <div style="flex: 2; margin:2%;"> | ||
<img style="width:100%; margin:2%;" src="https://static.igem.org/mediawiki/2020/6/6b/T--UCopenhagen--tf-comparison.png"> | <img style="width:100%; margin:2%;" src="https://static.igem.org/mediawiki/2020/6/6b/T--UCopenhagen--tf-comparison.png"> | ||
− | <figcaption>Higher affinity of a transcription factor to a promoter lead to a higher sensitivity and higher reporter concentrations. | + | <figcaption> |
− | </figcaption> | + | <b>Figure 10: Higher affinity of a transcription factor to a promoter lead to a higher sensitivity and higher reporter concentrations.</b> |
+ | </figcaption> | ||
+ | |||
</div> | </div> | ||
</div> | </div> | ||
Line 544: | Line 543: | ||
<div style="flex: 2; margin:2%;"> | <div style="flex: 2; margin:2%;"> | ||
<img style="width:100%; margin:2%;" style="width:100%;" src="https://static.igem.org/mediawiki/2020/5/5a/T--UCopenhagen--inflammation-change.png"> | <img style="width:100%; margin:2%;" style="width:100%;" src="https://static.igem.org/mediawiki/2020/5/5a/T--UCopenhagen--inflammation-change.png"> | ||
− | <figcaption>High levels of inflammation lead to a higher reporter concentration, if compared to baseline inflammation levels. | + | |
− | </figcaption> | + | <figcaption> |
+ | <b>Figure 11: High levels of inflammation lead to a higher reporter concentration, if compared to baseline inflammation levels.</b> | ||
+ | </figcaption> | ||
+ | |||
</div> | </div> | ||
<div style="display: flex; font-family: Avenir, Arial; font-size: 16px; flex: 2.5; align-items: center;"> | <div style="display: flex; font-family: Avenir, Arial; font-size: 16px; flex: 2.5; align-items: center;"> | ||
Line 551: | Line 553: | ||
In this case study, the design #3 was subjected to low and high | In this case study, the design #3 was subjected to low and high | ||
interleukin concentration (i.e. levels of inflammation). | interleukin concentration (i.e. levels of inflammation). | ||
− | High inflammation lead to a fast development of high reporter | + | High inflammation (c = 10 nM) lead to a fast development of high reporter |
− | concentrations, compared with low inflammation. Significant differences | + | concentrations, compared with low inflammation (c = 1 pM). Significant differences |
of the plotted curves suggested the design #3 could distinguish between | of the plotted curves suggested the design #3 could distinguish between | ||
low and high inflammation, as long as the measurements were done at the | low and high inflammation, as long as the measurements were done at the | ||
Line 558: | Line 560: | ||
<br><br><br> | <br><br><br> | ||
<b>This fits perfectly with our Human Practices considerations, where | <b>This fits perfectly with our Human Practices considerations, where | ||
− | multiple experts preferred a design that provided | + | multiple experts preferred a design that provided fast results, and |
patients were willing to wear the patch for up to 24 hours, according to | patients were willing to wear the patch for up to 24 hours, according to | ||
our survey results!</b> | our survey results!</b> | ||
Line 572: | Line 574: | ||
concentrations. There were many versions of TEV proteases available | concentrations. There were many versions of TEV proteases available | ||
with various levels of proteolytic activity. Protease choice had a | with various levels of proteolytic activity. Protease choice had a | ||
− | significant impact on the dynamic range of a biosensor. | + | significant impact on the dynamic range of a biosensor. See figure 12 for comparison of high proteolytic actvity (k_cat = 0.18 s^-1) and low proteolytic actvity (k_cat = 0.000907 s^-1). |
</div> | </div> | ||
<div style="flex: 2; margin:2%;"> | <div style="flex: 2; margin:2%;"> | ||
<img style="width:100%; margin:2%;" src="https://static.igem.org/mediawiki/2020/c/c8/T--UCopenhagen--tev-comparison.png"> | <img style="width:100%; margin:2%;" src="https://static.igem.org/mediawiki/2020/c/c8/T--UCopenhagen--tev-comparison.png"> | ||
− | <figcaption>Higher proteolytic activity of a TEV protease lead to a higher sensitivity of the biosensor. | + | <figcaption> |
− | </figcaption> | + | <b>Figure 12: Higher proteolytic activity of a TEV protease lead to a higher sensitivity of the biosensor.</b> |
+ | </figcaption> | ||
+ | |||
</div> | </div> | ||
</div> | </div> | ||
Line 589: | Line 593: | ||
<h5>Rapid aging lead to disorder in the cell</h5> | <h5>Rapid aging lead to disorder in the cell</h5> | ||
This case study examined the impact of noise increasing with time. Such a | This case study examined the impact of noise increasing with time. Such a | ||
− | phenomenon could be present in senescent cells, where | + | phenomenon could be present in senescent cells, where by-products of |
− | aging | + | aging obstruct the reliability of cellular processes. If such divergence of |
− | curves was observed in our sweat patches, we would have known the | + | curves was observed in our sweat patches, we would have known whether the |
dysfunctionality is primarily time dependent. | dysfunctionality is primarily time dependent. | ||
</div> | </div> | ||
<img src="https://static.igem.org/mediawiki/2020/4/40/T--UCopenhagen--aging-disorder.png"> | <img src="https://static.igem.org/mediawiki/2020/4/40/T--UCopenhagen--aging-disorder.png"> | ||
− | <figcaption>Amount of noise increasing with time (illustrative). | + | |
− | </figcaption> | + | <figcaption> |
+ | <b>Figure 13: Amount of noise increasing with time (illustrative).</b> | ||
+ | </figcaption> | ||
+ | |||
<div class="txt-btwn"> | <div class="txt-btwn"> | ||
− | <h5>Longer | + | <h5>Longer measurement allows the signal pathway to yield more consistent results</h5> |
This case study examined the impact of noise decreasing with time. Such a | This case study examined the impact of noise decreasing with time. Such a | ||
phenomenon could be observed in any cell, especially soon after triggering | phenomenon could be observed in any cell, especially soon after triggering | ||
the receptor system before the pathway processes overcame the magnitude of | the receptor system before the pathway processes overcame the magnitude of | ||
the biochemical background. If such divergence of curves was observed in our | the biochemical background. If such divergence of curves was observed in our | ||
− | sweat patches, we would | + | sweat patches, we would know that this is a perfectly expected behavior |
coming from the physical law that longer measurements yield more precise results. | coming from the physical law that longer measurements yield more precise results. | ||
</div> | </div> | ||
<img src="https://static.igem.org/mediawiki/2020/f/f2/T--UCopenhagen--convergence.png"> | <img src="https://static.igem.org/mediawiki/2020/f/f2/T--UCopenhagen--convergence.png"> | ||
− | <figcaption>Amount of noise decreasing with time (illustrative). | + | <figcaption> |
− | </figcaption> | + | <b>Figure 14: Amount of noise decreasing with time (illustrative).</b> |
+ | </figcaption> | ||
+ | |||
<div class="txt-btwn"> | <div class="txt-btwn"> | ||
− | <h5>Unknown process | + | <h5>Unknown process decreases and increases disorder in the cell</h5> |
− | This case study examined rather | + | This case study examined a rather exotic-looking hypothetical scenario where high levels |
of noise diminished and came back to initial levels afterwards. Such bizarre | of noise diminished and came back to initial levels afterwards. Such bizarre | ||
− | divergence of curves could | + | divergence of curves could indicate that the cellular pathway only worked under |
very specific conditions that were achieved only for a short time. | very specific conditions that were achieved only for a short time. | ||
</div> | </div> | ||
<img src="https://static.igem.org/mediawiki/2020/3/3c/T--UCopenhagen--decrease-and-increase.png"> | <img src="https://static.igem.org/mediawiki/2020/3/3c/T--UCopenhagen--decrease-and-increase.png"> | ||
− | <figcaption>Decreasing and increasing amount of noise (illustrative) | + | <figcaption> |
− | </figcaption> | + | <b>Figure 15: Decreasing and increasing amount of noise (illustrative)</b> |
+ | </figcaption> | ||
+ | |||
<div class="txt-btwn"> | <div class="txt-btwn"> | ||
<h5>Natural variability in measurements</h5> | <h5>Natural variability in measurements</h5> | ||
Line 627: | Line 638: | ||
</div> | </div> | ||
<img src="https://static.igem.org/mediawiki/2020/3/32/T--UCopenhagen--natural-variability.png"> | <img src="https://static.igem.org/mediawiki/2020/3/32/T--UCopenhagen--natural-variability.png"> | ||
− | <figcaption>Uniform amount of noise (illustrative) | + | <figcaption> |
− | </figcaption> | + | <b>Figure 16: Uniform amount of noise (illustrative)</b> |
+ | </figcaption> | ||
+ | |||
<div class="txt-btwn"> | <div class="txt-btwn"> | ||
<h5>Receptors are not specific enough</h5> | <h5>Receptors are not specific enough</h5> | ||
Line 638: | Line 651: | ||
</div> | </div> | ||
<img src="https://static.igem.org/mediawiki/2020/7/71/T--UCopenhagen--non-specific-receptors.png"> | <img src="https://static.igem.org/mediawiki/2020/7/71/T--UCopenhagen--non-specific-receptors.png"> | ||
− | <figcaption>Noise applied on the receptors (illustrative). | + | <figcaption> |
− | </figcaption> | + | <b>Figure 17: Noise applied on the receptors (illustrative).</b> |
+ | </figcaption> | ||
+ | |||
<div class="txt-btwn"> | <div class="txt-btwn"> | ||
<h5>Receptor disrupts cellular processes</h5> | <h5>Receptor disrupts cellular processes</h5> | ||
− | This case study examined impact of noise increase tied to reporter | + | This case study examined the impact of noise increase tied to reporter |
concentration. If such divergence of curves was observed in our sweat | concentration. If such divergence of curves was observed in our sweat | ||
− | patches, we would have known the reporter was causing significant | + | patches, we would have known that the reporter was causing significant |
− | harm to reliability of cellular processes. | + | harm to the reliability of cellular processes. |
</div> | </div> | ||
<img src="https://static.igem.org/mediawiki/2020/6/6f/T--UCopenhagen--reporter-disrupts.png"> | <img src="https://static.igem.org/mediawiki/2020/6/6f/T--UCopenhagen--reporter-disrupts.png"> | ||
− | <figcaption>Amount of noise increasing with reporter concentration (illustrative). | + | <figcaption> |
− | </figcaption> | + | <b>Figure 18: Amount of noise increasing with reporter concentration (illustrative).</b> |
+ | </figcaption> | ||
+ | |||
<div class="txt-btwn"> | <div class="txt-btwn"> | ||
<h5>Variability in Interleukin Levels</h5> | <h5>Variability in Interleukin Levels</h5> | ||
− | This case study examined impact of variability in interleukin concentration. | + | This case study examined the impact of variability in interleukin concentration. |
This behavior could be expected in our sweat patch as the baseline interleukin | This behavior could be expected in our sweat patch as the baseline interleukin | ||
− | concentrations | + | concentrations likely varies among people. If such divergence of curves was |
− | observed in our sweat patches, we would have known this behavior was expected | + | observed in our sweat patches, we would have known that this behavior was expected |
and we could deduce variability among people. | and we could deduce variability among people. | ||
</div> | </div> | ||
<img src="https://static.igem.org/mediawiki/2020/3/34/T--UCopenhagen--variable-interleukin.png"> | <img src="https://static.igem.org/mediawiki/2020/3/34/T--UCopenhagen--variable-interleukin.png"> | ||
− | <figcaption>Noise applied on interleukin concentrations (illustrative). | + | <figcaption> |
− | </figcaption> | + | <b>Figure 19: Noise applied on interleukin concentrations (illustrative).</b> |
+ | </figcaption> | ||
<div class="txt-btwn"> | <div class="txt-btwn"> | ||
<h5>Diminishing disorder of cellular processes</h5> | <h5>Diminishing disorder of cellular processes</h5> | ||
− | This case study examined effect of initially high but diminishing noise. | + | This case study examined the effect of initially high but diminishing noise. |
This exotic scenario could represent the cellular process stabilizing at | This exotic scenario could represent the cellular process stabilizing at | ||
− | random states. Such a behavior | + | random states. Such a behavior is not expected in a functioning biosensor |
and could indicate there was a serious issue with part of the signaling pathway, | and could indicate there was a serious issue with part of the signaling pathway, | ||
− | or perhaps the cells had been compromised. | + | or perhaps that the cells had been compromised. |
</div> | </div> | ||
<img src="https://static.igem.org/mediawiki/2020/b/bd/T--UCopenhagen--stabilization.png"> | <img src="https://static.igem.org/mediawiki/2020/b/bd/T--UCopenhagen--stabilization.png"> | ||
− | <figcaption>Decreasing amount of noise (illustrative). | + | <figcaption> |
− | </figcaption> | + | <b>Figure 20: Decreasing amount of noise (illustrative).</b> |
+ | </figcaption> | ||
+ | |||
<div class="txt-btwn"> | <div class="txt-btwn"> | ||
− | <h2 id="protein">G-alpha and Structural Modelling | + | <h2 id="protein">G-alpha and Structural Modelling with Rosetta</h2> |
<div style="display:flex; flex-direction:row;"> | <div style="display:flex; flex-direction:row;"> | ||
<div style="display: flex; font-family: Avenir, Arial; font-size: 16px; flex: 2.5; align-items: center;"> | <div style="display: flex; font-family: Avenir, Arial; font-size: 16px; flex: 2.5; align-items: center;"> | ||
Line 680: | Line 700: | ||
natural function when there was no signal (i.e. no cytokine), and also dissociated | natural function when there was no signal (i.e. no cytokine), and also dissociated | ||
appropriately in the presence of a signal.<br><br> | appropriately in the presence of a signal.<br><br> | ||
− | Following the recommendation from team Aalto-Helsinki, we looked into Rosetta | + | Following the recommendation from team Sinisens (iGEM Aalto-Helsinki 2020), we looked into Rosetta |
Software Suite. Alanine scan was a Rosetta script available via an old version | Software Suite. Alanine scan was a Rosetta script available via an old version | ||
of Robetta server, and we used it to calculate how Gibbs free energy of the | of Robetta server, and we used it to calculate how Gibbs free energy of the | ||
two proteins (G-alpha and G-beta) changed upon changing a single residue to | two proteins (G-alpha and G-beta) changed upon changing a single residue to | ||
− | alanine. Changes in Gibbs free energy | + | alanine. Changes in Gibbs free energy correlate with binding affinity.<br><br> |
We submitted a list of residues selected according to our previous analyses. | We submitted a list of residues selected according to our previous analyses. | ||
− | + | A crystal structure of the bovine G protein complex served as a template for spatial | |
structural positioning of our proteins.<br> | structural positioning of our proteins.<br> | ||
</div> | </div> | ||
<div style="flex: 2; margin:2%;"> | <div style="flex: 2; margin:2%;"> | ||
<img style="width: 100%; margin:2%;" src="https://static.igem.org/mediawiki/2020/1/11/T--UCopenhagen--Alanine_Scan.png"> | <img style="width: 100%; margin:2%;" src="https://static.igem.org/mediawiki/2020/1/11/T--UCopenhagen--Alanine_Scan.png"> | ||
− | <figcaption> | + | <figcaption> |
+ | <b>Figure 21: Results of the alanine scan for select residues in G protein alpha subunit.</b> | ||
</figcaption> | </figcaption> | ||
+ | |||
</div> | </div> | ||
</div> | </div> | ||
− | Thanks to Kurt V. Mikkelsen, we got the opportunity to further our simulations | + | Thanks to Kurt V. Mikkelsen, quantum chemistry researcher from Niels Bohr Institute in Copenhagen, and |
+ | Andreas Erbs Hillers-Bendtsen, MSc student of quantum chemistry, we got the opportunity to further our simulations | ||
with Rosetta. After several engineering cycles, we crafted seven G-alpha mutants. | with Rosetta. After several engineering cycles, we crafted seven G-alpha mutants. | ||
We wanted to find out which one was the best candidate for our purposes. We wanted | We wanted to find out which one was the best candidate for our purposes. We wanted | ||
− | to identify a candidate mutant G-alpha with binding affinity | + | to identify a candidate mutant G-alpha with binding affinity to G-beta similar to wild type G-alpha, |
− | but with lowest possible binding affinity of its fragments after cleavage.<br> | + | but with the lowest possible binding affinity to G-beta of its fragments after cleavage by TEV protease.<br> |
Thus, we needed to simulate binding affinities (or change in Gibbs free energy | Thus, we needed to simulate binding affinities (or change in Gibbs free energy | ||
− | – delta delta G) of seven G- | + | – delta delta G) of seven G-alpha proteins and nine fragments formed after cleavage. |
− | We eventually landed on using a Rosetta script flex_ddg, as this script | + | We eventually landed on using a Rosetta script called flex_ddg, as this script has been |
shown to be less computationally demanding than other scripts, while making great | shown to be less computationally demanding than other scripts, while making great | ||
predictions of molecular dynamics<a href="#rosetta" aria-describedby="footnote-label" id="rosetta-ref"> </a>. | predictions of molecular dynamics<a href="#rosetta" aria-describedby="footnote-label" id="rosetta-ref"> </a>. | ||
<br> | <br> | ||
<br> | <br> | ||
− | We particularly | + | We particularly found the script that was written with a capability to include |
− | refinement of every protein structure many times | + | refinement of every protein structure many times to be useful. We followed recommendations |
− | and went with 35 000 | + | from script creators and went with 35 000 independent refinements per structure. Although the script was initially developed for |
− | mutagenesis simulations, we found | + | mutagenesis simulations, we found an easy way to adjust it for our own use. |
For each of our 16 proteins (mutants + fragments), we generated 35 independently | For each of our 16 proteins (mutants + fragments), we generated 35 independently | ||
optimized structures to probe the conformational space (this took 20 days on 16 | optimized structures to probe the conformational space (this took 20 days on 16 | ||
CPUs with 64 GBs of RAM). We used the scoring function REF2015 as it was shown | CPUs with 64 GBs of RAM). We used the scoring function REF2015 as it was shown | ||
− | to evaluate results of flex_ddg better than others. | + | to evaluate results of flex_ddg better than others (https://pubs.acs.org/doi/10.1021/acscentsci.8b00717). |
</div> | </div> | ||
<div style="display: flex; flex-direction: row; max-width: 95%; margin: 0 auto; justify-content: space-between;"> | <div style="display: flex; flex-direction: row; max-width: 95%; margin: 0 auto; justify-content: space-between;"> | ||
<img style="max-width: 45%; align-self: center;" src="https://static.igem.org/mediawiki/2020/f/fe/T--UCopenhagen--fragments-affinity.png"> | <img style="max-width: 45%; align-self: center;" src="https://static.igem.org/mediawiki/2020/f/fe/T--UCopenhagen--fragments-affinity.png"> | ||
+ | <figcaption> | ||
+ | <b>Figure 22: Change in Gibbs free energy of post-cleavage fragments of G-alpha informs us about their binding affinity. Number of residues on X-axis serves as a convenient way to observe whether there is a correlation.</b> | ||
+ | </figcaption> | ||
+ | |||
<img style="max-width: 45%; align-self: center;" src="https://static.igem.org/mediawiki/2020/b/ba/T--UCopenhagen--mutants-affinity.png"> | <img style="max-width: 45%; align-self: center;" src="https://static.igem.org/mediawiki/2020/b/ba/T--UCopenhagen--mutants-affinity.png"> | ||
+ | <figcaption> | ||
+ | <b>Figure 23: Change in Gibbs free energy of G-alpha mutants (labeled according to what locations contain a cleavage site).</b> | ||
+ | </figcaption> | ||
</div> | </div> | ||
<div class="txt-btwn"> | <div class="txt-btwn"> | ||
− | Results of the simulation were | + | Results of the simulation were surprising as they revealed our |
reasoning might have been flawed. Essentially, our assumption had been that | reasoning might have been flawed. Essentially, our assumption had been that | ||
− | cleaved fragments dissociate from the protein, yet our results | + | cleaved fragments dissociate from the protein, yet our results, as can be observed in figures 22 and 23, |
− | showed that majority of protein fragments had a higher affinity than parent proteins. | + | showed that the majority of protein fragments had a higher affinity than the parent proteins. |
We also realized that our simulation was not directly answering our true question | We also realized that our simulation was not directly answering our true question | ||
– could cleavage of our mutant G-alpha trigger Ste5 recruitment by G-beta? | – could cleavage of our mutant G-alpha trigger Ste5 recruitment by G-beta? | ||
Directly simulating that could include up to seven simultaneously interacting | Directly simulating that could include up to seven simultaneously interacting | ||
− | proteins, | + | proteins, which translated into a very demanding computational task.<br><br> |
After analyzing the results, we decided to run another simulation with the same | After analyzing the results, we decided to run another simulation with the same | ||
script but different input, to focus on fragments of the best-scoring protein | script but different input, to focus on fragments of the best-scoring protein | ||
– M124 (GPA mutant with TEV sites numbered 1, 2, and 4). What made this simulation | – M124 (GPA mutant with TEV sites numbered 1, 2, and 4). What made this simulation | ||
− | different and more precise that the previous one | + | different and more precise that the previous one was a ground reference for the |
calculation of Gibbs free energy – this time optimized for a particular protein | calculation of Gibbs free energy – this time optimized for a particular protein | ||
whose spatial coordinates were obtained through the previous Rosetta simulation. | whose spatial coordinates were obtained through the previous Rosetta simulation. | ||
Due to a lack of time, only a single structure per protein fragment was generated. | Due to a lack of time, only a single structure per protein fragment was generated. | ||
− | Results were analyzed in the same way (see | + | Results were analyzed in the same way (see figure 24) and further confirmed that |
the fragments were likely to exhibit higher affinity. Spatial coordinates of the | the fragments were likely to exhibit higher affinity. Spatial coordinates of the | ||
− | fragments were compared in PyMol and it turned out their position clashed, which | + | fragments were compared in PyMol (see figure 25) and it turned out their position clashed, which |
meant that real protein fragments did not have the option to stay in their optimal | meant that real protein fragments did not have the option to stay in their optimal | ||
position. This exposed another limitation of simulating interaction only of a single | position. This exposed another limitation of simulating interaction only of a single | ||
protein fragment at a time. Furthermore, there was still the option that fragments | protein fragment at a time. Furthermore, there was still the option that fragments | ||
could get digested by naturally present enzymes due to their clash and probable | could get digested by naturally present enzymes due to their clash and probable | ||
− | dysfunctionality | + | dysfunctionality, although that could compromise the whole G protein complex.<br> |
</div> | </div> | ||
<div style="display: flex; flex-direction: row; max-width: 95%; margin: 0 auto; justify-content: space-between;"> | <div style="display: flex; flex-direction: row; max-width: 95%; margin: 0 auto; justify-content: space-between;"> | ||
<img style="max-width: 45%; align-self: center;" src="https://static.igem.org/mediawiki/2020/9/99/T--UCopenhagen--m124.png"> | <img style="max-width: 45%; align-self: center;" src="https://static.igem.org/mediawiki/2020/9/99/T--UCopenhagen--m124.png"> | ||
+ | <figcaption> | ||
+ | <b>Figure 24: Change in Gibbs free energy of fragments produced by cleavage of G-alpha mutant m124. | ||
+ | The labels on the x-axis signify the position of each fragment in the original protein.</b> | ||
+ | </figcaption> | ||
+ | |||
<img style="max-width: 45%; align-self: center;" src="https://static.igem.org/mediawiki/2020/5/5d/T--UCopenhagen--m124-clash.png"> | <img style="max-width: 45%; align-self: center;" src="https://static.igem.org/mediawiki/2020/5/5d/T--UCopenhagen--m124-clash.png"> | ||
+ | <figcaption> | ||
+ | <b>Figure 25: An example of clashing positions of the fragments.</b> | ||
+ | </figcaption> | ||
</div> | </div> | ||
+ | |||
+ | <h2>Future Vision</h2> | ||
+ | |||
<div class="txt-btwn"> | <div class="txt-btwn"> | ||
In the future, simulating molecular dynamics via force fields could be applied | In the future, simulating molecular dynamics via force fields could be applied | ||
− | to obtain insights into interactions of protein fragments. Furthermore, this | + | to obtain insights into interactions of the protein fragments of G-alpha. Furthermore, this |
could inform us whether G-alpha cleavage lead to any conformational changes in | could inform us whether G-alpha cleavage lead to any conformational changes in | ||
G-beta, which could shed light on its ability or inability to recruit Ste5 | G-beta, which could shed light on its ability or inability to recruit Ste5 |
Revision as of 11:49, 26 October 2020
Introduction
Our dry lab efforts this year expanded to several different tracks, each
focused on different problems that needed solving by different methods and tools.
Our focus has primarily been on:
- ✧ Structural modeling and mutation of the alpha subunit of yeast G protein.
- ✧ Developing and comparing ODE models for our three iterative designs.
- ✧ Developing and comparing SDE models for our three iterative designs.
Differential Equations for Signaling Pathways
Since we were developing a biosensor which needed to sense changes in rather
low concentrations, it was crucial for us to have a good understanding of
the signal transduction pathway we were coupling it to, and what reactions
made up that pathway (please refer to our project design page
for more information on our biosensor design). Through our dry lab work,
we set out to answer these questions:
To answer these questions, we made use of ordinary differential equations (ODEs): Equations that illustrate how the state of a variable changes with time.
After working with ODEs for a while, we realized that the limits of using deterministic mathematics to model something as unpredictable as a cell. Therefore, we chose to utilize stochastic differential equations (SDEs).
Results from such a model contained everything from the ODE model and much more, thanks to expanding from single values to ranges of probabilities. With this tool, we started asking questions like – how much can the reporter concentration vary given the same biomarker concentration? How well does the biosensor function if the cells face senescence? How does the behavior of the cell change if one of the metabolites is toxic?
These are the kinds of biologically relevant questions that could be answered with stochastic modeling, partially because SDEs can represent inherent randomness present in all existing processes. Additionally, the SDEs allowed us to take our lack of knowledge of fine bioprocesses into account.
We needed to identify the values of our parameters of interest for the models to allow for realistic simulations to run. These values were found by researching literature, other models, and past iGEM wiki-pages. The parameters were selected to be consistent among our different designs. The parameters found were taken as rough estimates (see figure 1), because the performance of biological systems can vary for several reasons that have yet to be elucidated. However, rough estimates are enough to get great insights into the mechanics of cellular pathways. Optimally, producing more precise laboratory data, and fitting the models hereto, would be the next step to obtain more realistic modeling results.
- 1) How do we make sure that our biosensor can react to the low concentrations of cytokines found in sweat, within a reasonable time limit set by our Human Practices research?
- 2) Where are the bottlenecks in the pathway, and can we fix them?
- 3) Which signal transduction steps are robust, and which are sensitive?
To answer these questions, we made use of ordinary differential equations (ODEs): Equations that illustrate how the state of a variable changes with time.
After working with ODEs for a while, we realized that the limits of using deterministic mathematics to model something as unpredictable as a cell. Therefore, we chose to utilize stochastic differential equations (SDEs).
Results from such a model contained everything from the ODE model and much more, thanks to expanding from single values to ranges of probabilities. With this tool, we started asking questions like – how much can the reporter concentration vary given the same biomarker concentration? How well does the biosensor function if the cells face senescence? How does the behavior of the cell change if one of the metabolites is toxic?
These are the kinds of biologically relevant questions that could be answered with stochastic modeling, partially because SDEs can represent inherent randomness present in all existing processes. Additionally, the SDEs allowed us to take our lack of knowledge of fine bioprocesses into account.
How We Model
We generated two classes of models based on ODEs and SDEs, respectively. First, we'll walk through the details of ODE implementation, after which we built upon them by introducing features unique to SDEs. Some hints about alternative modeling frameworks were added at the bottom of this page for future iGEMers!We needed to identify the values of our parameters of interest for the models to allow for realistic simulations to run. These values were found by researching literature, other models, and past iGEM wiki-pages. The parameters were selected to be consistent among our different designs. The parameters found were taken as rough estimates (see figure 1), because the performance of biological systems can vary for several reasons that have yet to be elucidated. However, rough estimates are enough to get great insights into the mechanics of cellular pathways. Optimally, producing more precise laboratory data, and fitting the models hereto, would be the next step to obtain more realistic modeling results.
Ordinary Differential Equations (ODEs)
To build our ODEs, the below-stated three steps were followed:- 1) Formulate a diagram that specifies the key players (state variables) and describes all possible ways these variables might interact with each other.
- 2) Explicitly write out the quantitative form of all interactions.
- 3) Convert the explicit quantitative interactions into a mathematical framework.
Step 1 - Making a Diagram
In order to make a diagram, we needed to have an idea of what our design should look like. For this, one could assume that:- a) The presence of a biomolecule in proximity to the cellular membrane lead to the binding of our receptor of interest to its ligand/biomolecule, and then the binding of a coreceptor to this complex.
- b) Our receptor and coreceptor were both transmembrane proteins, with intracellular domains that associate intracellularly following the extracellular association.
- c) Their intracellular association was facilitated by the natural affinity of the two halves of a split-ubiquitin molecule (for reference, see design page). This association and reconstitution of a full ubiquitin prompted its cleavage by naturally occurring deubiquitinases.
- d) Such a cleaving event released a transcription factor that was bound to the C-terminal half of the split ubiquitin.
- e) The freed transcription factor traveled to the nucleus and initiated transcription of the gene into mRNA.
- f) mRNA was translated into a reporter protein.
Step 2 - Writing Interactions in a Quantitative Form
The next step was to find out how much product was produced from a given amount of reactant. As we were dealing with signaling pathways, our concentrations were changing in accordance with the association and dissociation of the variables. The equation derived hereof can be found below:k1[A][B]=>[C]
where k1 is a kinetic parameter, [A] and [B] are reactant concentrations, and [C] is a product concentration. Note that such a reaction proceeds only in one direction. However, the reality of many reactions was better represented by breakdown into a forward reaction and a backward reaction. This was because it was empirically observed that chemical reactions tend to reach a state of equilibrium with nonzero concentrations of reactants and products. This was understood as a situation described by:k_on[A][B]<=>k_off[C]
where k_on and k_off are kinetic parameters (binding and unbinding constant, respectively). Note that the reaction proceeds in both directions. After some time, equilibrium is reached, which is described by:k=k_on/k_off=[C]/([A][B]).
In practice, it was much easier to deduce k compared to its constituents k_on and k_off, because k could be calculated in a straightforward manner from observed concentrations. This was first addressed later. The equilibrium was now written as:
k_on * [biomolecule] * [receptor] * [coreceptor] <=> k_off * [complex]
However, the next reaction where deubiquitinase cleaved the reassociated split-ubiquitin proceeded only in a forward direction, as it is irreversible – once the transcription factor becomes released, it is never again part of the transmembrane protein.k_on * [complex] * [deubiquitinase] => [transcription factor]
The next step was transcription, which was a slightly more complicated reaction, only proceeding in the forward direction. The kinetics of transcription follows a simple form of Hill equation that was written as [A]/(Kd+[A])=>[B] where [A] and [B] are concentrations and Kd is a dissociation constant (a measure of how tightly the transcription factor bound a promoter). This equation further needed to be multiplied by a specific rate of transcription r=g*R/L where g is a number of copies of a given gene, R is a general rate of transcription and L is a gene length.rate * [transcription factor] / (Kd + [transcription factor]) => [mRNA]
Translation rate was a simpler case, as knowing the concentration of mRNA and the specific translation rate was enough to determine the translation rate. The specific translation rate is r=R/L where R is a general rate of translation and L is the length of mRNA.rate * [mRNA] => [protein]
Step 3 - Converting quantified interactions into a mathematical framework.
Concerting quantified interactions into a mathematical framework is the first step that was model dependent. Here, we wanted to create a model based on ordinary differential equations. Such equations described how values of given variables change over time, or rather how they changed during each time step. The ODE-based models worked by computing a set of equations iteratively – i.e. output of a previous computation was used as an input for a next run. Turning quantified interactions into differential equations was done by rewriting them from a standpoint of variables. Each differential equation needed to capture all changes that happen to a single variable. This included both additions and subtractions – if a certain molecule was a product of one reaction but a reactant in other reaction, the concentration of the molecule could both rise or drop.In the following reaction, the first term captured emergence of the complex by association of the biomarker, the receptor and the coreceptor. The second term captured the backward reaction where the complex dissociated into its constituent parts. The third term captured spending of the complex by catalytic activity of deubiquitinase.
D(complex) = k_on * [biomolecule] * [receptor] * [coreceptor] – k_off * [complex] – k_cat * [complex] * [deubiquitinase]
We had now made a system of ordinary differential equations that described kinetics of our signaling pathway. Next, we needed to translate our mathematical model into a programming language of choice, e.g. Matlab. Each variable was named ‘c’ and an arbitrary number. Parameters and whole equations were handled in a similar fashion:d4 = k1 * c1 * c2 * c3 – k2 * c4 – k3* c4 * k4;
Except from providing initial conditions, we did not need to worry about all the c’s and d’s any longer. The k’s, however, informed about the speed and order by which our system worked.As scientific work must be translatable at all times, we paid special attention to satisfying standards like MIRIAM. Knowing that many various software tools were used for pathway modeling, we utilized Moccasin to turn our Matlab code into SBML files that could be opened with a wide variety of programs. Below are snippets of code from our models:
Click for our parameters, plot, and differential equations
Stochastic Differential Equations (SDEs)
This part built upon the previous section concerning ODEs. In fact, the ODEs described in the previous section were so similar to the SDEs that we took our implemented system of differential equations and copied it in the right place of our code for the SDE-based model. Before that, we simplified the equation of our mathematical model to the form where M denoted all the deterministic parts. In other words, M was all the contribution coming from the signaling pathway of our interest.D(complex) = k_on * [biomolecule] * [receptor] * [coreceptor] – k_off * [complex] – k_cat * [complex] * [deubiquitinase] = M
Now, we needed a term that was driven entirely by randomness. We called it N. If we decided that N stood for a stochastic part of the equation, we got an equation like this:
D(complex) = M + N
In the context of SDEs, N represents a Wiener process. N is seen as a variable, whose value changes randomly but without any discontinuity occuring in its graph.
This allowed us to model the effects of randomness and probabilities. Going back to the mathematical formula, this was expanded further, so that we obtained a stochastic differential equation that we later implemented in Matlab.
D(complex) = M + N = k1 * [biomolecule] * [receptor] * [coreceptor] - k2 * [complex] - k3 * [complex] * [deubiquitinase] + N
We implemented our SDE-based model in Matlab. Yet, there was a limited built-in support for SDEs in Matlab, so we used a package called SDETools that was available online. It separated the deterministic and stochastic parts of each equation so that the deterministic parts stayed grouped together, which was why we could easily copy and paste our system of ODEs. The stochastic part was then defined separately.This was no issue because Wiener processes can be defined in any number of dimensions, i.e. if we had ten differential equations, we simply defined N as a 10-dimensional process – one dimension per equation.
As mentioned earlier, our intention was to make the code readable and translatable, so we followed the same sets of standards for our SDE models as for our ODE models. One exception was that no SBML file was generated as there was no software support for SDETools due to it being an unofficial package. For the structure of our Matlab code, click here for contact information, here for parameters, here for deterministic equations, Wiener process, vectors and options, and finally, running the code.
To Future iGEMers!
There are some limitations to SDETools. Using a programming
language called Julia together with available packages could give a better SDE
modeling experience (our recommendation for future teams).
There are more ways to model that are quite different from differential equations. Some of them are rules-based or stochastic, such as: Gillespie algorithm, Markov chains, Petri nets.
Different algorithms come with different limitations or assumptions (that may not be immediately visible) so it could be a good idea to model the same system in several ways to make sure that one's predictions are meaningful.
Figure 5: Comparison of amplification provided by a one-to-one relationship versus one-to-many relationship.
Note On Alternative Modeling Paradigms
Apart from ODEs and SDEs, there were many other modeling approaches that are useful for modeling of signaling pathways. The first one to mention is partial differential equations (PDEs). They are very similar to ODEs, but they also allow a variable change to be dependent on other variables, and not only on time. This extension was essential for modeling physical processes that are located in space (e.g. diffusion).There are more ways to model that are quite different from differential equations. Some of them are rules-based or stochastic, such as: Gillespie algorithm, Markov chains, Petri nets.
Different algorithms come with different limitations or assumptions (that may not be immediately visible) so it could be a good idea to model the same system in several ways to make sure that one's predictions are meaningful.
Pheromone Cascade – Attractive but Complicated
Taking inspiration from last year’s UCopenhagen iGEM team, Ovulaid, we modeled the
yeast pheromone pathway in a similar way to what the did (using only two equations to
abstract the whole cascade) and provided it with a parameter fixed for
experimental data found in literature. After some experimentation with our
models, we became convinced that the reactions, in order to lead to empirically observed hypersensitivity, were necessarily
more complex compared to what some papers claimed.
Notably, one-to-one reactions could not provide significant amplification from a theorical standpoint. After reading through literature, we found more detailed descriptions of the pheromone cascade which included a series of complexes and one-to-many interactions. In other words, a lot of amplification. As such, we became convinced that our model now more accurately represented the complexity found in the pheromone pathway.
Notably, one-to-one reactions could not provide significant amplification from a theorical standpoint. After reading through literature, we found more detailed descriptions of the pheromone cascade which included a series of complexes and one-to-many interactions. In other words, a lot of amplification. As such, we became convinced that our model now more accurately represented the complexity found in the pheromone pathway.
We found a brilliant curated model on biomodels.com. We extracted the important
part of the cascade and added it in our model. This complicated our model, however, making the SDE implementation run into RAM memory issues.
Impact on Wet Lab - Experiments
How our modeling was tied to wet lab experiments
Since the beginning of our iGEM journey, wet lab and dry lab worked closely together. Thanks to our supervisors, we always intended to use our dry lab activities to guide our wet lab decisions. As our project ended up focusing a great deal on protein analysis, the wet lab work with G-alpha would not have been possible without the simultaneous work in the dry lab.The dry lab guided the wet lab in finding out which of the several biosensor designs was the most viable, in shaping the design of the assays for our Saccharomyces cerevisiae cells, and also in engineering the proteins with the desired properties.
In all the following experiments, all values were held equal except from the one studied at a time.
Our ODE Models
Dynamic Ranges
Dynamic ranges showed how well different designs responded to a range of
interleukin concentrations. It was observed that the lack of pheromone
cascade in the design #1 and the design #2 lead to insufficient sensitivity.
The pheromone cascade increased sensitivity by several orders of magnitude.
These findings resulted in more effort allocated to the design #3 (with G-alpha and TEV protease).
Sensitivity Analysis
Sensitivity analysis showed how perturbations of parameter values impacted
the reporter concentration. It identified more robust parameters (i.e. resilient
to perturbations) and more sensitive parameters (i.e. even slight perturbations
lead to a noticeable impact). Additionally, it was observed how the importance
of certain parameters changed over time. In the design #1, the affinity of a
transcription factor to a promoter (k15 in the figure) was a sensitive parameter,
so that the reporter concentration range could be adjusted by choice of a promoter.
The reporter concentration was shown to be sensitive to interleukin concentration (c_1 in the figure).
Correlation Among Pathway Proteins
Here we performed principal component analysis (PCA) and obtained a matrix of all correlations between variables in the design #3 (spectrum from blue to red represented perfect negative correlation to perfect positive correlation, respectively). The results allowed us to understand indirect interactions between proteins in the signaling pathway. Provided that laboratory tools and methods for advanced pathway engineering existed, this figure could guide us in pathway optimization.Impact of Promoter
Affinity of a transcription factor (TF) to a promoter was varied in
this case study. High-affinity TF-promoter pair (Kd = 1 nM) lead to an increased
sensitivity accompanied by significantly higher reporter concentrations,
compared with a low-affinity TF-promoter pair (Kd = 100 nM). The biosensor’s dynamic range
as well as reporter concentrations can be adjusted by the right choice of
a TF-promoter pair.
Inflammatory vs baseline reporter concentrations
In this case study, the design #3 was subjected to low and high
interleukin concentration (i.e. levels of inflammation).
High inflammation (c = 10 nM) lead to a fast development of high reporter
concentrations, compared with low inflammation (c = 1 pM). Significant differences
of the plotted curves suggested the design #3 could distinguish between
low and high inflammation, as long as the measurements were done at the
right time (preferably ~5 hours mark).
This fits perfectly with our Human Practices considerations, where multiple experts preferred a design that provided fast results, and patients were willing to wear the patch for up to 24 hours, according to our survey results!
This fits perfectly with our Human Practices considerations, where multiple experts preferred a design that provided fast results, and patients were willing to wear the patch for up to 24 hours, according to our survey results!
Impact of proteolytic activity
Proteolytic activity of TEV protease shifted the dynamic range of
design #3 by orders of magnitude, without a change in overall reporter
concentrations. There were many versions of TEV proteases available
with various levels of proteolytic activity. Protease choice had a
significant impact on the dynamic range of a biosensor. See figure 12 for comparison of high proteolytic actvity (k_cat = 0.18 s^-1) and low proteolytic actvity (k_cat = 0.000907 s^-1).
SDEs - hypothetical scenarios with effects of noise
Figures are illustrative. Use of SDE models was limited by computational requirements of the current implementation. Thus, desired phenomena were amplified to be observable on a short timescale. SDE models were an extension that was less of a design tool and more of an analytic tool for further improvement of the biosensor – it aided with identification of unknown behaviors.Rapid aging lead to disorder in the cell
This case study examined the impact of noise increasing with time. Such a phenomenon could be present in senescent cells, where by-products of aging obstruct the reliability of cellular processes. If such divergence of curves was observed in our sweat patches, we would have known whether the dysfunctionality is primarily time dependent.Longer measurement allows the signal pathway to yield more consistent results
This case study examined the impact of noise decreasing with time. Such a phenomenon could be observed in any cell, especially soon after triggering the receptor system before the pathway processes overcame the magnitude of the biochemical background. If such divergence of curves was observed in our sweat patches, we would know that this is a perfectly expected behavior coming from the physical law that longer measurements yield more precise results.Unknown process decreases and increases disorder in the cell
This case study examined a rather exotic-looking hypothetical scenario where high levels of noise diminished and came back to initial levels afterwards. Such bizarre divergence of curves could indicate that the cellular pathway only worked under very specific conditions that were achieved only for a short time.Natural variability in measurements
This case study examined natural variability of measurements that followed the normal distribution. By fitting data from experiments to the divergence of these curves, we could know how much noise was present in our biosensor or provided by the measurement technique.Receptors are not specific enough
This case study examined impact of noisy activation of signaling pathways that could be mediated by a receptor system that was not specific enough. If such divergence of curves was observed in our sweat patches, we would have known the problem was possibly caused by unspecific receptors.Receptor disrupts cellular processes
This case study examined the impact of noise increase tied to reporter concentration. If such divergence of curves was observed in our sweat patches, we would have known that the reporter was causing significant harm to the reliability of cellular processes.Variability in Interleukin Levels
This case study examined the impact of variability in interleukin concentration. This behavior could be expected in our sweat patch as the baseline interleukin concentrations likely varies among people. If such divergence of curves was observed in our sweat patches, we would have known that this behavior was expected and we could deduce variability among people.Diminishing disorder of cellular processes
This case study examined the effect of initially high but diminishing noise. This exotic scenario could represent the cellular process stabilizing at random states. Such a behavior is not expected in a functioning biosensor and could indicate there was a serious issue with part of the signaling pathway, or perhaps that the cells had been compromised.G-alpha and Structural Modelling with Rosetta
Just to refresh, our work on G-alpha hinged on finding a mutant that preserved
natural function when there was no signal (i.e. no cytokine), and also dissociated
appropriately in the presence of a signal.
Following the recommendation from team Sinisens (iGEM Aalto-Helsinki 2020), we looked into Rosetta Software Suite. Alanine scan was a Rosetta script available via an old version of Robetta server, and we used it to calculate how Gibbs free energy of the two proteins (G-alpha and G-beta) changed upon changing a single residue to alanine. Changes in Gibbs free energy correlate with binding affinity.
We submitted a list of residues selected according to our previous analyses. A crystal structure of the bovine G protein complex served as a template for spatial structural positioning of our proteins.
Following the recommendation from team Sinisens (iGEM Aalto-Helsinki 2020), we looked into Rosetta Software Suite. Alanine scan was a Rosetta script available via an old version of Robetta server, and we used it to calculate how Gibbs free energy of the two proteins (G-alpha and G-beta) changed upon changing a single residue to alanine. Changes in Gibbs free energy correlate with binding affinity.
We submitted a list of residues selected according to our previous analyses. A crystal structure of the bovine G protein complex served as a template for spatial structural positioning of our proteins.
Thus, we needed to simulate binding affinities (or change in Gibbs free energy – delta delta G) of seven G-alpha proteins and nine fragments formed after cleavage. We eventually landed on using a Rosetta script called flex_ddg, as this script has been shown to be less computationally demanding than other scripts, while making great predictions of molecular dynamics .
We particularly found the script that was written with a capability to include refinement of every protein structure many times to be useful. We followed recommendations from script creators and went with 35 000 independent refinements per structure. Although the script was initially developed for mutagenesis simulations, we found an easy way to adjust it for our own use. For each of our 16 proteins (mutants + fragments), we generated 35 independently optimized structures to probe the conformational space (this took 20 days on 16 CPUs with 64 GBs of RAM). We used the scoring function REF2015 as it was shown to evaluate results of flex_ddg better than others (https://pubs.acs.org/doi/10.1021/acscentsci.8b00717).
Results of the simulation were surprising as they revealed our
reasoning might have been flawed. Essentially, our assumption had been that
cleaved fragments dissociate from the protein, yet our results, as can be observed in figures 22 and 23,
showed that the majority of protein fragments had a higher affinity than the parent proteins.
We also realized that our simulation was not directly answering our true question
– could cleavage of our mutant G-alpha trigger Ste5 recruitment by G-beta?
Directly simulating that could include up to seven simultaneously interacting
proteins, which translated into a very demanding computational task.
After analyzing the results, we decided to run another simulation with the same script but different input, to focus on fragments of the best-scoring protein – M124 (GPA mutant with TEV sites numbered 1, 2, and 4). What made this simulation different and more precise that the previous one was a ground reference for the calculation of Gibbs free energy – this time optimized for a particular protein whose spatial coordinates were obtained through the previous Rosetta simulation. Due to a lack of time, only a single structure per protein fragment was generated. Results were analyzed in the same way (see figure 24) and further confirmed that the fragments were likely to exhibit higher affinity. Spatial coordinates of the fragments were compared in PyMol (see figure 25) and it turned out their position clashed, which meant that real protein fragments did not have the option to stay in their optimal position. This exposed another limitation of simulating interaction only of a single protein fragment at a time. Furthermore, there was still the option that fragments could get digested by naturally present enzymes due to their clash and probable dysfunctionality, although that could compromise the whole G protein complex.
After analyzing the results, we decided to run another simulation with the same script but different input, to focus on fragments of the best-scoring protein – M124 (GPA mutant with TEV sites numbered 1, 2, and 4). What made this simulation different and more precise that the previous one was a ground reference for the calculation of Gibbs free energy – this time optimized for a particular protein whose spatial coordinates were obtained through the previous Rosetta simulation. Due to a lack of time, only a single structure per protein fragment was generated. Results were analyzed in the same way (see figure 24) and further confirmed that the fragments were likely to exhibit higher affinity. Spatial coordinates of the fragments were compared in PyMol (see figure 25) and it turned out their position clashed, which meant that real protein fragments did not have the option to stay in their optimal position. This exposed another limitation of simulating interaction only of a single protein fragment at a time. Furthermore, there was still the option that fragments could get digested by naturally present enzymes due to their clash and probable dysfunctionality, although that could compromise the whole G protein complex.
Future Vision
In the future, simulating molecular dynamics via force fields could be applied
to obtain insights into interactions of the protein fragments of G-alpha. Furthermore, this
could inform us whether G-alpha cleavage lead to any conformational changes in
G-beta, which could shed light on its ability or inability to recruit Ste5
and thus trigger the pheromone cascade.
- https://pubmed.ncbi.nlm.nih.gov/29401388/ ↩