{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Simple transcription-translation system modeled by a deterministic ODE system"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## Biochemical reactions\n",
"\n",
"In a simple transcription-translation system the RNA of a certain gene is prodcued with a fixed rate $\\beta$ and decays with a rate proportional to the curren amount of RNA, a so called first order decay. The protein that our gene is coding for is produced with a rate that is proportional to the current amount of RNA and follows a first order decay as well. A reaction scheme of this small system looks like this:\n",
"\n",
"$$\\xrightarrow{k_{1}} RNA \\xrightarrow{d_{1}} o$$ $$RNA \\xrightarrow{k_{2}} protein$$ $$protein \\xrightarrow{d_{2}} o$$\n",
"\n",
"## Deterministic ODE system\n",
"\n",
"The above reaction scheme can be tranlated into two [differential equations](https://en.wikipedia.org/wiki/Differential_equation). For the sake of simplicity we have assumed [mass action kinetics](https://en.wikipedia.org/wiki/Law_of_mass_action) without any saturation effects.\n",
"\n",
"$$ \\frac{\\mathrm{d}RNA}{\\mathrm{d}t} = k_{1} - d_{1}\\cdot RNA$$ \n",
"\n",
"$$ \\frac{\\mathrm{d}protein}{\\mathrm{d}t} = k_{2}\\cdot RNA - d_{2} \\cdot protein$$\n",
"\n",
"This system of equations describes the rule how the systems state given by the concentration of RNA and protein molecules changes given the current state. The left hand side of the equations are the rate of change and the right hand side the rule how the change depends on the current state. Note that the concentration of RNA appears in both equations, thus the equations are coupled as the amount of proteins depends on the amount of RNA.\n",
"\n",
"To explicitly solve ODEs it is necessary to know the state of the system at one instant of time. Often refered to as the intial conditions and associated with time point $t = 0$. \n",
"\n",
"We assumed the transcription rate $\\beta$ to be constant over time. You can imagine this situation as if at time point zero a transcription factor binds the gene's promoter and activates transcription. For now we will investigate the temporal evolution of the RNA and protein assuming constant transcription, later we will consider the case when transcription is only transiently activated.\n",
"\n",
"### Steady state\n",
"\n",
"The transcription-translation reaches a steady state, i.e. RNA and protein concentrations that do not change over time because the respective production and decay processes balance each other. In steady state the left hand side of the differential equations equals zero and we are left with two algebraic equations that we can solve to calculate the steady state:\n",
"\n",
"**Exercise 1**\n",
"\n",
"Solve the algebraic equations for the steady state: $RNA_{ss}$ and $P_{ss}$.\n",
"\n",
"$$RNA_{ss} = \\frac{k_{1}}{d_{1}}$$\n",
"\n",
"$$P_{ss} = \\frac{k_{2}\\cdot RNA_{ss}}{d_{2}} = \\frac{k_{1}k_{2}}{d_{1}d_{2}}$$\n",
"\n",
"The steady state is fully determined by the biochemical reaction rates of our system and here it is a stable state. This means if we move the system out of steady state it will return.\n",
"\n",
"## Analytic solution\n",
"\n",
"The solution of the ODE system with the intital condition $RNA_{0} = protein_{0} = 0$ is given by:\n",
"\n",
"$$ \\text{RNA}(t)= RNA_{ss}\\big( 1- e^{-d_{1} t}\\big)$$\n",
"\n",
"$$\\text{protein}(t)= P_{ss}\\bigg(1 - \\frac{d_{1}e^{-d_{2}t} - d_{2}e^{-d_{2}t}}{d_{1} - d_{2}}\\bigg)$$\n",
"\n",
"For the simple ODE system given here it is possible to obtain a analytical solution, in more complex situations this usually is not possible and it is necessary to numerically integrate differential equations. Although we could work here with the analytical solutions we want to go the numerical way. A first impression on how an algorithm could approximate a solution of an ODE gives [Euler method](https://en.wikipedia.org/wiki/Euler_method)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# imports\n",
"import pylab as pl\n",
"import numpy as np\n",
"import scipy as sp\n",
"import matplotlib\n",
"from matplotlib import gridspec\n",
"from scipy.optimize import fsolve\n",
"from scipy.integrate import odeint\n",
"# inline displaying\n",
"%matplotlib inline\n",
"# settings for plots\n",
"matplotlib.rcParams['axes.labelsize'] = 16\n",
"matplotlib.rcParams['xtick.labelsize'] = 16\n",
"matplotlib.rcParams['ytick.labelsize'] = 16\n",
"matplotlib.rcParams['legend.fontsize'] = 14\n",
"matplotlib.rcParams['font.family'] = ['sans-serif']\n",
"line_width = 2\n",
"color_rna = 'dodgerblue'\n",
"color_protein = 'darkorange'\n",
"color_phase_space = 'dimgrey'\n",
"color_initial_state = 'crimson'\n",
"color_steady_state = 'saddlebrown'"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# functions\n",
"\n",
"# numerical solution of the ODE system\n",
"def solve_ODE_transcription_translation(x,t,param):\n",
" \n",
" RNA = x[0]\n",
" protein = x[1]\n",
" \n",
" # define the two coupled ODEs\n",
" dRNA = param[0] - param[1]*RNA\n",
" dprotein = param[2]*RNA - param[3]*protein \n",
" \n",
" # stack the results together into one array with the first row being the RNA and the second the protein\n",
" return np.hstack((dRNA,dprotein))\n",
"\n",
"# analytical solution based on the equations above\n",
"def analytic_solution(tt,param):\n",
" \n",
" # steady state\n",
" RNA_ss = param[0]/param[1]\n",
" p_ss = param[0]/param[1]*param[2]/param[3]\n",
" \n",
" # analytic ODE solution\n",
" RNA = RNA_ss*(1 - sp.exp(-param[1]*tt))\n",
" protein = p_ss*(1 - (param[1]*sp.exp(-param[3]*tt) - param[3]*sp.exp(-param[1]*tt))/(param[1] - param[3]))\n",
" \n",
" return np.vstack((RNA,protein))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# define parameters to solve the ODEs\n",
"k_1 = 1\n",
"d_1 = 0.15\n",
"k_2 = 0.15\n",
"d_2 = 0.1\n",
"param = [k_1,d_1,k_2,d_2]\n",
"\n",
"# calculate steady state\n",
"steady_state_RNA = k_1/d_1\n",
"steady_state_protein = k_2*k_1/d_1/d_2\n",
"\n",
"# define intial conditions\n",
"RNA_0 = 0\n",
"protein_0 = 0\n",
"y_0 = [RNA_0,protein_0]\n",
"\n",
"# set the time vector\n",
"tt = sp.linspace(0,100,200)\n",
"\n",
"# solve the ODE system with the defined parameters and initial condition\n",
"sol_deterministic = odeint(solve_ODE_transcription_translation,y_0,tt,args = (param,))\n",
"\n",
"# calculate the analytic solution\n",
"analytic = analytic_solution(tt,param)\n",
"\n",
"# plot results\n",
"fig,ax = pl.subplots(1,1, figsize = (10,6))\n",
"ax.plot(tt,sol_deterministic[:,0], label = 'RNA, numeric', lw = line_width,color = color_rna)\n",
"ax.plot(tt,analytic[0],'--',lw = 2*line_width, color = color_rna, label = 'RNA, analytic')\n",
"ax.plot(tt,sol_deterministic[:,1], label = 'Protein, numeric', lw = line_width, color = color_protein)\n",
"ax.plot(tt,analytic[1],'--',lw = 2*line_width, color = color_protein, label = 'Protein, analytic')\n",
"ax.plot(tt,np.zeros(len(tt)) + steady_state_RNA,'--',color = color_rna,label = 'steady state RNA')\n",
"ax.plot(tt,np.zeros(len(tt)) + steady_state_protein,'--',color = color_protein, label = 'steady state protein')\n",
"ax.legend(loc = 'best')\n",
"ax.set_xlabel('Time')\n",
"ax.set_ylabel('Concentration')\n",
"ax.set_title('Time course')\n",
"ax.set_ylim(0,steady_state_protein+steady_state_protein*0.1)"
]
},
{
"cell_type": "markdown",
"metadata": {
"code_folding": []
},
"source": [
"## Influence of the RNA and protein decay rates on the time to reach the steady state\n",
"\n",
"**Exercise 2**\n",
"\n",
"Investigate an effect of different RNA and protein rates on kinetics of the given transcription-translation system\n",
"\n",
"**Instructions**\n",
"- Solve ODE systems with the defined parameters and initial conditions\n",
"- Investigate the influence of different parameters on the system's behaviour by perturbing ```k_1, d_1, d_2``` (in II-IV conditions) \n",
"\n",
"**Questions**\n",
"- Which rate(s) influence the steady state?\n",
"- Which rate(s) influence the kinetics, i.e. how fast the steady state is reached?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# define intial conditions \n",
"RNA_0 = 0\n",
"protein_0 = 0\n",
"y_0 = [RNA_0,protein_0]\n",
"\n",
"# set the time vector\n",
"tt = sp.linspace(0,100,200)\n",
"\n",
"# define parameter vector to solve the ODEs\n",
"\n",
"# I) Reference parameters\n",
"k_1 = 1\n",
"d_1 = 0.15\n",
"k_2 = 0.15\n",
"d_2 = 0.1\n",
"param = [k_1,d_1,k_2,d_2]\n",
"# ODE system\n",
"sol_deterministic = odeint(solve_ODE_transcription_translation,y_0,tt,args = (param,))\n",
"#---------------------------------\n",
"\n",
"# II) Perturbed parameters (RNA decay)\n",
"k_1 = \n",
"d_1 = \n",
"k_2 =\n",
"d_2 = \n",
"param1 = \n",
"sol_deterministic1 = \n",
"#---------------------------------\n",
"\n",
"# III) Perturbed parameters (protein decay)\n",
"k_1 = \n",
"d_1 = \n",
"k_2 =\n",
"d_2 = \n",
"param2 = \n",
"sol_deterministic2 = \n",
"#---------------------------------\n",
"\n",
"# IV) Perturbed parameters (RNA synthesis)\n",
"k_1 = \n",
"d_1 = \n",
"k_2 = \n",
"d_2 = \n",
"param3 = \n",
"sol_deterministic3 = \n",
"#---------------------------------\n",
"\n",
"# plot results\n",
"fig,ax = pl.subplots(4,2, figsize = (17,17))\n",
"ax[0,0].plot(tt,sol_deterministic[:,0], label = 'RNA, condition 1', lw = line_width,color = color_rna)\n",
"\n",
"ax[0,0].legend(loc = 'best')\n",
"ax[0,0].set_xlabel('Time')\n",
"ax[0,0].set_ylabel('Concentration')\n",
"\n",
"ax[0,1].plot(tt,sol_deterministic[:,1], label = 'Protein, condition 1', lw = line_width, color = color_protein)\n",
"ax[0,1].legend(loc = 'best')\n",
"ax[0,1].set_xlabel('Time')\n",
"ax[0,1].set_ylabel('Concentration')\n",
"\n",
"ax[1,0].plot(tt,sol_deterministic1[:,0], label = 'RNA, condition 2', lw = line_width,color = color_rna)\n",
"ax[1,0].legend(loc = 'best')\n",
"ax[1,0].set_xlabel('Time')\n",
"ax[1,0].set_ylabel('Concentration')\n",
"\n",
"ax[1,1].plot(tt,sol_deterministic1[:,1], label = 'Protein, condition 2', lw = line_width, color = color_protein)\n",
"ax[1,1].legend(loc = 'best')\n",
"ax[1,1].set_xlabel('Time')\n",
"ax[1,1].set_ylabel('Concentration')\n",
"\n",
"ax[2,0].plot(tt,sol_deterministic2[:,0], label = 'RNA, condition 3', lw = line_width,color = color_rna)\n",
"ax[2,0].legend(loc = 'best')\n",
"ax[2,0].set_xlabel('Time')\n",
"ax[2,0].set_ylabel('Concentration')\n",
"\n",
"ax[2,1].plot(tt,sol_deterministic2[:,1], label = 'Protein, condition 3', lw = line_width, color = color_protein)\n",
"ax[2,1].legend(loc = 'best')\n",
"ax[2,1].set_xlabel('Time')\n",
"ax[2,1].set_ylabel('Concentration')\n",
"\n",
"ax[3,0].plot(tt,sol_deterministic3[:,0], label = 'RNA, condition 4', lw = line_width,color = color_rna)\n",
"ax[3,0].legend(loc = 'best')\n",
"ax[3,0].set_xlabel('Time')\n",
"ax[3,0].set_ylabel('Concentration')\n",
"\n",
"ax[3,1].plot(tt,sol_deterministic3[:,1], label = 'Protein, condition 4', lw = line_width, color = color_protein)\n",
"ax[3,1].legend(loc = 'best')\n",
"ax[3,1].set_xlabel('Time')\n",
"ax[3,1].set_ylabel('Concentration')\n",
"\n",
"\n",
"pl.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's test the effect of the protein decay rates on the systems kinetcs. This time we define a set of decay rates, keep the steady state fixed by adjusting the potein production rate $k$ and then iterate over the defined rates in a for loop."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# range of protein decay rates\n",
"D_2 = np.array([0.02,0.05,0.1,0.2,0.3,0.4,0.5])\n",
"K_2 = D_2 * steady_state_protein / steady_state_RNA\n",
"\n",
"k_1 = 1\n",
"d_1 = 0.15\n",
"\n",
"tt = sp.linspace(0,100,1000)\n",
"sol = []\n",
"for ii,d_2 in enumerate(D_2):\n",
" param = [k_1,d_1,K_2[ii],d_2]\n",
" s = odeint(solve_ODE_transcription_translation,y_0,tt,args = (param,))\n",
" sol.append(s)\n",
" \n",
"# plot results\n",
"fig,ax = pl.subplots(1,2,figsize = (17,6))\n",
"alpha = sp.linspace(0.2,1,len(sol))\n",
"for ii,ss in enumerate(sol):\n",
" ax[0].plot(tt,ss[:,0],lw = line_width, color = color_rna, alpha = alpha[ii], \n",
" label = '$d_{2} = $'+str(D_2[ii]))\n",
" \n",
" ax[1].plot(tt,ss[:,1],lw = line_width, color = color_protein, alpha = alpha[ii], \n",
" label = '$d_{2} = $'+str(D_2[ii]))\n",
"\n",
"ax[0].legend(loc = 'best')\n",
"ax[0].set_xlabel('Time')\n",
"ax[0].set_ylabel('RNA')\n",
"\n",
"ax[1].legend(loc = 'best')\n",
"ax[1].set_xlabel('Time')\n",
"ax[1].set_ylabel('Protein')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The time to reach half steady state provides a measure for the response time of the system. The response time of the RNA depends only on RNA decay and can be calculated as:\n",
"\n",
"$$\\tau_{1/2, RNA} = \\frac{\\ln 2}{d_{1}}$$\n",
"\n",
"As the plots above indicate the response time of the protein depends on both decay rates $\\gamma_{m}$ and $\\gamma_{p}$. The time to reach steady state is determined by\n",
"\n",
"$$ \\frac{d_{1}e^{-d_{2}t} - d_{2}e^{-d_{1}t}}{d_{1} - d_{2}} - \\frac{1}{2} = 0$$\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Transient activation of gene transcription\n",
"\n",
"Now we will investigate the case when gene transcription is only transiently activated by a transcription factor which in turn might me activated by a external stimulus. The biological scenario is that the transcription factor remains at the promoter only for a certain amount of time. Mathematically this means that the transcription rate is a function of time $k_1 = k_1(t)$. We assume the functional relation for $k_1$ to be:\n",
"\n",
"$$k_1(t) = \\beta\\cdot e^{-r_1 t}\\cdot (1 - e^{-r_{2}t})$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# functions\n",
"def TF_binding(param,tt):\n",
" r_1,r_2 = param[0],param[1]\n",
" return sp.exp(-k_1*tt)*(1-sp.exp(-k_2*tt))\n",
"\n",
"# ODE system with time varying transcription rate\n",
"def solve_ODE_time_dependent_transcription(x,t,param,param_binding):\n",
" RNA = x[0]\n",
" dRNA = param[0]*TF_binding(param_binding,t) - param[1]*RNA\n",
" return dRNA"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercise 3**\n",
"\n",
"Tune the RNA decay rate to create different behaviour of the RNA time courses. Try to resemble the following expression characteristics:\n",
"- fast and transient (cluster I)\n",
"- slow and sustained (cluster III)\n",
"- intermediate behaviour (cluster II)\n",
"as in the example from the paper [PMID:19198593]. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# define the binding pattern\n",
"r_1 = 1.5\n",
"r_2 = 0.1\n",
"param_binding = [k_1,k_2]\n",
"tt = sp.linspace(0,5,1000)\n",
"binding = TF_binding(param_binding,tt)\n",
"binding = binding/binding.max()\n",
"y_0 = 0\n",
"\n",
"# define parameters \n",
"k_1 = 1\n",
"d_1 = 10\n",
"param = [k_1,d_1]\n",
"\n",
"# solve the ODE system\n",
"sol_fast = odeint(solve_ODE_time_dependent_transcription,y_0,tt,args = (param,param_binding))\n",
"sol_fast = sol_fast/sol_fast[:,0].max()\n",
"\n",
"# plot the result\n",
"fig,ax = pl.subplots(1,1, figsize = (17,6))\n",
"\n",
"ax.plot(tt,sol_fast[:,0], label = 'RNA, decay rate', lw = line_width, color = color_initial_state)\n",
"ax.plot(tt,binding,'--',lw = line_width,color = color_phase_space,label = 'Binding pattern')\n",
"ax.legend(loc = 'best')\n",
"ax.set_xlabel('Time (h)')\n",
"ax.set_ylabel('Concentration')\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Correlation between RNA and protein abundance depends on the RNA decay rate\n",
"\n",
"**Exercise 4** Investigate the correlation between RNA and protein levels. \n",
"- Plot RNA and protein levels agains each other (use simulations you have created earlier (Exercise 2))\n",
"\n",
"**Questions**\n",
"- Under what conditions RNA and protein are strongly correlated?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"fig,ax = pl.subplots(1,1, figsize = (8,8))\n",
"\n",
"# enter the data to plot here\n",
"ax.plot()\n",
"ax.plot()\n",
"ax.set_ylabel()\n",
"ax.set_xlabel()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Circadian oscillations as an example of negative feedback\n",
"\n",
"Circadian oscillations are our daily rhythms: activity, metabolism, sleep, etc. They are controlled by multiple genes that influence each other. The simple transcription-translation system that we were investigating so far is not able to create oscillations. For sustained oscillation we need negative feedback. This means that the protein acts back on its own transcription as a supressor. A second important ingredient to create oscillations is a time delay, thus, the feedback of the protein (here _Per1_) acts not instantaneous but with some delay. The protein _Per1_ exists in two forms here indicated as $protein$ and $protein_p$ which supresses its own transcription but first it has to be created which generates the delay. \n",
"\n",
"Here are the differential equations describing the negative feedback loop of _Per1_ on its own transcription. The feedback on the RNA prodiction is mediated by a special input function, a so called Hill-Function. A sigmoidal function, i.e. the effect of the feedback saturates and does not become arbitrarily large. The Hill-Coefficient $n$ describes the steepness of the Hill-Function. We will see during the investigation of the system that $n$ plays a cruicial role and that it has to pass a minimal threshold to obtain undampened oscillations.\n",
"\n",
"$$\\frac{dRNA}{dt} = k_{1}\\frac{K_{i}^{n}}{K_{i}^{n} + protein_{p}^{n}} - k_{2} RNA $$\n",
"\n",
"$$\\frac{dprotein}{dt} = k_{3} RNA - k_4 protein$$\n",
"\n",
"$$\\frac{dprotein_{p}}{dt} = k_{5}protein - k_{6}protein_{p}$$\n",
"\n",
"**Exercise 5**\n",
"\n",
"Implement negative feedback model and plot Z(t) for the following parameters: \n",
"$$k_{1} = k_{3} = k_{5} = 1$$\n",
"$$k_{2} = k_{4} = k_{6} = 0.1$$\n",
"$$K_{i} = 1,$$\n",
"$$n = 10$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def circadian_ODEs(x,t,param):\n",
" k_1,k_2,k_3,k_4,k_5,k_6,K_i,n = param[0],param[1],param[2],param[3],param[4],param[5],param[6],param[7]\n",
" RNA = x[0]\n",
" protein = x[1]\n",
" protein_p = x[2]\n",
" \n",
" # write here the right hand side of the ODE\n",
" dRNA = \n",
" dprotein = \n",
" dprotein_p = \n",
" \n",
" return np.hstack((dRNA,dprotein,dprotein_p))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# define parameters\n",
"k_1 = \n",
"k_2 = \n",
"k_3 = \n",
"k_4 = \n",
"k_5 =\n",
"k_6 = \n",
"K_i = \n",
"n = \n",
"param = [k_1,k_2,k_3,k_4,k_5,k_6,K_i,n]\n",
"\n",
"# define intital states\n",
"RNA_0 = 0\n",
"protein_0 = 1\n",
"protein_p_0 = 4\n",
"y_0 = [RNA_0,protein_0,protein_p_0]\n",
"\n",
"tt = sp.linspace(0,500,1000)\n",
"\n",
"# generate simulations \n",
"sol_circ = odeint(circadian_ODEs,y_0,tt,args = (param,))\n",
"\n",
"# plot results \n",
"fig,ax = pl.subplots(1,1, figsize = (10,6))\n",
"ax.plot(tt,sol_circ[:,0], lw = line_width, color = color_rna, label = 'RNA')\n",
"ax.plot(tt,sol_circ[:,1], lw = line_width, color = color_protein, label = 'Protein')\n",
"ax.plot(tt,sol_circ[:,2], lw = line_width, color = color_initial_state, label = '$Protein_{p}$')\n",
"ax.legend(loc = 'best')\n",
"ax.set_ylabel('Concentration')\n",
"ax.set_xlabel('Time')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Effect of the protein decay rates and the cooperativity on the oscillation frequency\n",
"\n",
"**Exercise 6**\n",
"Change the protein decay rates to ```k_4 = k_6 = 0.2```\n",
"\n",
"**Exercise 7** \n",
"Change cooperativity factor to ```n = 3```\n",
"\n",
"**Questions**\n",
"- What is the effect on the oscillation period?\n",
"- What is the effect on the dampening of the oscillations?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Transcription-translation system modeled as stochastic system: single cell behaviour\n",
"\n",
"Biochemical interactions in single cells can be random. For instance, various transcription factors are present in cells only in hundreds to thousands of copies, which is not much considering the size of the genome and that each gene is usually present in two alleles only. Such low copy numbers can create strong stochastic effects on the single cell level. In this part of the lecture we want to briefly dig into its implications. We use a so called Monte Carlo algorithm to simulate stochastic realistations of our transcription-translation system. We will test the effects of different parameter values on temporal fluctuations."
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## The stochastic simulation algorithm\n",
"\n",
"Developed mainly by Daniel Gillespie. The algorithm simulates the stochastic biochemical interactions of the reacting species, here RNA and protein, based on the current amount molecules of each species and the rates of the possible reactions. The algorithm essentially is a loop and in each iteration the next reaction that happens and the time point of the reaction are simulated based on random number generation. A deeper explanation of the algorithm can be found on [wikipedia](https://en.wikipedia.org/wiki/Gillespie_algorithm) or in the original publications of [D. Gillespie](http://pubs.acs.org/doi/abs/10.1021/j100540a008). [Here](https://www.ncbi.nlm.nih.gov/pubmed/17037977) you can find a comprehensible review of the algorithm, it's foundations and possible approximate implementations.\n",
"\n",
"In contrast to the deterministic ODE system above the amount of RNA and protein now is not measured in concentrations but in the number of molecules of each species. The systems state is defined by the set of molecule counts. During each iteration of the algorithm the systems state is update according to the sampled reaction.\n",
"\n",
"The system consists of four different reactions:\n",
"\n",
"1. **RNA production:** $ RNA \\xrightarrow{\\beta} RNA + 1$\n",
"2. **RNA decay:** $RNA \\xrightarrow{\\gamma_{m}\\cdot RNA} RNA - 1$\n",
"3. **Protein production:** $protein \\xrightarrow{k\\cdot RNA} protein + 1$\n",
"4. **Protein decay:** $protein \\xrightarrow{\\gamma_{p}\\cdot protein} protein - 1$\n",
"\n",
"The likelihood of each reaction depends of the current state of the system, i.e. the number of present RNA and protein molecules.\n",
"\n",
"The **steady state** of the stochastic system is the same as for the deterministic ODE system. However, now we will observe fluctuations around the steady state. The characteristics of this fluctuations are controlled by the reaction rates. Time course measurements of fluctuations around the steady state allow to estimate the kinetic reaction rates which would not be for the ODE system."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# implementation of the Gillespie algorithm to simulate the stochastic transcription-tranlsation system\n",
"def stoch_sim_transcription_translation(param,initial_state,tf):\n",
" beta,gamma_m,k,gamma_p = param[0],param[1],param[2],param[3]\n",
" \n",
" # reaction matrix, the systems state includes the time points of the reactions in the first column\n",
" reactions = np.array([[0,1,0],\n",
" [0,-1,0],\n",
" [0,0,1],\n",
" [0,0,-1]])\n",
" \n",
" # initialise the systems state\n",
" state = np.zeros(3)\n",
" state[1] = initial_state[0]\n",
" state[2] = initial_state[1]\n",
" STATE = state\n",
" \n",
" tt = 0\n",
" while tt <= tf:\n",
" # sample two random numbers uniformly between 0 and 1\n",
" rr = sp.random.uniform(0,1,2)\n",
" \n",
" a_0 = beta + gamma_m*state[1] + k*state[1] + gamma_p*state[2] \n",
" a_s = np.array([beta,gamma_m*state[1],k*state[1],gamma_p*state[2]],dtype = float)\n",
" \n",
" # time step \n",
" tt = tt + 1. / a_0 * sp.log(1. / rr[0])\n",
" state[0] = tt\n",
" \n",
" # find the next reaction\n",
" prop = rr[1] * a_0\n",
" cum_a_s = np.cumsum(a_s)\n",
" \n",
" ind = np.where(prop <= cum_a_s)[0][0]\n",
" \n",
" # update the systems state\n",
" state = state+reactions[ind]\n",
" \n",
" STATE = np.vstack((STATE,state))\n",
" \n",
" return STATE"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# define parameters \n",
"k_1 = 2\n",
"d_1 = 0.15\n",
"k_2 = 0.2\n",
"d_2 = 0.1\n",
"param = [k_1,d_1,k_2,d_2]\n",
"\n",
"# calculate steady states\n",
"steady_state_RNA = k_1/d_1\n",
"steady_state_protein = k_2*k_1/d_1/d_2\n",
"\n",
"# define intial conditions\n",
"RNA_0 = 0\n",
"protein_0 = 0\n",
"initial_state = [RNA_0,protein_0]\n",
"\n",
"# simulation time\n",
"tf = 100\n",
"\n",
"sim = stoch_sim_transcription_translation(param,initial_state,tf)\n",
"\n",
"# determinisitic simulation\n",
"tt = sp.linspace(0,100,200)\n",
"y_0 = [RNA_0,protein_0]\n",
"sol_deterministic = odeint(solve_ODE_transcription_translation,y_0,tt,args = (param,))\n",
"\n",
"fig,ax = pl.subplots(1,2, figsize = (15,7))\n",
"ax[0].plot(sim[:,0],sim[:,1], label = 'RNA', color = color_rna, lw = line_width,drawstyle = 'steps')\n",
"ax[0].plot(tt,sol_deterministic[:,0], color = color_rna, lw = 2*line_width, label = 'RNA, ODE solution')\n",
"ax[0].plot(sim[:,0],sim[:,2], label = 'Protein', color = color_protein,lw = line_width,drawstyle = 'steps')\n",
"ax[0].plot(tt,sol_deterministic[:,1], color = color_protein, lw = 2*line_width, label = 'Protein, ODE solution')\n",
"ax[0].plot(sim[:,0],np.zeros(len(sim[:,0])) + steady_state_RNA,'--', label = 'steady state RNA',\n",
" color = color_rna)\n",
"ax[0].plot(sim[:,0],np.zeros(len(sim[:,0])) + steady_state_protein,'--', label = 'steady state protein',\n",
" color = color_protein)\n",
"ax[0].legend(loc = 'best')\n",
"ax[0].set_xlabel('Time')\n",
"ax[0].set_ylabel('Molecule number')\n",
"ax[0].set_ylim(0,2*steady_state_protein)\n",
"ax[0].set_xlim(0,130)\n",
"\n",
"ax[1].plot(sim[:,1],sim[:,2], label = 'state trajectory', lw = line_width,color = color_phase_space)\n",
"ax[1].plot(steady_state_RNA,steady_state_protein,'o',ms = 15,label = 'steady state', color = color_steady_state)\n",
"ax[1].plot(RNA_0,protein_0,'o',ms = 15,label = 'Initial state', color = color_initial_state)\n",
"ax[1].legend(loc = 'best')\n",
"ax[1].set_ylabel('Protein')\n",
"ax[1].set_xlabel('RNA')\n",
"ax[1].set_title('Phase space')"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"# Gene expression noise depends on transcription and translation rates\n",
"\n",
"Due to the stochastic nature of biochemical reactions noise in the expression level of proteins is inherent. The noise level can be tuned by a cell by controlling the transcription and translation rates of a gene. On our simple we want to exemplify the effects on two examples:\n",
"\n",
"1. High transcription rate and low tranlation rate\n",
"2. Low transcription rate and High translation rate\n",
"\n",
"The steady state protein level for both cases is the same, but the fluctuations around it are very different. \n",
"\n",
"**Exercise 8**\n",
"\n",
"Below you find a simulation of our system for the first case with high transcription rate and low translation rate. \n",
"- Create a second realisation with low transcription and high translation rates.\n",
"- Estimate the expression noise by means of the Coefficient of variation $CV= \\frac{STD}{Mean}$.\n",
"\n",
"**Questions**\n",
"- How should you change the rates to keep the protein steady state fixed?\n",
"- Which effect do the inverse changes of the transcription and translation rates have on temporal fluctuations? What did you expect on your previous experience?\n",
"- What are mechanisms for a cell to either promote or suppress expression noise?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# transcription rate\n",
"k_1 = 1\n",
"\n",
"# translation rate\n",
"k_2 = 0.1\n",
"\n",
"# decay rates\n",
"d_1 = 0.15\n",
"d_2 = 0.1\n",
"\n",
"param_k_1_high = [k_1,d_1,k_2,d_2]\n",
"\n",
"# simulation time\n",
"tf = 5000\n",
"\n",
"# set the initial state to the steady state\n",
"initial_state_k_1_high = [0,0]\n",
"\n",
"# run stochastic simulation\n",
"sim_k_1_high = stoch_sim_transcription_translation(param_k_1_high,initial_state_k_1_high,tf)\n",
"\n",
"# calculate coefficient of variation\n",
"CV_k_1_high = (sim_k_1_high[:,2].std()/sim_k_1_high[:,2].mean())\n",
"\n",
"fig,ax = pl.subplots(1,2, figsize = (17,5))\n",
"\n",
"ax[0].plot(sim_k_1_high[:,0],sim_k_1_high[:,2],lw = line_width,color = color_protein, \n",
" label = 'Transcription rate, $CV$ = '+str(np.round(CV_k_1_high,2))+', $k_1$ = '+str(k_1)+', $k_2$ = '+str(k_2))\n",
"ax[0].set_xlim(0,500)\n",
"ax[0].legend(loc = 'best')\n",
"ax[0].set_ylabel('Protein')\n",
"ax[0].set_xlabel('Time')\n",
"\n",
"ax[1].hist(sim_k_1_high[:,2], color = color_protein);\n",
"ax[1].set_ylabel('Frequency')\n",
"ax[1].set_xlabel('Protein')\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## Influence of the RNA decay rate on temporal fluctuations\n",
"\n",
"We will run stochastic simulations of the fluctuations around steady state and test the influence of the decay rates of RNA and protein while keeping the steady state constant. First lets start with the influence of the decay rate of the RNA. We define a range of RNA decay rates $\\gamma_{m}$ and keep the RNA steady state fixed by adjusting the RNA production rate."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# define kinetic parameters\n",
"steady_state_RNA = 10\n",
"steady_state_protein = 20\n",
"\n",
"# range of decay rates of the RNA\n",
"D_1 = np.array([0.01,0.05,0.1,0.5])\n",
"k_1 = steady_state_RNA * D_1\n",
"\n",
"d_2 = k_1 * steady_state_RNA / steady_state_protein\n",
"\n",
"k_2 = 0.2\n",
"\n",
"# set the inital conditions to the steady state\n",
"initial_state = np.round([steady_state_RNA,steady_state_protein])\n",
"\n",
"# set the simulation time longer than before for better statistics\n",
"tf = 500\n",
"\n",
"# iteration over the defined d_1\n",
"res = []\n",
"for ii,d_1 in enumerate(D_1):\n",
" param = [k_1[ii],d_1,k_2,d_2[ii]]\n",
" sim = stoch_sim_transcription_translation(param,initial_state,tf)\n",
" res.append(sim)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# visualize the results\n",
"fig,ax = pl.subplots(len(res),1, figsize = (9,8))\n",
"\n",
"color = 'dimgrey'\n",
"color_ss = 'darkorange'\n",
"xmin,xmax = 0,tf\n",
"ymin_m,ymax_m = 0,30\n",
"ymin_p,ymax_p = 0,40\n",
"\n",
"for ii,rr in enumerate(res):\n",
" ax[ii].plot(rr[:,0],rr[:,1],color = color, lw = 2,\n",
" label = '$\\gamma_{m}$ = '+str(D_1[ii]))\n",
" ax[ii].plot([0,tf],[steady_state_RNA,steady_state_RNA],'--',lw = 2, color = color_ss)\n",
" ax[ii].set_ylim(ymin_m,ymax_m)\n",
" ax[ii].set_xlim(xmin,xmax)\n",
" ax[ii].set_ylabel('RNA')\n",
" ax[ii].set_title('$\\gamma_{m}$ = '+str(D_1[ii]))\n",
" \n",
"ax[ii].set_xlabel('Time')\n",
" \n",
"pl.tight_layout()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"hide_input": false,
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.13"
}
},
"nbformat": 4,
"nbformat_minor": 2
}