{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "# imports\n", "import numpy as np\n", "import scipy as sp\n", "from scipy.integrate import odeint\n", "from scipy.optimize import fsolve\n", "import matplotlib\n", "from matplotlib import gridspec\n", "%matplotlib inline\n", "import pylab as pl\n", "\n", "# define some 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", "\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": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "# Introduction\n", "\n", "In this lecture you will learn some aspects of dynamic modeling applied to biological systems. We will investigate the properties of a simple transcription-translation system and its implications on real biological questions. The language of dynamical models are ordinary differential equations (ODEs). You will learn a bit what they are and how to solve them numerically using python libraries. We will make use of the concept of python functions to solve the same ODEs with differen parameter sets. A function usually looks like this:\n", "\n", "```python\n", "def function_name(param_1,param_2):\n", " some python code...\n", " to calculate result\n", " return result\n", "```\n", "\n", "In addition to plain python we rely on three libraries for scientific computing and visualisation: NumPy, SciPy and Matplotlib. Probably across this course you came across them already. " ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "# Simple transcription-translation system modeled by a deterministic ODE system" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true, "deletable": true, "editable": 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{\\beta} RNA \\xrightarrow{\\gamma_{m}} o$$ $$RNA \\xrightarrow{k} protein$$ $$protein \\xrightarrow{\\gamma_{p}} 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} = \\beta - \\gamma_{m}\\cdot RNA$$ \n", "\n", "$$ \\frac{\\mathrm{d}protein}{\\mathrm{d}t} = k\\cdot RNA - \\gamma_{p} \\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", "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^{-\\gamma_{m} t}\\big)$$\n", "\n", "$$\\text{protein}(t)= P_{ss}\\bigg(1 - \\frac{\\gamma_{m}e^{-\\gamma_{p}t} - \\gamma_{m}e^{-\\gamma_{p}t}}{\\gamma_{m} - \\gamma_{p}}\\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, "deletable": true, "editable": true }, "outputs": [], "source": [ "# functions\n", "\n", "# numerical solution of the ODE system\n", "def solve_ODE_transcription_translation(x,t,param):\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", " # 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": { "code_folding": [], "collapsed": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "# define parameters and time vector to solve the ODEs\n", "beta = 1\n", "gamma_m = 0.15\n", "k = 0.15\n", "gamma_p = 0.1\n", "param = [beta,gamma_m,k,gamma_p]\n", "\n", "# calculate steady state\n", "steady_state_RNA = beta/gamma_m\n", "steady_state_protein = k*beta/gamma_m/gamma_p\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 the result\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": [], "deletable": true, "editable": true }, "source": [ "## Influence of the RNA and protein decay rates on the time to reach the steady state\n", "\n", "**Exercise 2**\n", "\n", "Please try the effects of different combinations of RNA and protein rates on the kinetics of our small transcription ans translation system by changing the parameters ```gamma_m``` and ```gamma_p``` in the code example below. \n", "\n", "Reference parameters:\n", "\n", "- $\\beta = 1$\n", "- $\\gamma_{m} = 0.15$\n", "- $k = 0.15$\n", "- $\\gamma_{p} = 0.1$\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?\n", "\n", "**Instructions**\n", "- Below you see a plot of the same simulations as above, use them as references. Create a new simulation were you change the parameters ```beta, gamma_m``` and ```gamma_p```. Leave the translation rate ```k``` fixed.\n", "- Add your results into the plot below the reference simulation and investigate the influence of the different parameters\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "deletable": true, "editable": 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", "#--------------------------------\n", "\n", "# I) define parameters and time vector to solve the ODEs (Reference parameters)\n", "beta = 1\n", "gamma_m = 0.15\n", "k = 0.15\n", "gamma_p = 0.1\n", "param = [beta,gamma_m,k,gamma_p]\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", "#---------------------------------\n", "\n", "# II) define parameters and time vector to solve the ODEs (Perturbed parameters)\n", "# beta = 1\n", "# gamma_m = 0.15\n", "# k = 0.15\n", "# gamma_p = 0.1\n", "# param1 = [beta,gamma_m,k,gamma_p]\n", "\n", "# # solve the ODE system with the defined parameters and initial condition\n", "# sol_deterministic1 = odeint(solve_ODE_transcription_translation,y_0,tt,args = (param1,))\n", "\n", "#---------------------------------\n", "\n", "# plot the result\n", "fig,ax = pl.subplots(2,2, figsize = (17,12))\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')" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "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": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "# range of protein decay rates\n", "GAMMA_p = np.array([0.02,0.05,0.1,0.2,0.3,0.4,0.5])\n", "k = GAMMA_p * steady_state_protein / steady_state_RNA\n", "\n", "beta = 1\n", "gamma_m = 0.15\n", "\n", "tt = sp.linspace(0,100,1000)\n", "sol = []\n", "for ii,gamma_p in enumerate(GAMMA_p):\n", " param = [beta,gamma_m,k[ii],gamma_p]\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 = '$\\gamma_{p} = $'+str(GAMMA_p[ii]))\n", " \n", " ax[1].plot(tt,ss[:,1],lw = line_width, color = color_protein, alpha = alpha[ii], \n", " label = '$\\gamma_{p} = $'+str(GAMMA_p[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": { "deletable": true, "editable": true }, "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}{\\gamma_{m}}$$\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{\\gamma_{m}e^{-\\gamma_{p}t} - \\gamma_{p}e^{-\\gamma_{m}t}}{\\gamma_{m} - \\gamma_{p}} - \\frac{1}{2} = 0$$\n", "\n", "Unfortunately this equation cannont be solve for $t$ analytically, thus, we have to turn to a numerics using the ```fsolve``` function from scipy." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "deletable": true, "editable": true }, "outputs": [], "source": [ "# function that calculates the RHS of the above equation for a given set of gamma_m & gamma_p & t\n", "def response_time(t,decay):\n", " return (decay[0]*sp.exp(-decay[1]*t) - decay[1]*sp.exp(-decay[0]*t))/(decay[0] - decay[1]) - 0.5" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "# we want to find the value t for a given pair gamma_m & gamma_p that solves the above equation, \n", "# i.e. the RHS equals 0\n", "\n", "decay = [0.02,0.01]\n", "resp = fsolve(response_time,1,args = (decay,))\n", "print(resp)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "# define sets of decay rates and estimate the protein response time using a nested for loop\n", "steps = 15\n", "GAMMA_m = sp.logspace(-3,0,steps)\n", "GAMMA_p = sp.logspace(-3,0,steps)\n", "\n", "initial_guess = 1\n", "res = np.zeros((len(GAMMA_m),len(GAMMA_p)))\n", "for ii,gamma_m in enumerate(GAMMA_m):\n", " for jj,gamma_p in enumerate(GAMMA_p):\n", " decay = [gamma_m,gamma_p]\n", " t_12 = fsolve(response_time,initial_guess,args = (decay,))\n", " res[ii,jj] = t_12[0]\n", "\n", "res = sp.log10(res)\n", " \n", "# visualise the results as a colormap\n", "fig,ax = pl.subplots(1,1, figsize = (8,8))\n", "cax = ax.imshow(res)\n", "cbar = fig.add_axes([0.9, 0.123, 0.04, 0.755])\n", "fig.colorbar(cax, cax = cbar, orientation = \"vertical\")\n", "cbar.set_ylabel('$log_{10}(t_{1/2})$')\n", "ticks = np.arange(0,steps,2)\n", "ax.set_yticks(ticks)\n", "ax.set_xticks(ticks)\n", "ax.set_yticklabels(np.round(GAMMA_m[ticks],3))\n", "ax.set_xticklabels(np.round(GAMMA_p[ticks],3))\n", "ax.set_ylabel('RNA decay')\n", "ax.set_xlabel('Protein decay')" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "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 $\\beta = \\beta(t)$. We assume the functional relation for $\\beta$ to be:\n", "\n", "$$\\beta(t) = \\beta\\cdot e^{-k_1 t}\\cdot (1 - e^{-k_{2}t})$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "# functions\n", "def TF_binding(param,tt):\n", " k_1,k_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", " \n", " dRNA = param[0]*TF_binding(param_binding,t) - param[1]*RNA\n", " \n", " return dRNA" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "# visualize the TF binding pattern\n", "k_1 = 1.5\n", "k_2 = 0.1\n", "param_binding = [k_1,k_2]\n", "tt = sp.linspace(0,12,1000)\n", "binding = TF_binding(param_binding,tt)\n", "binding = binding/binding.max()\n", "\n", "fig,ax = pl.subplots(1,1, figsize = (6,4))\n", "ax.plot(tt,binding, lw = line_width, color = color_phase_space)\n", "ax.set_ylabel('TF binding')\n", "ax.set_xlabel('Time (h)')" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "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", "\n", "as in the example from the paper. Use the example given below to create new simulations and add the results to the plot." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "# define the binding pattern\n", "k_1 = 1.5\n", "k_2 = 0.1\n", "param_binding = [k_1,k_2]\n", "tt = sp.linspace(0,5,1000)\n", "\n", "binding = TF_binding(param_binding,tt)\n", "binding = binding/binding.max()\n", "\n", "y_0 = 0\n", "\n", "# global transcription rate\n", "beta = 1\n", "\n", "#-----------------------------\n", "# fast RNA decay\n", "gamma_m = 10\n", "param = [beta,gamma_m]\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", "\n", "# plot the result\n", "fig,ax = pl.subplots(1,1, figsize = (17,6))\n", "\n", "ax.plot(tt,sol_fast[:,0], label = 'RNA, fast decay', 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 = 'lower right')\n", "ax.set_xlabel('Time (h)')\n", "ax.set_ylabel('Concentration')\n", "ax.set_title('Time course')\n", "# ax[0].set_ylim(0,steady_state_protein+steady_state_protein*0.1)\n" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "## Correlation between RNA and protein abundance depends on the RNA decay rate\n", "\n", "**Exercise 4**\n", "Investigate the correlation between RNA and protein levels. \n", "- Plot RNA and protein levels agains each other\n", "- You can use simulations you have created earlier for a step like increase in RNA synthesis (Exercise 2)\n", "\n", "**Questions**\n", "- Under what conditions RNA and protein are strongly correlated" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "deletable": true, "editable": 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", "\n", "ax.set_ylabel('Protein')\n", "ax.set_xlabel('RNA')\n", "ax.legend(loc = 'best')" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "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", "On the example of the function implementing the simple transcription-translation ODE system try to write your own function to numerically solve the three ODEs resembling a circadian oscillator. \n", "\n", "- Name it 'circadian_ODEs'\n", "- it needs three input values: \n", " - an array _x_ for the actual values of the RNA and proteins\n", " - the time vector _t_ \n", " - a vector _param_ with the kinetic parameters. \n", "- You can solve the ODE system as before by calling the ```odeint``` function, you have to hand over the initial conditions, the predefined time vector and the parameter vector." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "deletable": true, "editable": 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", " \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": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "# define parameters\n", "k_1 = 1\n", "k_2 = 0.1\n", "k_3 = 1\n", "k_4 = 0.1\n", "k_5 = 1\n", "k_6 = 0.1\n", "K_i = 1\n", "n = 10\n", "\n", "param = [k_1,k_2,k_3,k_4,k_5,k_6,K_i,n]\n", "\n", "# define an intital state\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", "# use odeint and your function to generate simulations \n", "sol_circ = odeint(circadian_ODEs,y_0,tt,args = (param,))\n", "\n", "fig,ax = pl.subplots(1,1, figsize = (10,6))\n", "\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": { "deletable": true, "editable": true }, "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": { "deletable": true, "editable": true }, "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, "deletable": true, "editable": 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": false, "deletable": true, "editable": 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": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "# define parameters and time vector to solve the ODEs\n", "beta = 2\n", "gamma_m = 0.15\n", "k = 0.2\n", "gamma_p = 0.1\n", "param = [beta,gamma_m,k,gamma_p]\n", "\n", "# calculate steady state\n", "steady_state_RNA = beta/gamma_m\n", "steady_state_protein = k*beta/gamma_m/gamma_p\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", "# determinisitc 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, "deletable": true, "editable": 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", "- In addition, visualise the variability of the protein level by plotting a histogram of it.\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": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "# high transcription rate\n", "beta_high = 1\n", "\n", "# low translation rate\n", "k_low = 0.1\n", "\n", "# decay rates\n", "gamma_m = 0.15\n", "gamma_p = 0.1\n", "\n", "param_beta_high = [beta_high,gamma_m,k_low,gamma_p]\n", "\n", "# simulation time\n", "tf = 5000\n", "\n", "# set the initial state to the steady state\n", "initial_state_beta_high = [0,0]\n", "\n", "# run stochastic simulation\n", "sim_beta_high = stoch_sim_transcription_translation(param_beta_high,initial_state_beta_high,tf)\n", "\n", "# calculate coefficient of variation\n", "CV_beta_high = (sim_beta_high[:,2].std()/sim_beta_high[:,2].mean())\n", "\n", "fig,ax = pl.subplots(1,2, figsize = (17,5))\n", "\n", "ax[0].plot(sim_beta_high[:,0],sim_beta_high[:,2],lw = line_width,color = color_protein, \n", " label = 'High transcription, low translation, $CV$ = '+str(np.round(CV_beta_high,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_beta_high[:,2], color = color_protein);\n", "ax[1].set_ylabel('Frequency')\n", "ax[1].set_xlabel('Protein')" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true, "deletable": true, "editable": 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": false, "deletable": true, "editable": 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", "GAMMA_m = np.array([0.01,0.05,0.1,0.5])\n", "beta = steady_state_RNA * GAMMA_m\n", "\n", "gamma_p = k * steady_state_RNA / steady_state_protein\n", "\n", "k = 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 gamma_m\n", "res = []\n", "for ii,gamma_m in enumerate(GAMMA_m):\n", " param = [beta[ii],gamma_m,k,gamma_p]\n", " sim = stoch_sim_transcription_translation(param,initial_state,tf)\n", " \n", " res.append(sim)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "deletable": true, "editable": 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(GAMMA_m[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(GAMMA_m[ii]))\n", " \n", "ax[ii].set_xlabel('Time')\n", " \n", "pl.tight_layout()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "deletable": true, "editable": true }, "outputs": [], "source": [] } ], "metadata": { "hide_input": false, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.0" } }, "nbformat": 4, "nbformat_minor": 2 }