Line 291: | Line 291: | ||
<div class="econ_box"> | <div class="econ_box"> | ||
<p class="para">In addition, the project would be moot without a consideration of the insulin market as a whole. </p> <br> <p class="para">Modelling helped us to gain insight into the global insulin market, which informed our approach towards entrepeneurship. | <p class="para">In addition, the project would be moot without a consideration of the insulin market as a whole. </p> <br> <p class="para">Modelling helped us to gain insight into the global insulin market, which informed our approach towards entrepeneurship. | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> |
</div> | </div> | ||
</div> | </div> | ||
Line 328: | Line 328: | ||
<div class="cell_model"> | <div class="cell_model"> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> The notation of Weisse et al has been adopted as follows: </p> |
<table style="width:70%;margin:auto;"> | <table style="width:70%;margin:auto;"> | ||
<caption> Notation for biochemical species' considered in Weisse et al. model</caption> | <caption> Notation for biochemical species' considered in Weisse et al. model</caption> | ||
Line 435: | Line 435: | ||
<div class="cyto_model"> | <div class="cyto_model"> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;">Our first expression system was cytoplasmic expression in E. coli. We modelled the rate of change of five biochemical species in the cell: |
</p> | </p> | ||
<table style="width:70%;margin:auto;"> | <table style="width:70%;margin:auto;"> | ||
Line 466: | Line 466: | ||
<br> | <br> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> A diagram showing the species we modelled and notation used is shown below:</p> |
<img class='crispy center-block' src="https://static.igem.org/mediawiki/2017/6/6a/T--Sydney_Australia--cyto_model_schematic.png" style="width:80%; margin:auto;"/> | <img class='crispy center-block' src="https://static.igem.org/mediawiki/2017/6/6a/T--Sydney_Australia--cyto_model_schematic.png" style="width:80%; margin:auto;"/> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> The following reactions were considered </p> |
<table style="width:70%;margin:auto;"> | <table style="width:70%;margin:auto;"> | ||
<caption> List of reactions</caption> | <caption> List of reactions</caption> | ||
Line 519: | Line 519: | ||
</tr> | </tr> | ||
</table> | </table> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> Here, \(\omega_p(a)\), the rate of transcription, is an energy dependent process.</p> |
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;">We used the transcription rate form used in Weisse et al to denote the amount being transcribed (\(\omega_p(a)\)). That is,</p> |
\[\color{#3e3f3f}{\omega_p(a)=w_p \frac{a}{\theta_p+a} } \] | \[\color{#3e3f3f}{\omega_p(a)=w_p \frac{a}{\theta_p+a} } \] | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;">Where \(w_p\) is the maximal rate of transcription, dependent on the speed of transcriptional elongation, as well as the gene length, induction and copy number. \(a\) is the energy in the cell such as ATP (transcription is an energy dependent process), and |
\(\theta_p\) is the transcriptional threshold of the recombinant protein.</p> | \(\theta_p\) is the transcriptional threshold of the recombinant protein.</p> | ||
<br> | <br> | ||
<br> | <br> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;">In addition, we used the form in Weisse et al. for the translation rate term </p> |
\[\color{#3e3f3f}{\upsilon_p(c_p,a)=c_p \frac{\gamma(a)}{n_p} } \] | \[\color{#3e3f3f}{\upsilon_p(c_p,a)=c_p \frac{\gamma(a)}{n_p} } \] | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> Where \(n_p\) is the length of recombinant protein, and \(\gamma(a)\) is an expression for the rate of transcriptional elongation: </p> |
\[\color{#3e3f3f}{\gamma(a)=\frac{\gamma_{max} a}{K_{\gamma} + a} } \] | \[\color{#3e3f3f}{\gamma(a)=\frac{\gamma_{max} a}{K_{\gamma} + a} } \] | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> Where \(\gamma_{max}\) is the maximal rate of translation, \(K_{\gamma}\) is the translational elongation threshold, and \(a\) is the energy in the cell.</p> |
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;">For the model of inclusion body aggregation, we assumed first order deposition of monomers of unfolded protein, dependent on the concentration of unfolded protein. as in Hoffmann et al (2001). |
</p> | </p> | ||
<br> | <br> | ||
Line 537: | Line 537: | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> Using the law of mass action kinetics we can derive a set of ordinary differential equations from these reactions.</p> |
</p> | </p> | ||
<h1 style="font-family:'Quicksand';"> Summary of Cytoplasmic Expression Model</h1> | <h1 style="font-family:'Quicksand';"> Summary of Cytoplasmic Expression Model</h1> | ||
Line 547: | Line 547: | ||
<h1 style="font-family:'Quicksand';"> Parametrising the model </h1> | <h1 style="font-family:'Quicksand';"> Parametrising the model </h1> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;">The following table shows the parameters we needed to find for our model, and the values we used</p> |
<table style="width:90%;margin:auto;"> | <table style="width:90%;margin:auto;"> | ||
<caption> Cytoplasmic Expression Model Parameters. * Set to 0 as degradation is dominated by the rate of dilution due to cell division for stable proteins [3]</caption> | <caption> Cytoplasmic Expression Model Parameters. * Set to 0 as degradation is dominated by the rate of dilution due to cell division for stable proteins [3]</caption> | ||
Line 644: | Line 644: | ||
<div class="peri_model"> | <div class="peri_model"> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> Next we looked at periplasmic expression of our recombinant protein in E. coli. We modelled the rate of change of 6 species </p> |
<table style="width:70%;margin:auto;"> | <table style="width:70%;margin:auto;"> | ||
<caption> Cytoplasmic Expression Model Variables</caption> | <caption> Cytoplasmic Expression Model Variables</caption> | ||
Line 677: | Line 677: | ||
</table> | </table> | ||
<br> | <br> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> A diagram showing the species we modelled and notation used is shown below:</p> |
<img class='crispy center-block' src="https://static.igem.org/mediawiki/2017/3/3f/T--Sydney_Australia--peri_model_schematic.png" style="width:80%; margin:auto; padding:5px;"/> | <img class='crispy center-block' src="https://static.igem.org/mediawiki/2017/3/3f/T--Sydney_Australia--peri_model_schematic.png" style="width:80%; margin:auto; padding:5px;"/> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> The following reactions were considered </p> |
<table style="width:70%;margin:auto;"> | <table style="width:70%;margin:auto;"> | ||
<caption> List of reactions</caption> | <caption> List of reactions</caption> | ||
Line 736: | Line 736: | ||
</table> | </table> | ||
<br> | <br> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> here, \(\omega(a)\) and \(\upsilon(c_p,a)\) are as in the cytoplasmic reactions. |
The amount being transported is found with the term \(\tau_p(p_t,a)\). | The amount being transported is found with the term \(\tau_p(p_t,a)\). | ||
Protein translocation to the periplasm occurs via an ATP-dependent motor protein, secA [4]. | Protein translocation to the periplasm occurs via an ATP-dependent motor protein, secA [4]. | ||
Line 747: | Line 747: | ||
</figure> | </figure> | ||
<br> | <br> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> Following the logic used to derive the translation rate in [1], we derive the net rate of translocating a protein \(p\) by defining \(K_p:=\frac{k_1k_2}{k_{-1}+k_2}\). This leads to </p> |
\[\color{#3e3f3f}{\tau_p(p_t,a)=p_t\Big(\frac{n_p}{50}\Big(\frac{1}{K_pa}+\frac{1}{k_2} \Big)+\frac{1}{k_t}\Big)^{-1}}\] | \[\color{#3e3f3f}{\tau_p(p_t,a)=p_t\Big(\frac{n_p}{50}\Big(\frac{1}{K_pa}+\frac{1}{k_2} \Big)+\frac{1}{k_t}\Big)^{-1}}\] | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;">If we assume the final termination step is fast, so \(\frac{1}{k_t}<< \frac{n_p}{50}\Big(\frac{1}{K_pa}+\frac{1}{k_2} \Big) \), this is approximately equal to |
</p> | </p> | ||
\[\color{#3e3f3f}{\tau_p(p_t,a)\approx 50p_t \frac{\epsilon(a)}{n_p}\qquad \epsilon(a):=\frac{\epsilon_{max}a}{K_{\epsilon}+a} }\] | \[\color{#3e3f3f}{\tau_p(p_t,a)\approx 50p_t \frac{\epsilon(a)}{n_p}\qquad \epsilon(a):=\frac{\epsilon_{max}a}{K_{\epsilon}+a} }\] | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;">Where \(\epsilon_{max}\) is the maximal translocation rate, \(K_{\epsilon}\) is the threshold, and \(n_p\) is the length of the protein in amino acids |
</p> | </p> | ||
<br> | <br> | ||
<h3> Parametrising Translocation</h3> | <h3> Parametrising Translocation</h3> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> To find the parameters for translocation (\(\epsilon(a)\)) and (\(K_{\epsilon}\)), we used kinetic parameters determined in [5]. They measured translocation of a 346aa protein proOmpA and found the apparent Km of SecA was 50nM, and the threshold was 2.7 proOmpa/site/min. This converts to 2.7 \(\cdot\) 346 proOmpA/site/min aa/proOmpa \(\rightarrow\) 934.2 aa/molec/min |
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> Using the law of mass action kinetics we can derive a set of ordinary differential equations from these reactions.</p> |
</p> | </p> | ||
Line 771: | Line 771: | ||
<h1 style="font-family:'Quicksand';"> Parametrising the model </h1> | <h1 style="font-family:'Quicksand';"> Parametrising the model </h1> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;">The following table shows the parameters we needed to find for our model, and the values we used</p> |
<table style="width:90%;margin:auto;"> | <table style="width:90%;margin:auto;"> | ||
<caption> Periplasmic Expression Model Parameters. † Doubled relative to cytoplasmic folding rate to reflect the effect of an oxidising environment on disulfide bond formation. * Set to 0 as degradation is dominated by the rate of dilution due to cell division for stable proteins [3]</caption> | <caption> Periplasmic Expression Model Parameters. † Doubled relative to cytoplasmic folding rate to reflect the effect of an oxidising environment on disulfide bond formation. * Set to 0 as degradation is dominated by the rate of dilution due to cell division for stable proteins [3]</caption> | ||
Line 994: | Line 994: | ||
<div class="bacillus_model"> | <div class="bacillus_model"> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> We also developed a model of our secretory protein expression system in bacillus subtilis. The model included 6 species </p> |
<table style="width:70%;margin:auto;"> | <table style="width:70%;margin:auto;"> | ||
<caption> Cytoplasmic Expression Model Variables</caption> | <caption> Cytoplasmic Expression Model Variables</caption> | ||
Line 1,027: | Line 1,027: | ||
</table> | </table> | ||
<br> | <br> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> A diagram showing the species we modelled and notation used is shown below:</p> |
<img class='crispy center-block' src="https://static.igem.org/mediawiki/2017/b/be/T--Sydney_Australia--sec_model_schematic.png" style="width:80%; margin:auto; padding:5px;"/> | <img class='crispy center-block' src="https://static.igem.org/mediawiki/2017/b/be/T--Sydney_Australia--sec_model_schematic.png" style="width:80%; margin:auto; padding:5px;"/> | ||
<br> | <br> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> Structurally, this is the same process as the periplasmic expression system, so the equations' structure is the same. However the parameters are different, reflecting the different environment of bacillus and medium and its effect on expression of recombinant protein </p> |
<h1 style="font-family:'Quicksand';"> Summary of Bacillus Secretory Expression Model </h1> | <h1 style="font-family:'Quicksand';"> Summary of Bacillus Secretory Expression Model </h1> | ||
Line 1,042: | Line 1,042: | ||
\[\color{#3e3f3f}{\frac{d}{dt}{p}_f=k_fp_u-(k_d+\lambda)p_f}\] | \[\color{#3e3f3f}{\frac{d}{dt}{p}_f=k_fp_u-(k_d+\lambda)p_f}\] | ||
− | We were unfortunately unable to parametrise the bacillus model, so for our in silico experiments we focused on comparing cytoplasmic and periplasmic E. coli expression. | + | <p>We were unfortunately unable to parametrise the bacillus model, so for our in silico experiments we focused on comparing cytoplasmic and periplasmic E. coli expression.</p> |
</div> | </div> | ||
Line 1,049: | Line 1,049: | ||
<div class="insilico_exp"> | <div class="insilico_exp"> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> Once we had modelled our different expression systems for recombinant insulin, we integrated them into the whole cell model developed in [1]. The methodology behind integrating models of our expression system into this model was to more accurately reflect reality. Recombinant protein expression occurs within a complex cellular environment with finite resources. A model which ignores the actiities of the host cells would ignore important host-circuit interactions. Ignoring the finite resources of the cell may skew our prediction of the yield of our expression systems. |
</p> | </p> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> We then interrogated these models for insights into how to optimise the expression of insulin </p> |
<h1 style="font-family:'Quicksand';"> Comparing Cytoplasmic and Periplasmic Expression </h1> | <h1 style="font-family:'Quicksand';"> Comparing Cytoplasmic and Periplasmic Expression </h1> | ||
Line 1,063: | Line 1,063: | ||
<div class="row physiological_modelling"> | <div class="row physiological_modelling"> | ||
<div class="phys_box" style="width:80%;margin:auto;"> | <div class="phys_box" style="width:80%;margin:auto;"> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;">For our physiological modelling, we used a model of subcutaneous insulin absorption developed in [1] and used it to develop a model relating the free energy of insulin hexamer formation and insulin dynamics (the time of peak of action, and the duration of action). </p> |
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> |
The authors of [1] developed a system of partial differential equations to describe the insulin infusion process. They modelled the change in three species: | The authors of [1] developed a system of partial differential equations to describe the insulin infusion process. They modelled the change in three species: | ||
</p> | </p> | ||
Line 1,086: | Line 1,086: | ||
</tr> | </tr> | ||
</table> | </table> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> They modelled the conversion between hexameric and dimeric insulin as follows </p> |
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> Insulin<sub>Hexamer </sub> \(\color{#3e3f3f}{\rightleftharpoons}\) Insulin<sub>Dimer</sub></p> |
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;">Where the forward rate was called \(\color{#3e3f3f}{P}\) and the reverse rate was \(\color{#3e3f3f}{PQ}\) where we can interpret \(\color{#3e3f3f}{P}\) as the production rate and \(\color{#3e3f3f}{Q}\) as the equilibrium constant.</p> |
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> The final model was as follows </p> |
\[\color{#3e3f3f}{\eqalignno{{\partial c_{d}(t,r)\over\partial t}=&\,P\left(c_{h}(t,r)-Qc_{d}(t,r)^{3}\right)-B_{d}c_{d}(t,r)\cr&+D\nabla^{2}c_{d}(t,r),\cr{\partial c_{h}(t,r)\over\partial t}=&\,-P\left(c_{h}(t,r)-Qc_{d}(t,r)^{3}\right)\cr&+D\nabla^{2}c_{h}(t,r)}}\] | \[\color{#3e3f3f}{\eqalignno{{\partial c_{d}(t,r)\over\partial t}=&\,P\left(c_{h}(t,r)-Qc_{d}(t,r)^{3}\right)-B_{d}c_{d}(t,r)\cr&+D\nabla^{2}c_{d}(t,r),\cr{\partial c_{h}(t,r)\over\partial t}=&\,-P\left(c_{h}(t,r)-Qc_{d}(t,r)^{3}\right)\cr&+D\nabla^{2}c_{h}(t,r)}}\] | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> Where \(\color{#3e3f3f}{P, Q, B_d, D}\) are parameters, and exogenous insulin flow is obtained by integrating the expression denoting the amount of insulin dimer entering the bloodstream: |
\[\color{#3e3f3f}{I_{ex}(t)=B_{d}\int\limits_{V_{sc}}c_{d}(t,r)dV.}\] | \[\color{#3e3f3f}{I_{ex}(t)=B_{d}\int\limits_{V_{sc}}c_{d}(t,r)dV.}\] | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> The parameters found for the different insulin analogues and their impact on insulin dynamics are </p> |
<table style="width:70%;margin:auto;"> | <table style="width:70%;margin:auto;"> | ||
<caption> Table 2. Parameter Values and resultant Dynamics for different insulin analogues. Values of parameters from [1] Table IV. Insulin dynamics taken from [1] Fig. 6</caption> | <caption> Table 2. Parameter Values and resultant Dynamics for different insulin analogues. Values of parameters from [1] Table IV. Insulin dynamics taken from [1] Fig. 6</caption> | ||
Line 1,141: | Line 1,141: | ||
</table> | </table> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;">Since the parameter \(\color{#3e3f3f}{Q}\) seemed to have the most impact on insulin dynamics, we tried to see if there was a relationship between the two (figure 1).</p> |
<figure> | <figure> | ||
<img class='crispy center-block' src="https://static.igem.org/mediawiki/2017/9/9c/T--Sydney_Australia--insulin_dynamics.png" style="width:70%; margin:auto; padding:5px;"/> | <img class='crispy center-block' src="https://static.igem.org/mediawiki/2017/9/9c/T--Sydney_Australia--insulin_dynamics.png" style="width:70%; margin:auto; padding:5px;"/> | ||
Line 1,147: | Line 1,147: | ||
</figure> | </figure> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;">Now, since \(\color{#3e3f3f}{Q}\) in the model formed in [1] is the equilbrium constant of the reaction</p> |
<br> | <br> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> Insulin<sub>Hexamer </sub> \(\color{#3e3f3f}{\rightleftharpoons}\) Insulin<sub>Dimer</sub>, |
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;">it is related to the Gibbs free energy of the reaction by the expression \(\color{#3e3f3f}{\Delta G^{o}=-RT\ln{Q}}\), where \(\color{#3e3f3f}{R=8.314472 J K^{-1} mol^{-1}}\) is the gas constant and \(T\) is the temperature in kelvins.</p> |
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> Therefore if we know the Gibbs free energy of insulin hexamer formation, we can use this to find some qualitative information on the dynamics of insulin absorption using the model from [1], and thus estimate the time of peak insulin action and the duration of insulin action from thermodynamic information.</p> |
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> We then wanted to see if we could estimate the dynamics our insulin analogue (winsulin) would have without experimental work, which we did not have time for. This would require a computational estimate of the \(\Delta G\) of hexamer formation for our analogue. |
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> The <a href="https://www.ncbi.nlm.nih.gov/research/mutabind/index.fcgi/">Mutabind</a> tool [2] was used to model the effects of our variations to proinsulin's sequence on protein-protein interactions within the insulin hexamer. We used PDB file <a href="https://www.rcsb.org/pdb/explore.do?structureId=3AIY">3AIY</a> and inputted the sequence variants we had designed our winsulin with. Since the B chains are buried at the center of the insulin hexamer (see <a href="https://www.rcsb.org/pdb/explore.do?structureId=3AIY">3AIY</a>), we analysed the effects of the mutations in this chain on hexamer stability</p> |
<figure> | <figure> | ||
Line 1,163: | Line 1,163: | ||
</figure> | </figure> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> Mutabind predicted that all of our sequence changes would decrease the stability of the insulin hexamer for our analogue. Results are shown in table 3 </p> |
<table style="width:70%;margin:auto;"> | <table style="width:70%;margin:auto;"> | ||
Line 1,186: | Line 1,186: | ||
</table> | </table> | ||
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> Where \(\color{#3e3f3f}{\Delta\Delta G_{bind} (kcal mol^{-1})}\) is the predicted change in binding affinity induced by a mutation. A positive result corresponds to destabilising mutations, so the hexamer formation of winsulin will be less stable than that of human insulin.</p> |
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> This corresponds to a decrease in \(\color{#3e3f3f}{Q}\), meaning we predict our winsulin analogue will be relatively fast acting, compared to regular human insulin. Human insulin activity peaks in 2-4 hours and lasts for 6-8 hours [3] so this would make winsulin a rapid-acting analogue.</p> |
− | < | + | <h1 style="font-family:'Quicksand';font-size:16px;"> Although this is a crude estimate, it does give us some qualitative information on the action profile of our novel winsulin</p> |
<h1 style="font-family:'Quicksand';"> References </h1> | <h1 style="font-family:'Quicksand';"> References </h1> | ||
<ol> | <ol> | ||
− | <li>Tarin, C., Teufel, E., Pico, J., Bondia, J., Pfleiderer, H.J. (2005). Comprehensive pharmacokinetic model of insulin Glargine and other insulin formulations. IEEE Transactions on Biomedical Engineering, vol. 52, no. 12, pp. 1994-2005</li> | + | <li style="font-family:'Quicksand';">Tarin, C., Teufel, E., Pico, J., Bondia, J., Pfleiderer, H.J. (2005). Comprehensive pharmacokinetic model of insulin Glargine and other insulin formulations. IEEE Transactions on Biomedical Engineering, vol. 52, no. 12, pp. 1994-2005</li> |
− | <li>Li, M., Simonetti, F. L., Goncearenco, A., & Panchenko, A. R. (2016). MutaBind estimates and interprets the effects of sequence variants on protein–protein interactions. Nucleic Acids Research, 44(Web Server issue), W494–W501. http://doi.org/10.1093/nar/gkw374 </li> | + | <li style="font-family:'Quicksand';">Li, M., Simonetti, F. L., Goncearenco, A., & Panchenko, A. R. (2016). MutaBind estimates and interprets the effects of sequence variants on protein–protein interactions. Nucleic Acids Research, 44(Web Server issue), W494–W501. http://doi.org/10.1093/nar/gkw374 </li> |
− | <li>Diabetes Education Online. 2017. Types of Insulin. [ONLINE] Available at: https://dtc.ucsf.edu/types-of-diabetes/type2/treatment-of-type-2-diabetes/medications-and-therapies/type-2-insulin-rx/types-of-insulin/.</li> | + | <li style="font-family:'Quicksand';">Diabetes Education Online. 2017. Types of Insulin. [ONLINE] Available at: https://dtc.ucsf.edu/types-of-diabetes/type2/treatment-of-type-2-diabetes/medications-and-therapies/type-2-insulin-rx/types-of-insulin/.</li> |
</ol> | </ol> | ||
Revision as of 11:07, 26 October 2017
The aim of modelling in synthetic biology is to simulate the behaviour of your project to gain insight into how to best improve it. For our project, we saw three levels at which modelling could aid in the pursuit of its central aims.
As has been explored on our integrated human practices and applied design pages, the problem of insulin accessibility is complex and multi-faceted. As such, we decided it was not enough to consider our project as a problem whose solution could be found solely in a test tube. Distilled down, our project can be viewed as three sequential aims which we believe together can be used to address insulin accessibility.
Our modelling efforts were split into three branches, which reflected these major aspects of our project
Difficulties optimising production of recombinant are a key issue in the state of its accessibility.
in silico experiments to simulate how best to optimise expression led to theoretical insights which informed the direction of our efforts.
It is imperative to test the feasibility of our recombinant insulin as a therapeutic for diabetics.
Modelling the effects of changes to insulin’s biochemical makeup on its therapeutic effects supplement our wet-lab efforts to characterise our molecule
In addition, the project would be moot without a consideration of the insulin market as a whole.
Modelling helped us to gain insight into the global insulin market, which informed our approach towards entrepeneurship.
Whole Cell Model
The notation of Weisse et al has been adopted as follows:
Symbol | Meaning |
---|---|
\[ \color{#3e3f3f}{ s_i}\] | Internal Nutrient |
\[\color{#3e3f3f}{a}\] | energy, such as ATP |
\[\color{#3e3f3f}{r}\] | Ribosomes |
\[\color{#3e3f3f}{e_t}\] | A Transporter enzyme |
\[\color{#3e3f3f}{e_m}\] | A Metabolic Enzyme |
\[\color{#3e3f3f}{q}\] | House-keeping proteins |
\[\color{#3e3f3f}{m_x \textrm{ with } x\in\{r,t,m,q\}}\] | free mRNAs of the four species of proteins |
\[\color{#3e3f3f}{c_p\textrm{ with } x\in\{r,t,m,q\}}\] | ribosome-bound mRNA of the four species of proteins |
Symbol | Meaning |
---|---|
\[ \color{#3e3f3f}{ \upsilon_{imp}}\] | Rate of nutrient import |
\[\color{#3e3f3f}{\upsilon_{cat}}\] | Rate of nutrient metabolism |
\[\color{#3e3f3f}{\lambda}\] | Growth Rate |
\[\color{#3e3f3f}{n_x\textrm{ with } x\in\{r,t,m,q\}}\] | length of proteins of different species |
\[\color{#3e3f3f}{\upsilon_x \textrm{ with } x\in\{r,t,m,q\}}\] | Rate of translating protein species' |
\[\color{#3e3f3f}{k_b}\] | mRNA ribosome binding rate |
\[\color{#3e3f3f}{k_u}\] | mRNA ribosome unbinding rate |
\[\color{#3e3f3f}{\omega_x \textrm{ with } x\in\{r,t,m,q\}}\] | Transcription rates of the four species of proteins |
\[\color{#3e3f3f}{d_m}\] | mRNA degradation rate |
Symbol | Meaning |
---|---|
\[\color{#3e3f3f}{n_x\textrm{ with } x\in\{r,t,m,q\}}\] | length of proteins of different species |
\[\color{#3e3f3f}{d_m}\] | ribosome-bound mRNA of the four species of proteins |
Cytoplasmic Expression Model
Our first expression system was cytoplasmic expression in E. coli. We modelled the rate of change of five biochemical species in the cell:
Symbol | Meaning |
---|---|
\[\color{#3e3f3f}{m_p}\] | free mRNA of recombinant protein |
\[\color{#3e3f3f}{c_p}\] | ribosome-bound mRNA of recombinant protein |
\[\color{#3e3f3f}{p_u}\] | Unfolded recombinant protein |
\[\color{#3e3f3f}{p_f}\] | Folded recombinant protein |
\[\color{#3e3f3f}{p_a}\] | recombinant protein aggregated in inclusion bodies |
A diagram showing the species we modelled and notation used is shown below:
The following reactions were considered
Process | Reaction | Rate |
---|---|---|
Transcription | \[\color{#3e3f3f}{\varnothing\rightarrow m_p}\] | \[\color{#3e3f3f}{\omega_p(a)}\] |
Dilution and degradation of mRNA | \[\color{#3e3f3f}{m_p\rightarrow\varnothing}\] | \[\color{#3e3f3f}{\lambda+d_m}\] |
ribosome binding | \[\color{#3e3f3f}{r+m_p\rightleftharpoons c_p}\] | \[\color{#3e3f3f}{\textrm{forward: } k_b \textrm{, reverse: } k_u}\] |
Dilution of ribosome-bound protein | \[\color{#3e3f3f}{c_p\rightarrow\varnothing}\] | \[\color{#3e3f3f}{\lambda}\] |
Translation | \[\color{#3e3f3f}{n_pa+c_p\rightarrow m_p+p_u+r}\] | \[\color{#3e3f3f}{\upsilon_p(c_p,a)}\] |
Aggregation | \[\color{#3e3f3f}{p_u\rightarrow p_a}\] | \[\color{#3e3f3f}{k_a}\] |
Folding | \[\color{#3e3f3f}{p_u\rightarrow p_f}\] | \[\color{#3e3f3f}{k_f}\] |
Dilution and degradation of folded protein | \[\color{#3e3f3f}{p_f\rightarrow \varnothing} \] | \[\color{#3e3f3f}{\lambda+k_d}\] |
Here, \(\omega_p(a)\), the rate of transcription, is an energy dependent process.
We used the transcription rate form used in Weisse et al to denote the amount being transcribed (\(\omega_p(a)\)). That is,
\[\color{#3e3f3f}{\omega_p(a)=w_p \frac{a}{\theta_p+a} } \]
Where \(w_p\) is the maximal rate of transcription, dependent on the speed of transcriptional elongation, as well as the gene length, induction and copy number. \(a\) is the energy in the cell such as ATP (transcription is an energy dependent process), and
\(\theta_p\) is the transcriptional threshold of the recombinant protein.
In addition, we used the form in Weisse et al. for the translation rate term
\[\color{#3e3f3f}{\upsilon_p(c_p,a)=c_p \frac{\gamma(a)}{n_p} } \]
Where \(n_p\) is the length of recombinant protein, and \(\gamma(a)\) is an expression for the rate of transcriptional elongation:
\[\color{#3e3f3f}{\gamma(a)=\frac{\gamma_{max} a}{K_{\gamma} + a} } \]
Where \(\gamma_{max}\) is the maximal rate of translation, \(K_{\gamma}\) is the translational elongation threshold, and \(a\) is the energy in the cell.
For the model of inclusion body aggregation, we assumed first order deposition of monomers of unfolded protein, dependent on the concentration of unfolded protein. as in Hoffmann et al (2001).
Using the law of mass action kinetics we can derive a set of ordinary differential equations from these reactions.
Summary of Cytoplasmic Expression Model
\[ \color{#3e3f3f}{\frac{d}{dt}{m}_p=\omega_p(a)+\upsilon_p(c_p,a)+k_uc_p-(\lambda +d_m)m_p-k_brm_p} \]
\[ \color{#3e3f3f}{\frac{d}{dt}{c}_p=k_brm_p-\lambda c_p-k_uc_p-\upsilon_p(c_p,a)}\]
\[ \color{#3e3f3f}{\frac{d}{dt}{p}_u=\upsilon_p(c_p,a)-(k_f+k_a+\lambda)p_u}\]
\[ \color{#3e3f3f}{\frac{d}{dt}{p}_a=k_ap_u-\lambda p_a}\]
\[ \color{#3e3f3f}{\frac{d}{dt}{p}_f=k_fp_u-(k_d+\lambda) p_f}\]
Parametrising the model
The following table shows the parameters we needed to find for our model, and the values we used
Where \(w_p\) is the maximal rate of transcription, dependent on the speed of transcriptional elongation, as well as the gene length, induction and copy number. \(a\) is the energy in the cell such as ATP (transcription is an energy dependent process), and
\(\theta_p\) is the transcriptional threshold of the recombinant protein.
In addition, we used the form in Weisse et al. for the translation rate term
\[\color{#3e3f3f}{\upsilon_p(c_p,a)=c_p \frac{\gamma(a)}{n_p} } \]
Where \(n_p\) is the length of recombinant protein, and \(\gamma(a)\) is an expression for the rate of transcriptional elongation:
\[\color{#3e3f3f}{\gamma(a)=\frac{\gamma_{max} a}{K_{\gamma} + a} } \]
Where \(\gamma_{max}\) is the maximal rate of translation, \(K_{\gamma}\) is the translational elongation threshold, and \(a\) is the energy in the cell.
For the model of inclusion body aggregation, we assumed first order deposition of monomers of unfolded protein, dependent on the concentration of unfolded protein. as in Hoffmann et al (2001).
Using the law of mass action kinetics we can derive a set of ordinary differential equations from these reactions.
Summary of Cytoplasmic Expression Model
\[ \color{#3e3f3f}{\frac{d}{dt}{m}_p=\omega_p(a)+\upsilon_p(c_p,a)+k_uc_p-(\lambda +d_m)m_p-k_brm_p} \]
\[ \color{#3e3f3f}{\frac{d}{dt}{c}_p=k_brm_p-\lambda c_p-k_uc_p-\upsilon_p(c_p,a)}\]
\[ \color{#3e3f3f}{\frac{d}{dt}{p}_u=\upsilon_p(c_p,a)-(k_f+k_a+\lambda)p_u}\]
\[ \color{#3e3f3f}{\frac{d}{dt}{p}_a=k_ap_u-\lambda p_a}\]
\[ \color{#3e3f3f}{\frac{d}{dt}{p}_f=k_fp_u-(k_d+\lambda) p_f}\]
Parametrising the model
The following table shows the parameters we needed to find for our model, and the values we used
Where \(n_p\) is the length of recombinant protein, and \(\gamma(a)\) is an expression for the rate of transcriptional elongation:
\[\color{#3e3f3f}{\gamma(a)=\frac{\gamma_{max} a}{K_{\gamma} + a} } \]
Where \(\gamma_{max}\) is the maximal rate of translation, \(K_{\gamma}\) is the translational elongation threshold, and \(a\) is the energy in the cell.
For the model of inclusion body aggregation, we assumed first order deposition of monomers of unfolded protein, dependent on the concentration of unfolded protein. as in Hoffmann et al (2001).
Using the law of mass action kinetics we can derive a set of ordinary differential equations from these reactions.
Summary of Cytoplasmic Expression Model
\[ \color{#3e3f3f}{\frac{d}{dt}{m}_p=\omega_p(a)+\upsilon_p(c_p,a)+k_uc_p-(\lambda +d_m)m_p-k_brm_p} \]
\[ \color{#3e3f3f}{\frac{d}{dt}{c}_p=k_brm_p-\lambda c_p-k_uc_p-\upsilon_p(c_p,a)}\]
\[ \color{#3e3f3f}{\frac{d}{dt}{p}_u=\upsilon_p(c_p,a)-(k_f+k_a+\lambda)p_u}\]
\[ \color{#3e3f3f}{\frac{d}{dt}{p}_a=k_ap_u-\lambda p_a}\]
\[ \color{#3e3f3f}{\frac{d}{dt}{p}_f=k_fp_u-(k_d+\lambda) p_f}\]
Parametrising the model
The following table shows the parameters we needed to find for our model, and the values we used
For the model of inclusion body aggregation, we assumed first order deposition of monomers of unfolded protein, dependent on the concentration of unfolded protein. as in Hoffmann et al (2001).
Using the law of mass action kinetics we can derive a set of ordinary differential equations from these reactions.
Summary of Cytoplasmic Expression Model
\[ \color{#3e3f3f}{\frac{d}{dt}{m}_p=\omega_p(a)+\upsilon_p(c_p,a)+k_uc_p-(\lambda +d_m)m_p-k_brm_p} \]
\[ \color{#3e3f3f}{\frac{d}{dt}{c}_p=k_brm_p-\lambda c_p-k_uc_p-\upsilon_p(c_p,a)}\]
\[ \color{#3e3f3f}{\frac{d}{dt}{p}_u=\upsilon_p(c_p,a)-(k_f+k_a+\lambda)p_u}\]
\[ \color{#3e3f3f}{\frac{d}{dt}{p}_a=k_ap_u-\lambda p_a}\]
\[ \color{#3e3f3f}{\frac{d}{dt}{p}_f=k_fp_u-(k_d+\lambda) p_f}\]
Parametrising the model
The following table shows the parameters we needed to find for our model, and the values we used
Summary of Cytoplasmic Expression Model
\[ \color{#3e3f3f}{\frac{d}{dt}{m}_p=\omega_p(a)+\upsilon_p(c_p,a)+k_uc_p-(\lambda +d_m)m_p-k_brm_p} \] \[ \color{#3e3f3f}{\frac{d}{dt}{c}_p=k_brm_p-\lambda c_p-k_uc_p-\upsilon_p(c_p,a)}\] \[ \color{#3e3f3f}{\frac{d}{dt}{p}_u=\upsilon_p(c_p,a)-(k_f+k_a+\lambda)p_u}\] \[ \color{#3e3f3f}{\frac{d}{dt}{p}_a=k_ap_u-\lambda p_a}\] \[ \color{#3e3f3f}{\frac{d}{dt}{p}_f=k_fp_u-(k_d+\lambda) p_f}\]Parametrising the model
The following table shows the parameters we needed to find for our model, and the values we used
Symbol | Meaning | Default value | Units | Source |
---|---|---|---|---|
\[\color{#3e3f3f}{w_p}\] | Maximal rate of transcription | <10^3 | mRNAs/min | Proportional to induction level. Varied around realistic values as recommended by[1] |
\[\color{#3e3f3f}{\theta_p}\] | transcriptional threshold of the recombinant protein | 4.38 | [molecs/cell] | [1] |
\[\color{#3e3f3f}{n_p}\] | Length of recombinant protein | 312/255 | [aa/molecs] | Length of cytoplasmic proinsulin/winsulin gblock *link to design/parts page?* |
\[\color{#3e3f3f}{\gamma_{max}}\] | Maximal rate of translation | 1260 | [aa/ min molecs] | [1] |
\[\color{#3e3f3f}{K_{\gamma}}\] | Translational elongation threshold | 7 | [molecs/ cell] | [1] |
\[\color{#3e3f3f}{k_u}\] | Rate of unbinding of mRNA and ribosomes | 1 | [/min] | [1] |
\[\color{#3e3f3f}{k_b}\] | Rate of binding of mRNA and ribosomes | 1 | [cell/ min molecs] | [1] |
\[\color{#3e3f3f}{d_m}\] | degradation rate of mRNA | 0.1 | [/min] | [1] |
\[\color{#3e3f3f}{k_f}\] | Rate of protein folding | 0.14 | [/min] | adapted to fit units from [2] |
\[\color{#3e3f3f}{k_a}\] | Rate of protein aggregation | 0.21 | [/min] | adapted to fit units from [2] |
\[\color{#3e3f3f}{k_d}\] | Rate of protein degradation | 0 | * |
Periplasmic Expression Model
Next we looked at periplasmic expression of our recombinant protein in E. coli. We modelled the rate of change of 6 species
Symbol | Meaning |
---|---|
\[\color{#3e3f3f}{m_p}\] | free mRNA of recombinant protein |
\[\color{#3e3f3f}{c_p}\] | ribosome-bound mRNA of recombinant protein |
\[\color{#3e3f3f}{p_c}\] | Unfolded recombinant protein in the cytoplasm |
\[\color{#3e3f3f}{p_t}\] | Unfolded recombinant protein bound to transporter |
\[\color{#3e3f3f}{p_u}\] | Unfolded recombinant protein in the periplasm |
\[\color{#3e3f3f}{p_f}\] | Folded recombinant protein in the periplasm |
A diagram showing the species we modelled and notation used is shown below:
The following reactions were considered
Process | Reaction | Rate |
---|---|---|
Transcription | \[\color{#3e3f3f}{\varnothing\rightarrow m_p}\] | \[\color{#3e3f3f}{\omega_p(a)}\] |
Dilution and degradation of mRNA | \[\color{#3e3f3f}{m_p\rightarrow\varnothing}\] | \[\color{#3e3f3f}{\lambda+d_m}\] |
ribosome binding | \[\color{#3e3f3f}{r+m_p\rightleftharpoons c_p}\] | \[\color{#3e3f3f}{\textrm{forward: } k_b \textrm{, reverse: } k_u}\] |
Dilution of ribosome-bound protein | \[\color{#3e3f3f}{c_p\rightarrow\varnothing}\] | \[\color{#3e3f3f}{\lambda}\] |
Translation | \[\color{#3e3f3f}{n_pa+c_p\rightarrow m_p+p_u+r}\] | \[\color{#3e3f3f}{\upsilon_p(c_p,a)}\] |
Translocator binding | \[\color{#3e3f3f}{p_c+t\rightarrow p_t}\] | \[\color{#3e3f3f}{k_bt}\] |
Translocation | \[\color{#3e3f3f}{p_t\rightarrow p_u}\] where \(t\) refers to the amount of translocons | \[\color{#3e3f3f}{\tau(p_t,a)}\] |
Folding | \[\color{#3e3f3f}{p_u\rightarrow p_f}\] | \[\color{#3e3f3f}{k_f}\] |
Dilution and degradation of folded protein | \[\color{#3e3f3f}{p_f\rightarrow \varnothing} \] | \[\color{#3e3f3f}{\lambda+k_d}\] |
here, \(\omega(a)\) and \(\upsilon(c_p,a)\) are as in the cytoplasmic reactions.
The amount being transported is found with the term \(\tau_p(p_t,a)\).
Protein translocation to the periplasm occurs via an ATP-dependent motor protein, secA [4].
Post-translational translocation uses ATP as a stepwisesource of energy to drive the protein through the membrane.
It follows mechanism illustrated in Figure 1 [4].
Following the logic used to derive the translation rate in [1], we derive the net rate of translocating a protein \(p\) by defining \(K_p:=\frac{k_1k_2}{k_{-1}+k_2}\). This leads to
\[\color{#3e3f3f}{\tau_p(p_t,a)=p_t\Big(\frac{n_p}{50}\Big(\frac{1}{K_pa}+\frac{1}{k_2} \Big)+\frac{1}{k_t}\Big)^{-1}}\]
If we assume the final termination step is fast, so \(\frac{1}{k_t}<< \frac{n_p}{50}\Big(\frac{1}{K_pa}+\frac{1}{k_2} \Big) \), this is approximately equal to
\[\color{#3e3f3f}{\tau_p(p_t,a)\approx 50p_t \frac{\epsilon(a)}{n_p}\qquad \epsilon(a):=\frac{\epsilon_{max}a}{K_{\epsilon}+a} }\]
Where \(\epsilon_{max}\) is the maximal translocation rate, \(K_{\epsilon}\) is the threshold, and \(n_p\) is the length of the protein in amino acids
Parametrising Translocation
To find the parameters for translocation (\(\epsilon(a)\)) and (\(K_{\epsilon}\)), we used kinetic parameters determined in [5]. They measured translocation of a 346aa protein proOmpA and found the apparent Km of SecA was 50nM, and the threshold was 2.7 proOmpa/site/min. This converts to 2.7 \(\cdot\) 346 proOmpA/site/min aa/proOmpa \(\rightarrow\) 934.2 aa/molec/min
Using the law of mass action kinetics we can derive a set of ordinary differential equations from these reactions.
Summary of Periplasmic Expression Model
\[ \color{#3e3f3f}{\frac{d}{dt}{m}_p=\omega_p(a)+\upsilon_p(c_p,a)+k_uc_p-(\lambda +d_m)m_p-k_brm_p} \]
\[ \color{#3e3f3f}{\frac{d}{dt}p=k_brm_p-\lambda c_p-k_uc_p-\upsilon_p(c_p,a)}\]
\[ \color{#3e3f3f}{\frac{d}{dt}c=\upsilon_p(c_p,a)-(k_{bt}t+\lambda)p_c}\]
\[\color{#3e3f3f}{\frac{d}{dt}t=k_{bt}tp_c-\tau_p(p_t,a)-\lambda p_t}\]
\[\color{#3e3f3f}{\frac{d}{dt}u=\tau_p(p_t,a)-(k_f+\lambda) p_u}\]
\[\color{#3e3f3f}{\frac{d}{dt}f=k_fp_u-(k_d+\lambda)p_f}\]
Parametrising the model
The following table shows the parameters we needed to find for our model, and the values we used
If we assume the final termination step is fast, so \(\frac{1}{k_t}<< \frac{n_p}{50}\Big(\frac{1}{K_pa}+\frac{1}{k_2} \Big) \), this is approximately equal to
\[\color{#3e3f3f}{\tau_p(p_t,a)\approx 50p_t \frac{\epsilon(a)}{n_p}\qquad \epsilon(a):=\frac{\epsilon_{max}a}{K_{\epsilon}+a} }\]
Where \(\epsilon_{max}\) is the maximal translocation rate, \(K_{\epsilon}\) is the threshold, and \(n_p\) is the length of the protein in amino acids
Parametrising Translocation
To find the parameters for translocation (\(\epsilon(a)\)) and (\(K_{\epsilon}\)), we used kinetic parameters determined in [5]. They measured translocation of a 346aa protein proOmpA and found the apparent Km of SecA was 50nM, and the threshold was 2.7 proOmpa/site/min. This converts to 2.7 \(\cdot\) 346 proOmpA/site/min aa/proOmpa \(\rightarrow\) 934.2 aa/molec/min
Using the law of mass action kinetics we can derive a set of ordinary differential equations from these reactions.
Summary of Periplasmic Expression Model
\[ \color{#3e3f3f}{\frac{d}{dt}{m}_p=\omega_p(a)+\upsilon_p(c_p,a)+k_uc_p-(\lambda +d_m)m_p-k_brm_p} \]
\[ \color{#3e3f3f}{\frac{d}{dt}p=k_brm_p-\lambda c_p-k_uc_p-\upsilon_p(c_p,a)}\]
\[ \color{#3e3f3f}{\frac{d}{dt}c=\upsilon_p(c_p,a)-(k_{bt}t+\lambda)p_c}\]
\[\color{#3e3f3f}{\frac{d}{dt}t=k_{bt}tp_c-\tau_p(p_t,a)-\lambda p_t}\]
\[\color{#3e3f3f}{\frac{d}{dt}u=\tau_p(p_t,a)-(k_f+\lambda) p_u}\]
\[\color{#3e3f3f}{\frac{d}{dt}f=k_fp_u-(k_d+\lambda)p_f}\]
Parametrising the model
The following table shows the parameters we needed to find for our model, and the values we used
Parametrising Translocation
To find the parameters for translocation (\(\epsilon(a)\)) and (\(K_{\epsilon}\)), we used kinetic parameters determined in [5]. They measured translocation of a 346aa protein proOmpA and found the apparent Km of SecA was 50nM, and the threshold was 2.7 proOmpa/site/min. This converts to 2.7 \(\cdot\) 346 proOmpA/site/min aa/proOmpa \(\rightarrow\) 934.2 aa/molec/min
Using the law of mass action kinetics we can derive a set of ordinary differential equations from these reactions.
Summary of Periplasmic Expression Model
\[ \color{#3e3f3f}{\frac{d}{dt}{m}_p=\omega_p(a)+\upsilon_p(c_p,a)+k_uc_p-(\lambda +d_m)m_p-k_brm_p} \]
\[ \color{#3e3f3f}{\frac{d}{dt}p=k_brm_p-\lambda c_p-k_uc_p-\upsilon_p(c_p,a)}\]
\[ \color{#3e3f3f}{\frac{d}{dt}c=\upsilon_p(c_p,a)-(k_{bt}t+\lambda)p_c}\]
\[\color{#3e3f3f}{\frac{d}{dt}t=k_{bt}tp_c-\tau_p(p_t,a)-\lambda p_t}\]
\[\color{#3e3f3f}{\frac{d}{dt}u=\tau_p(p_t,a)-(k_f+\lambda) p_u}\]
\[\color{#3e3f3f}{\frac{d}{dt}f=k_fp_u-(k_d+\lambda)p_f}\]
Parametrising the model
The following table shows the parameters we needed to find for our model, and the values we used
Summary of Periplasmic Expression Model
\[ \color{#3e3f3f}{\frac{d}{dt}{m}_p=\omega_p(a)+\upsilon_p(c_p,a)+k_uc_p-(\lambda +d_m)m_p-k_brm_p} \] \[ \color{#3e3f3f}{\frac{d}{dt}p=k_brm_p-\lambda c_p-k_uc_p-\upsilon_p(c_p,a)}\] \[ \color{#3e3f3f}{\frac{d}{dt}c=\upsilon_p(c_p,a)-(k_{bt}t+\lambda)p_c}\] \[\color{#3e3f3f}{\frac{d}{dt}t=k_{bt}tp_c-\tau_p(p_t,a)-\lambda p_t}\] \[\color{#3e3f3f}{\frac{d}{dt}u=\tau_p(p_t,a)-(k_f+\lambda) p_u}\] \[\color{#3e3f3f}{\frac{d}{dt}f=k_fp_u-(k_d+\lambda)p_f}\]Parametrising the model
The following table shows the parameters we needed to find for our model, and the values we used
Symbol | Meaning | Default value | Units | Source |
---|---|---|---|---|
\[\color{#3e3f3f}{w_p}\] | Maximal rate of transcription | <10^3 | mRNAs/min | Proportional to induction level. Varied around realistic values as recommended by[1] |
\[\color{#3e3f3f}{\theta_p}\] | transcriptional threshold of the recombinant protein | 4.38 | [molecs/cell] | [1] |
\[\color{#3e3f3f}{n_p}\] | Length of recombinant protein | 312/255 | [aa/molecs] | Length of cytoplasmic proinsulin/winsulin gblock *link to design/parts page?* |
\[\color{#3e3f3f}{\gamma_{max}}\] | Maximal rate of translation | 1260 | [aa/ min molecs] | [1] |
\[\color{#3e3f3f}{K_{\gamma}}\] | Translational elongation threshold | 7 | [molecs/ cell] | [1] |
\[\color{#3e3f3f}{k_u}\] | Rate of unbinding of mRNA and ribosomes | 1 | [/min] | [1] |
\[\color{#3e3f3f}{k_b}\] | Rate of binding of mRNA and ribosomes | 1 | [cell/ min molecs] | [1] |
\[\color{#3e3f3f}{d_m}\] | degradation rate of mRNA | 0.1 | [/min] | [1] |
\[\color{#3e3f3f}{t}\] | Number of translocons in a cell | 500 | [/cell] | [5] |
\[\color{#3e3f3f}{k_{bt}}\] | Rate of protein binding to translocon | 1 | [cell /min molecs] | [1] |
\[\color{#3e3f3f}{\epsilon_{max}}\] | Maximal translocation rate | 934.2 | [aa /min molecs] | [5] |
\[\color{#3e3f3f}{K_{\epsilon}}\] | Translocational threshold | [molecs/ cell] | [5] | |
\[\color{#3e3f3f}{k_f}\] | Rate of protein folding | 0.28 | [/min] | † |
\[\color{#3e3f3f}{k_d}\] | Rate of protein degradation | 0 | * |
Symbol | Meaning | Default value | Units | Source |
---|---|---|---|---|
\[\color{#3e3f3f}{w_p}\] | Maximal rate of transcription | <10^3 | mRNAs/min | Proportional to induction level. Varied around realistic values as recommended by[1] |
\[\color{#3e3f3f}{\theta_p}\] | transcriptional threshold of the recombinant protein | 4.38 | [molecs/cell] | [1] |
\[\color{#3e3f3f}{n_p}\] | Length of recombinant protein | 312/255 | [aa/molecs] | Length of cytoplasmic proinsulin/winsulin gblock *link to design/parts page?* |
\[\color{#3e3f3f}{\gamma_{max}}\] | Maximal rate of translation | 1260 | [aa/ min molecs] | [1] |
\[\color{#3e3f3f}{K_{\gamma}}\] | Translational elongation threshold | 7 | [molecs/ cell] | [1] |
\[\color{#3e3f3f}{k_u}\] | Rate of unbinding of mRNA and ribosomes | 1 | [/min] | [1] |
\[\color{#3e3f3f}{k_b}\] | Rate of binding of mRNA and ribosomes | 1 | [cell/ min molecs] | [1] |
\[\color{#3e3f3f}{d_m}\] | degradation rate of mRNA | 0.1 | [/min] | [1] |
\[\color{#3e3f3f}{k_{bt}}\] | Rate of protein binding to translocon | [cell /min molecs] | ||
\[\color{#3e3f3f}{\epsilon_{max}}\] | Maximal translocation rate | [aa /min molecs] | ||
\[\color{#3e3f3f}{K_{\epsilon}}\] | Translocational threshold | [molecs/ cell] | ||
\[\color{#3e3f3f}{k_f}\] | Rate of protein folding | 0.14 | [/min] | adapted to fit units from [2] |
\[\color{#3e3f3f}{k_d}\] | Rate of protein degradation | 0 | * |
Bacillus Secretory Expression Model
We also developed a model of our secretory protein expression system in bacillus subtilis. The model included 6 species
Symbol | Meaning |
---|---|
\[\color{#3e3f3f}{m_p}\] | free mRNA of recombinant protein |
\[\color{#3e3f3f}{c_p}\] | ribosome-bound mRNA of recombinant protein |
\[\color{#3e3f3f}{p_c}\] | Unfolded recombinant protein in the cytoplasm |
\[\color{#3e3f3f}{p_t}\] | Unfolded recombinant protein bound to transporter |
\[\color{#3e3f3f}{p_u}\] | Unfolded recombinant protein in the medium |
\[\color{#3e3f3f}{p_f}\] | Folded recombinant protein in the medium |
A diagram showing the species we modelled and notation used is shown below:
Structurally, this is the same process as the periplasmic expression system, so the equations' structure is the same. However the parameters are different, reflecting the different environment of bacillus and medium and its effect on expression of recombinant protein
Summary of Bacillus Secretory Expression Model
\[ \color{#3e3f3f}{\frac{d}{dt}{m}_p=\omega_p(a)+\upsilon_p(c_p,a)+k_uc_p-(\lambda +d_m)m_p-k_brm_p} \]
\[ \color{#3e3f3f}{\frac{d}{dt}{c}_p=k_brm_p-\lambda c_p-k_uc_p-\upsilon_p(c_p,a)}\]
\[ \color{#3e3f3f}{\frac{d}{dt}{p}_c=\upsilon_p(c_p,a)-(k_{bt}t+\lambda)p_c}\]
\[\color{#3e3f3f}{\frac{d}{dt}{p}_t=k_{bt}tp_c-\tau_p(p_t,a)-\lambda p_t}\]
\[\color{#3e3f3f}{\frac{d}{dt}{p}_u=\tau_p(p_t,a)-(k_f+\lambda) p_u}\]
\[\color{#3e3f3f}{\frac{d}{dt}{p}_f=k_fp_u-(k_d+\lambda)p_f}\]
Summary of Bacillus Secretory Expression Model
\[ \color{#3e3f3f}{\frac{d}{dt}{m}_p=\omega_p(a)+\upsilon_p(c_p,a)+k_uc_p-(\lambda +d_m)m_p-k_brm_p} \] \[ \color{#3e3f3f}{\frac{d}{dt}{c}_p=k_brm_p-\lambda c_p-k_uc_p-\upsilon_p(c_p,a)}\] \[ \color{#3e3f3f}{\frac{d}{dt}{p}_c=\upsilon_p(c_p,a)-(k_{bt}t+\lambda)p_c}\] \[\color{#3e3f3f}{\frac{d}{dt}{p}_t=k_{bt}tp_c-\tau_p(p_t,a)-\lambda p_t}\] \[\color{#3e3f3f}{\frac{d}{dt}{p}_u=\tau_p(p_t,a)-(k_f+\lambda) p_u}\] \[\color{#3e3f3f}{\frac{d}{dt}{p}_f=k_fp_u-(k_d+\lambda)p_f}\]We were unfortunately unable to parametrise the bacillus model, so for our in silico experiments we focused on comparing cytoplasmic and periplasmic E. coli expression.
In Silico Experiments
Once we had modelled our different expression systems for recombinant insulin, we integrated them into the whole cell model developed in [1]. The methodology behind integrating models of our expression system into this model was to more accurately reflect reality. Recombinant protein expression occurs within a complex cellular environment with finite resources. A model which ignores the actiities of the host cells would ignore important host-circuit interactions. Ignoring the finite resources of the cell may skew our prediction of the yield of our expression systems.
We then interrogated these models for insights into how to optimise the expression of insulin
Comparing Cytoplasmic and Periplasmic Expression
Comparing Cytoplasmic and Periplasmic Expression
References
For our physiological modelling, we used a model of subcutaneous insulin absorption developed in [1] and used it to develop a model relating the free energy of insulin hexamer formation and insulin dynamics (the time of peak of action, and the duration of action).
The authors of [1] developed a system of partial differential equations to describe the insulin infusion process. They modelled the change in three species:
Symbol | Meaning |
---|---|
\[\color{#3e3f3f}{c_d}\] | Insulin in dimeric form |
\[\color{#3e3f3f}{c_h}\] | Insulin in hexamer form |
\[\color{#3e3f3f}{c_b}\] | Insulin in bound form |
They modelled the conversion between hexameric and dimeric insulin as follows
InsulinHexamer \(\color{#3e3f3f}{\rightleftharpoons}\) InsulinDimer
Where the forward rate was called \(\color{#3e3f3f}{P}\) and the reverse rate was \(\color{#3e3f3f}{PQ}\) where we can interpret \(\color{#3e3f3f}{P}\) as the production rate and \(\color{#3e3f3f}{Q}\) as the equilibrium constant.
The final model was as follows
\[\color{#3e3f3f}{\eqalignno{{\partial c_{d}(t,r)\over\partial t}=&\,P\left(c_{h}(t,r)-Qc_{d}(t,r)^{3}\right)-B_{d}c_{d}(t,r)\cr&+D\nabla^{2}c_{d}(t,r),\cr{\partial c_{h}(t,r)\over\partial t}=&\,-P\left(c_{h}(t,r)-Qc_{d}(t,r)^{3}\right)\cr&+D\nabla^{2}c_{h}(t,r)}}\]
Where \(\color{#3e3f3f}{P, Q, B_d, D}\) are parameters, and exogenous insulin flow is obtained by integrating the expression denoting the amount of insulin dimer entering the bloodstream:
\[\color{#3e3f3f}{I_{ex}(t)=B_{d}\int\limits_{V_{sc}}c_{d}(t,r)dV.}\]
The parameters found for the different insulin analogues and their impact on insulin dynamics are
Where the forward rate was called \(\color{#3e3f3f}{P}\) and the reverse rate was \(\color{#3e3f3f}{PQ}\) where we can interpret \(\color{#3e3f3f}{P}\) as the production rate and \(\color{#3e3f3f}{Q}\) as the equilibrium constant.
The final model was as follows
\[\color{#3e3f3f}{\eqalignno{{\partial c_{d}(t,r)\over\partial t}=&\,P\left(c_{h}(t,r)-Qc_{d}(t,r)^{3}\right)-B_{d}c_{d}(t,r)\cr&+D\nabla^{2}c_{d}(t,r),\cr{\partial c_{h}(t,r)\over\partial t}=&\,-P\left(c_{h}(t,r)-Qc_{d}(t,r)^{3}\right)\cr&+D\nabla^{2}c_{h}(t,r)}}\]
Where \(\color{#3e3f3f}{P, Q, B_d, D}\) are parameters, and exogenous insulin flow is obtained by integrating the expression denoting the amount of insulin dimer entering the bloodstream:
\[\color{#3e3f3f}{I_{ex}(t)=B_{d}\int\limits_{V_{sc}}c_{d}(t,r)dV.}\]
The parameters found for the different insulin analogues and their impact on insulin dynamics are
Where \(\color{#3e3f3f}{P, Q, B_d, D}\) are parameters, and exogenous insulin flow is obtained by integrating the expression denoting the amount of insulin dimer entering the bloodstream:
\[\color{#3e3f3f}{I_{ex}(t)=B_{d}\int\limits_{V_{sc}}c_{d}(t,r)dV.}\]
The parameters found for the different insulin analogues and their impact on insulin dynamics are
Insulin Analogue | \(\color{#3e3f3f}{Q}\) | \(\color{#3e3f3f}{D}\) | \(\color{#3e3f3f}{B_d}\) | Time of peak Insulin action (hours) | Duration of Insulin action (hours) |
---|---|---|---|---|---|
Lispro, Humalog, NovoRapid | \[\color{#3e3f3f}{4.75\cdot 10^{-4}}\] | \[\color{#3e3f3f}{3.36\cdot 10^{-4}}\] | \[\color{#3e3f3f}{2.36\cdot 10^{-2}}\] | \[\color{#3e3f3f}{0.25}\] | \[\color{#3e3f3f}{4}\] |
Actrapid | \[\color{#3e3f3f}{1.9\cdot 10^{-3}}\] | \[\color{#3e3f3f}{8.4\cdot 10^{-5}}\] | \[\color{#3e3f3f}{1.18\cdot 10^{-2}}\] | \[\color{#3e3f3f}{0.75}\] | \[\color{#3e3f3f}{8}\] |
Semilente | \[\color{#3e3f3f}{7.6\cdot 10^{-2}}\] | \[\color{#3e3f3f}{8.4\cdot 10^{-5}}\] | \[\color{#3e3f3f}{1.18\cdot 10^{-2}}\] | \[\color{#3e3f3f}{1.3}\] | \[\color{#3e3f3f}{11}\] |
NPH | \[\color{#3e3f3f}{3.04}\] | \[\color{#3e3f3f}{8.4\cdot 10^{-5}}\] | \[\color{#3e3f3f}{1.18\cdot 10^{-2}}\] | \[\color{#3e3f3f}{4.5}\] | \[\color{#3e3f3f}{16}\] |
Since the parameter \(\color{#3e3f3f}{Q}\) seemed to have the most impact on insulin dynamics, we tried to see if there was a relationship between the two (figure 1).
Now, since \(\color{#3e3f3f}{Q}\) in the model formed in [1] is the equilbrium constant of the reaction
InsulinHexamer \(\color{#3e3f3f}{\rightleftharpoons}\) InsulinDimer,
it is related to the Gibbs free energy of the reaction by the expression \(\color{#3e3f3f}{\Delta G^{o}=-RT\ln{Q}}\), where \(\color{#3e3f3f}{R=8.314472 J K^{-1} mol^{-1}}\) is the gas constant and \(T\) is the temperature in kelvins.
Therefore if we know the Gibbs free energy of insulin hexamer formation, we can use this to find some qualitative information on the dynamics of insulin absorption using the model from [1], and thus estimate the time of peak insulin action and the duration of insulin action from thermodynamic information.
We then wanted to see if we could estimate the dynamics our insulin analogue (winsulin) would have without experimental work, which we did not have time for. This would require a computational estimate of the \(\Delta G\) of hexamer formation for our analogue.
The Mutabind tool [2] was used to model the effects of our variations to proinsulin's sequence on protein-protein interactions within the insulin hexamer. We used PDB file 3AIY and inputted the sequence variants we had designed our winsulin with. Since the B chains are buried at the center of the insulin hexamer (see 3AIY), we analysed the effects of the mutations in this chain on hexamer stability
Mutabind predicted that all of our sequence changes would decrease the stability of the insulin hexamer for our analogue. Results are shown in table 3
InsulinHexamer \(\color{#3e3f3f}{\rightleftharpoons}\) InsulinDimer,
it is related to the Gibbs free energy of the reaction by the expression \(\color{#3e3f3f}{\Delta G^{o}=-RT\ln{Q}}\), where \(\color{#3e3f3f}{R=8.314472 J K^{-1} mol^{-1}}\) is the gas constant and \(T\) is the temperature in kelvins.
Therefore if we know the Gibbs free energy of insulin hexamer formation, we can use this to find some qualitative information on the dynamics of insulin absorption using the model from [1], and thus estimate the time of peak insulin action and the duration of insulin action from thermodynamic information.
We then wanted to see if we could estimate the dynamics our insulin analogue (winsulin) would have without experimental work, which we did not have time for. This would require a computational estimate of the \(\Delta G\) of hexamer formation for our analogue.
The Mutabind tool [2] was used to model the effects of our variations to proinsulin's sequence on protein-protein interactions within the insulin hexamer. We used PDB file 3AIY and inputted the sequence variants we had designed our winsulin with. Since the B chains are buried at the center of the insulin hexamer (see 3AIY), we analysed the effects of the mutations in this chain on hexamer stability
Mutabind predicted that all of our sequence changes would decrease the stability of the insulin hexamer for our analogue. Results are shown in table 3
Therefore if we know the Gibbs free energy of insulin hexamer formation, we can use this to find some qualitative information on the dynamics of insulin absorption using the model from [1], and thus estimate the time of peak insulin action and the duration of insulin action from thermodynamic information.
We then wanted to see if we could estimate the dynamics our insulin analogue (winsulin) would have without experimental work, which we did not have time for. This would require a computational estimate of the \(\Delta G\) of hexamer formation for our analogue.
The Mutabind tool [2] was used to model the effects of our variations to proinsulin's sequence on protein-protein interactions within the insulin hexamer. We used PDB file 3AIY and inputted the sequence variants we had designed our winsulin with. Since the B chains are buried at the center of the insulin hexamer (see 3AIY), we analysed the effects of the mutations in this chain on hexamer stability
Mutabind predicted that all of our sequence changes would decrease the stability of the insulin hexamer for our analogue. Results are shown in table 3
The Mutabind tool [2] was used to model the effects of our variations to proinsulin's sequence on protein-protein interactions within the insulin hexamer. We used PDB file 3AIY and inputted the sequence variants we had designed our winsulin with. Since the B chains are buried at the center of the insulin hexamer (see 3AIY), we analysed the effects of the mutations in this chain on hexamer stability
Mutabind predicted that all of our sequence changes would decrease the stability of the insulin hexamer for our analogue. Results are shown in table 3
Mutation | \(\color{#3e3f3f}{\Delta\Delta G_{bind} (kcal mol^{-1})}\) |
---|---|
\[\color{#3e3f3f}{H10D}\] | \[\color{#3e3f3f}{0.51}\] |
\[\color{#3e3f3f}{T27S}\] | \[\color{#3e3f3f}{0.57}\] |
\[\color{#3e3f3f}{K295}\] | \[\color{#3e3f3f}{0.62}\] |