\n",
"## Click here for text recap of video

\n",
"\n",
"We can run this simulation many times and gather empirical distributions of open/closed states. Alternatively, we can formulate the exact same system probabilistically, keeping track of the probability of being in each state.\n",
"\n",
"\n",
"\n",
"(see diagram in lecture)\n",
"\n",
"The same system of transitions can then be formulated using a vector of 2 elements as the state vector and a dynamics matrix $\\mathbf{A}$. The result of this formulation is a *state transition matrix*:\n",
"\n",
"\\begin{equation}\n",
"\\left[ \\begin{array}{c} C \\\\ O \\end{array} \\right]_{k+1} = \\mathbf{A} \\left[ \\begin{array}{c} C \\\\ O \\end{array} \\right]_k = \\left[ \\begin{array} & 1-\\mu_{\\text{c2o}} & \\mu_{\\text{o2c}} \\\\ \\mu_{\\text{c2o}} & 1-\\mu_{\\text{o2c}} \\end{array} \\right] \\left[ \\begin{array}{c} C \\\\ O \\end{array} \\right]_k.\n",
"\\end{equation}\n",
"\n",
"Each transition probability shown in the matrix is as follows:\n",
"1. $1-\\mu_{\\text{c2o}}$, the probability that the closed state remains closed. \n",
"2. $\\mu_{\\text{c2o}}$, the probability that the closed state transitions to the open state.\n",
"3. $\\mu_{\\text{o2c}}$, the probability that the open state transitions to the closed state. \n",
"4. $1-\\mu_{\\text{o2c}}$, the probability that the open state remains open. \n",
"\n",
"

