notebooks/Theory - Introduction to Brownian Motion simulation.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:984b804c3a73ab80cfe5540492cecafd5d4d1ebce6300cc8ce294883273a6cdc"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## 1.1 The Weiner Process\n",
      "\n",
      "A [**Weiner process**](http://en.wikipedia.org/wiki/Wiener_process), $W(t)\\,$, is a continuos-time stocastic process. By definition, the value of the process at time $t$ is:\n",
      "\n",
      "$$W(t) \\sim \\mathcal{N}(0,t)$$\n",
      "\n",
      "where $\\mathcal{N}(0,t)$ is a Normally-distributed random variable (RV) with $\\mu$=0 and $\\sigma^2=t$. \n",
      "\n",
      "From the definition follows that $W(0)=0$.\n",
      "\n",
      "Also, for any time instant $t$ and time delay $\\tau$ >0 the following is true:\n",
      "\n",
      "$$W(t+\\tau)-W(t) \\sim \\mathcal{N}(0,\\tau)$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## 1.2 Brownian Motion\n",
      "\n",
      "In **Brownian motion** of a freely diffusing particle, the mean squared displacement of a particle $\\langle|\\vec{r}(t)-\\vec{r}(t+\\tau)|^2\\rangle$ is proportional to the time interval $\\tau$ according to\n",
      "\n",
      "$$\\langle|\\vec{r}(t)-\\vec{r}(t+\\tau)|^2\\rangle = 2 D N \\tau$$\n",
      "\n",
      "- $\\vec{r}(t)$ position at time $t$\n",
      "- $N$ number of dimensions ($N$=3 for 3D simulations)\n",
      "- $D$ diffusion coefficient\n",
      "- $\\tau$ time interval."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## 1.3 Brownian Motion as a Weiner Process\n",
      "\n",
      "\n",
      "When using a Weiner process to describe a Brownian motion we must set a physical link between \n",
      "the variance of the Weineer process and the diffusion coefficient.\n",
      "\n",
      "Remembering that \n",
      "\n",
      "$$k \\mathcal{N}(\\mu,\\sigma^2) = \\mathcal{N}(\\mu,k\\sigma^2)$$\n",
      "\n",
      "if we build a process in which \"dispacements\" are normally distributed with variance equal to $2DN\\tau$:\n",
      "\n",
      "$$W(t+\\tau)-W(t) \\sim \\mathcal{N}(0,2DN\\tau)$$\n",
      "\n",
      "than we are describing the Brownian motion of a particle with diffusion coefficient $D$. To simulate this process we must choose the times at which to evauate the position. For example we can sample the time at uniform intervals with step = $\\Delta t$.\n",
      "\n",
      "###How to choose the simulation step $\\Delta t$\n",
      "The choice of the step depends on which properties we want to simulate. For example, let assume we want to simulate a colloidal particle diffusing through a confocal excitation volume of lateral dimension $S_{spot}$. In order to gather significant information we want to sample the particle position may times during the average diffusion time. The average diffusion time can be estimated setting the standard deviation of the displacement ($W(t+\\tau)-W(t)\\quad$) equal to $S_{spot}$ and solving for $\\tau$\n",
      "\n",
      "$$S_{spot} = \\sqrt{2DN\\tau_{spot}} \\quad  \\rightarrow  \\quad \\tau_{spot} = \\frac{S_{spot}^2}{2ND}$$\n",
      "\n",
      "so we want our simulation step to be $<< \\tau_{spot}$.\n",
      "\n",
      "Although $\\tau_{spot}$ can be derived theorically from $D$ and from the knowledge of the PSF, we know that for typical biomolecules of few nanometers, diffusing through a diffraction limited exciation spot (of visible light), the diffusion\n",
      "time is of the order of 1ms. Therefore we can safely set the simulation step to 0.5-1\u03bcs."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###The diffusion coefficient $D$\n",
      "The diffusion coefficient $D$ is given by:\n",
      "$$ D = \\frac{k_B T}{3 \\pi \\eta d} $$\n",
      "\n",
      "- $k_B$   Boltzman constant\n",
      "- $T$   temperature in Kelvin\n",
      "- $\\eta$   viscosity (in SI units: Pa/s)\n",
      "- $d$   radius of the particle in meters\n",
      "\n",
      "> *See also [Theory - On Browniam motion and Diffusion coefficient](Theory - On Browniam motion and Diffusion coefficient.ipynb)*\n",
      "\n",
      "Note that the units of $D$ are $\\mathrm{m}^2/\\mathrm{s}\\quad$. Using some reasonable number we obtain:\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "d = 5e-9        # particle radius in meters\n",
      "eta = 1.0e-3    # viscosity of water in SI units (Pascal-seconds) at 293 K\n",
      "kB = 1.38e-23   # Boltzmann constant\n",
      "T = 293         # Temperature in degrees Kelvin\n",
      "\n",
      "D = kB*T/(3*np.pi*eta*d) # [m^2 / s]\n",
      "D"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 2,
       "text": [
        "8.580361291970262e-11"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Or expressing $D$ in $\\textrm{nm}^2 /\\mu s$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Du = D*(1e9)**2/(1e6) # [nm^2 / us]\n",
      "Du"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "85.80361291970262"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can also estimate $D$ experimentally from the knowledge of the PSF and the diffusion time $\\tau_{spot}$\n",
      "\n",
      "$$S_{spot} = \\sqrt{2DN\\tau_{spot}} \\quad  \\rightarrow  \\quad D = \\frac{S_{spot}^2}{2N\\tau_{spot}}$$\n",
      "\n",
      "Putting some reasonable number we obtain:\n",
      "\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "S_spot = 0.8e-6\n",
      "N = 3\n",
      "tau_spot = 1e-3\n",
      "\n",
      "D = S_spot**2 / (2*N*tau_spot) # [m^2 / s]\n",
      "D"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "1.0666666666666665e-10"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Du = D*(1e6)**2/(1e3) # [um^2 / ms]\n",
      "Du"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "0.10666666666666665"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "not very different from what we obtained before from the viscosity model."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##Examples"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- *How far we travel in X seconds (hint standard deviation of displacement)?*"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "time = 10.  # seconds\n",
      "sigma = np.sqrt(2*D*3*time)\n",
      "print 'Displacement (std_dev): %.2f um' % (sigma*1e6)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Displacement (std_dev): 80.00 um\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- *How long we need to diffuse an X distance?*"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "space = 0.1e-3  # m\n",
      "time = 1.*space**2/(2*D*3)\n",
      "print 'Time for %.1f um displacement: %.1f s' % (space*1e6, time)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Time for 100.0 um displacement: 15.6 s\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    }
   ],
   "metadata": {}
  }
 ]
}