\n",
"\n",
"_Notice_ that this system is written as a discrete step in time, and $\\mathbf{A}$ describes the transition, mapping the state from step $k$ to step $k+1$. This is different from what we did in the exercises above where $\\mathbf{A}$ had described the function from the state to the time derivative of the state."
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"## Coding Exercise 2: Probability Propagation\n",
"\n",
"*Referred to in video as exercise 2B*\n",
"\n",
"Complete the code below to simulate the propagation of probabilities of closed/open of the ion channel through time. A variable called `x_kp1` (short for, $x$ at timestep $k$ plus 1) should be calculated per each step *k* in the loop. However, you should plot $x$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {}
},
"outputs": [],
"source": [
"def simulate_prob_prop(A, x0, dt, T):\n",
" \"\"\" Simulate the propagation of probabilities given the transition matrix A,\n",
" with initial state x0, for a duration of T at timestep dt.\n",
"\n",
" Args:\n",
" A (ndarray): state transition matrix\n",
" x0 (ndarray): state probabilities at time 0\n",
" dt (scalar): timestep of the simulation\n",
" T (scalar): total duration of the simulation\n",
"\n",
" Returns:\n",
" ndarray, ndarray: `x` for all simulation steps and the time `t` at each step\n",
" \"\"\"\n",
"\n",
" # Initialize variables\n",
" t = np.arange(0, T, dt)\n",
" x = x0 # x at time t_0\n",
"\n",
" # Step through the system in time\n",
" for k in range(len(t)-1):\n",
" ###################################################################\n",
" ## TODO: Insert your code here to compute x_kp1 (x at k plus 1)\n",
" raise NotImplementedError(\"Student exercise: need to implement simulation\")\n",
" ## hint: use np.dot(a, b) function to compute the dot product\n",
" ## of the transition matrix A and the last state in x\n",
" ## hint 2: use np.vstack to append the latest state to x\n",
" ###################################################################\n",
"\n",
" # Compute the state of x at time k+1\n",
" x_kp1 = ...\n",
"\n",
" # Stack (append) this new state onto x to keep track of x through time steps\n",
" x = ...\n",
"\n",
" return x, t\n",
"\n",
"\n",
"# Set parameters\n",
"T = 500 # total Time duration\n",
"dt = 0.1 # timestep of our simulation\n",
"\n",
"# same parameters as above\n",
"# c2o: closed to open rate\n",
"# o2c: open to closed rate\n",
"c2o = 0.02\n",
"o2c = 0.1\n",
"A = np.array([[1 - c2o*dt, o2c*dt],\n",
" [c2o*dt, 1 - o2c*dt]])\n",
"\n",
"# Initial condition: start as Closed\n",
"x0 = np.array([[1, 0]])\n",
"\n",
"# Simulate probabilities propagation\n",
"x, t = simulate_prob_prop(A, x0, dt, T)\n",
"\n",
"# Visualize\n",
"plot_state_probabilities(t,x)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"execution": {}
},
"source": [
"[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/main//tutorials/W2D2_LinearSystems/solutions/W2D2_Tutorial2_Solution_c5188a2d.py)\n",
"\n",
"*Example output:*\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"Here, we simulated the propagation of probabilities of the ion channel's state changing through time. Using this method is useful in that we can **run the simulation once** and see **how the probabilities propagate throughout time**, rather than re-running and empirically observing the telegraph simulation over and over again. \n",
"\n",
"Although the system started initially in the Closed ($x=0$) state, over time, it settles into a equilibrium distribution where we can predict what fraction of time it is Open as a function of the $\\mu$ parameters. We can say that the plot above shows this _relaxation towards equilibrium_.\n",
"\n",
"Re-calculating our value of the probability of $c2o$ again with this method, we see that this matches the simulation output from the telegraph process! \n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {}
},
"outputs": [],
"source": [
"print(\"Probability of state c2o: %.3f\"%(c2o / (c2o + o2c)))\n",
"x[-1,:]"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"---\n",
"# Section 3: Equilibrium of the telegraph process\n",
"\n",
"*Estimated timing to here from start of tutorial: 30 min*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Video 3: Continuous vs. Discrete Time Formulation\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {},
"tags": [
"remove-input"
]
},
"outputs": [],
"source": [
"# @title Video 3: Continuous vs. Discrete Time Formulation\n",
"from ipywidgets import widgets\n",
"from IPython.display import YouTubeVideo\n",
"from IPython.display import IFrame\n",
"from IPython.display import display\n",
"\n",
"\n",
"class PlayVideo(IFrame):\n",
" def __init__(self, id, source, page=1, width=400, height=300, **kwargs):\n",
" self.id = id\n",
" if source == 'Bilibili':\n",
" src = f'https://player.bilibili.com/player.html?bvid={id}&page={page}'\n",
" elif source == 'Osf':\n",
" src = f'https://mfr.ca-1.osf.io/render?url=https://osf.io/download/{id}/?direct%26mode=render'\n",
" super(PlayVideo, self).__init__(src, width, height, **kwargs)\n",
"\n",
"\n",
"def display_videos(video_ids, W=400, H=300, fs=1):\n",
" tab_contents = []\n",
" for i, video_id in enumerate(video_ids):\n",
" out = widgets.Output()\n",
" with out:\n",
" if video_ids[i][0] == 'Youtube':\n",
" video = YouTubeVideo(id=video_ids[i][1], width=W,\n",
" height=H, fs=fs, rel=0)\n",
" print(f'Video available at https://youtube.com/watch?v={video.id}')\n",
" else:\n",
" video = PlayVideo(id=video_ids[i][1], source=video_ids[i][0], width=W,\n",
" height=H, fs=fs, autoplay=False)\n",
" if video_ids[i][0] == 'Bilibili':\n",
" print(f'Video available at https://www.bilibili.com/video/{video.id}')\n",
" elif video_ids[i][0] == 'Osf':\n",
" print(f'Video available at https://osf.io/{video.id}') \n",
" display(video)\n",
" tab_contents.append(out)\n",
" return tab_contents\n",
"\n",
"\n",
"video_ids = [('Youtube', 'csetTTauIh8'), ('Bilibili', 'BV1di4y1g7Yc')]\n",
"tab_contents = display_videos(video_ids, W=730, H=410)\n",
"tabs = widgets.Tab()\n",
"tabs.children = tab_contents\n",
"for i in range(len(tab_contents)):\n",
" tabs.set_title(i, video_ids[i][0])\n",
"display(tabs)"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"Since we have now modeled the propagation of probabilities by the transition matrix $\\mathbf{A}$ in Section 2, let's connect the behavior of the system at equilibrium with the eigendecomposition of $\\mathbf{A}$.\n",
"\n",
"As introduced in the lecture video, the eigenvalues of $\\mathbf{A}$ tell us about the stability of the system, specifically in the directions of the corresponding eigenvectors."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {}
},
"outputs": [],
"source": [
"# compute the eigendecomposition of A\n",
"lam, v = np.linalg.eig(A)\n",
"\n",
"# print the 2 eigenvalues\n",
"print(\"Eigenvalues:\",lam)\n",
"\n",
"# print the 2 eigenvectors\n",
"eigenvector1 = v[:,0]\n",
"eigenvector2 = v[:,1]\n",
"print(\"Eigenvector 1:\", eigenvector1)\n",
"print(\"Eigenvector 2:\", eigenvector2)"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"## Think! 3: Finding a stable state\n",
"\n",
"1. Which of these eigenvalues corresponds to the **stable** (equilibrium) solution? \n",
"2. What is the eigenvector of this eigenvalue? \n",
"3. How does that explain the equilibrium solutions in simulation in Section 2 of this tutorial?\n",
"\n",
"_hint_: our simulation is written in terms of probabilities, so they must sum to 1. Therefore, you may also want to rescale the elements of the eigenvector such that they also sum to 1. These can then be directly compared with the probabilities of the states in the simulation. "
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"execution": {}
},
"source": [
"[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/main//tutorials/W2D2_LinearSystems/solutions/W2D2_Tutorial2_Solution_37abbdad.py)\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"---\n",
"# Summary\n",
"\n",
"*Estimated timing of tutorial: 45 minutes*\n",
"\n",
"\n",
"In this tutorial, we learned:\n",
"\n",
"* The definition of a Markov process with history dependence.\n",
"* The behavior of a simple 2-state Markov process -- the telegraph process--can be simulated either as a state-change simulation or as a propagation of probability distributions.\n",
"* The relationship between the stability analysis of a dynamical system expressed either in continuous or discrete time.\n",
"* The equilibrium behavior of a telegraph process is predictable and can be understood using the same strategy as for deterministic systems in Tutorial 1: by taking the eigendecomposition of the A matrix."
]
}
],
"metadata": {
"colab": {
"collapsed_sections": [],
"include_colab_link": true,
"name": "W2D2_Tutorial2",
"provenance": [],
"toc_visible": true
},
"kernel": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"kernelspec": {
"display_name": "Python 3",
"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.7.13"
}
},
"nbformat": 4,
"nbformat_minor": 0
}