Tutorial 1#

Dynamical Systems II: Continuous Dynamical Systems & Differential Equations

[insert your name]

Important reminders: Before starting, click “File -> Save a copy in Drive”. Produce a pdf for submission by “File -> Print” and then choose “Save to PDF”.

To complete this tutorial, you should have watched Videos 4.1-4.4.

Imports

# @markdown Imports

# Imports
import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets  # interactive display

Plotting functions

# @markdown Plotting functions
import numpy
from numpy.linalg import inv, eig
from math import ceil
from matplotlib import pyplot, ticker, get_backend, rc
from mpl_toolkits.mplot3d import Axes3D
from itertools import cycle


%config InlineBackend.figure_format = 'retina'
plt.style.use("https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/nma.mplstyle")

classic = 'k'

def plot_derivative_vs_variable(x, dx):

  fig, ax = plt.subplots()

  ax.plot([-1, 1], [0, 0],'k')
  ax.plot(x, dx, 'b')
  ax.set(xlabel=r'P$_o$', ylabel=r'$\frac{dP_o}{dt}$')

def plot_trajectories(t, numerical_P, init_Po, alpha, beta):

  # Get analytical solution
  analytical_t = np.arange(0, t[-1], .01)
  C = init_Po - alpha/(alpha + beta)
  analytical_P = alpha/(alpha + beta) + C*np.exp(-(alpha+beta)*analytical_t)

  fig, ax = plt.subplots()
  ax.plot(analytical_t, analytical_P,'#33D396', linewidth=3)
  ax.plot(t, numerical_P, '#1F59EB')
  
  ax.set(xlabel='Time (ms)', ylabel = 'Probability of channel open', ylim=[-1, 1])
  plt.legend(['Analytical','Numerical'])
  

Key concept review & coding tips#

Dynamical systems#

System concerned with movement of point over time in some relevant geometrical space
The state of a dynamical system is a complete snapshot of the system at that point in time.
The state space is the relevant geometrical space in which to model the state evolving (the set of all possible states)
Are discrete or continuous depending on how the evolution rules are defined
Discrete dyanmical systems have step-like update rules: $x_{n+1} = F(x_n)$
Continuous dyanmical systems have differential equations governing the evolution: $\frac{dx}{dt} = F(x, t)$

Differential equations#

Equations that relate functions and their derivatives ($\frac{dx}{dt} = F(x, t)$)
Can sometimes compute the analytical solution $x(t)$ by integration and solving for constants using the initial condition. This is not always possible.
Can numerically approximate specific trajectories (corresponding to different initial conditions) using Euler’s method, which is based on the update rule $$ x(t+\Delta t) = F(x, t)\Delta t + x(t)$$
Can understand the system graphically by plotting a direction field (tangent lines on x vs t graph)
Can also understand the system graphically by plotting $\frac{dx}{dt}$ vs x and analyzing critical/equilibrium points
Critical/equilbrium points occur at x for which $\frac{dx}{dt} = 0$. They are either stable (attractors) or unstable (repellers) depening on the derivative sign to either side

Exercise 1: Defining the differential equation#

Note: please try this out first but ask for help if you get stuck! Don’t spend too long on this exercise - it’s a bit tricky and we’re happy to take you through it. Also, when you have a solution to part D, call us in to check it’s correct as the rest of the tutorial depends on it

Throughout this tutorial, we will be working on a specific dynamical system concerning ion channel kinetics.

Let’s say that we have a population of ion channels of the same type, which are individually either open or closed. There is some rate at which the channels open, dependent on the voltage of the cell, $\alpha(V)$ and some rate at which they close, $\beta(V)$. These would be in units 1/ms (essentially the number of times in a ms that a gate which is open closes and vice versa) and can be found empirically. We suddenly step the voltage to one that opens channels and keep it at that voltage. Although $\alpha$ and $\beta$ are normally voltage-dependent, we will treat them as constants from here on out since we are holding the voltage steady.

We need to introduce two more variables before writing out a differential equation that defines our system: $P_o$ is the proportion of channels that are open (or equivalently the probability that a particular channel is open) and $P_c$ is the the proportion of closed channels (or equivalently the probability that a particular channel is closed).

We will now write out a differential equation to model what happens when we step to a certain voltage. In particular, we want to figure out how the proportion of open channels, $P_o$, changes over time. This means we want to find a differential equation $\frac{dP_o}{dt} = ...$. We will take you through this process step by step.

Fun fact: this system is very important for the Hodgkin-Huxley model, although we need to add in voltage dependence for that.

A) Rate of change of the proportion of open channels#

Let’s say $P(opening)$ is the proportion of channels opening over some tiny time step and $P(closing)$ is the proportion of channels closing.

How does the rate of change of the proportion of open channels ($\frac{dP_o}{dt}$) relate to $P(opening)$ and $P(closing)$?

Answer#

Fill in your math answer here (if in latex) so it is green. Use

for line breaks

B) Proportion of channels opening#

Write out an equation for the proportion of channels opening over some small time step. (Assume we’ll later multiply with dt so no need to write a small time step here)

\[ P(opening) = ...\]

Hints: Note that $P(opening)$ is different from $P(open)$. What state does the channel need to be in at the start of the time step for it to open? Relevant variables might be $\alpha$, $\beta$, $P_o$, and $P_c$.

Answer#

Fill in your math answer here (if in latex) so it is green. Use

for line breaks

C) Proportion of channels closing#

Write an equation for the $P(closing)$. Should resemble that in part A (essentially a sort of inverse of part A).

\[P(closing) = \]

Answer#

Fill in your math answer here (if in latex) so it is green

D) Getting dP_o/dt#

Substitute what you equated $P(opening)$ and $P(closing)$ to in part B/C in the equation for $\frac{dP_o}{dt}$ from part A.

Use the fact that the proportion of closed channels equals 1 minus the proportion of open channels ($P_c = 1 - P_o$) to further simplify. Write out the equation as simply as possible.

Your equation should be a first order linear differential equation.

Answer#

Fill in your math answer here (if in latex) so it is green. Use

for line breaks

Exercise 2: Graphical solution: assessing critical points and stability#

We’ll now try and understand what’s happening in our system through graphical methods.

A) Plotting derivative vs variable#

We will first plot a graph of $\frac{dP_o}{dt}$ vs $P_o$ to determine critical points. So far, we have used general variables $\alpha$ and $\beta$ for representing the rates of opening and closing, but let’s define them numerically now:

\[\begin{split} \alpha = 0.3 \\ \beta = 0.1 \end{split}\]

Note that below we’re allowing $P_o$ to be negative - this doesn’t really make sense as it is a probability but we want to fully understand the associated dynamical system.

# Set parameters
alpha = .3
beta = .1

# We will evalute dP_o/dt on a span of possible P_o values
Po = np.arange(-1, 1, .1)

# Fill in your code here to compute dPo/dt for all Po above
dPo_dt = ...

# Uncomment below to use plotting function plot_derivative_vs_variable
# plot_derivative_vs_variable(Po, dPo_dt)

B) Identifying equilibrium points and stability#

i) What is the critical/equilibrium point for this system with these values of $\alpha$ and $\beta$? Why?

ii) What is the stability of this critical point? Why?

iii) What does the equilibrium point correspond to in this system? (aka bring it back to the neuro)

Answer#

Fill in your text answer here so it is green. Use

for line breaks

C) Describing potential trajectories#

Based on the answers in part B, let’s predict what some trajectories in this system might look like. Be specific - will it rise or fall and where will it end up after a long time. You are welcome to describe these in words or sketch out a labeled graph.

i) What will the trajectory look like if the initial $P_o$ is negative?

ii) What will the trajectory look like if the initial $P_o$ is above 0 but below the critical point?

iii) What will the trajectory look like if the initial $P_o$ is above the critical point?

Answer#

Fill in your text answer here so it is green. Use

for line breaks

(Optional) D) Getting the critical point for any value of alpha and beta#

We used specific values of $\alpha$ an $\beta$ and got a numerical critical point.

i) We could write the critical point in terms of $\alpha$ and $\beta$ though - what would it be?

ii) If the rate of opening was faster, would the critical point be higher or lower? Does this make sense intuitively?

Answer#

Fill in your text/math answer here so it is green. Use

for line breaks

Exercise 3: Numerical solution: implementing Euler’s method and examining trajectories#

We’ll now look at some specific trajectories by numerically approximation the system using Euler’s method. This sytem does have an analytical solution so we will compare the accuracy of our numerical approximation to the analytical solution. If the numerical approximation is good enough, it can be hard to see the analytical solution line, although it is a bit wider to help. If you are interested in finding out how to compute this analytical solution try Exercise 4!

A) Complete the function below to integrate our differential equation numerically#

def prob_open_euler_integration(init_Po, dt, T, alpha=0.3, beta=0.1):
  """Numerically approximate the solution of the differential equation governing 
  ion channel kinetics given an initial condition for a duration T. 

  Args:
    init_Po (scalar): value of Po at start of simulation (at time = 0)
    dt (scalar): timestep of the simulation
    T (scalar): total duration of the simulation
    alpha (scalar): rate of change of ion channel opening
    beta (scalar): rate of change of ion channel closing

  Returns:
    ndarray, ndarray: Po and the time t at each step
  """

  # Initialize variables
  t = np.arange(0, T, dt)
  P = np.zeros_like(t)
  P[0] = init_Po # This is Po at time 0

  # Step through time
  for n in range(0, len(t)-1):

    # for each point in time, compute dP_dt from P[n]
    dP_dt = ...

    # Update P[n+1] based on P[n] and dP_dt
    P[n+1] = ...

  return P, t


init_Po = .4
dt = .01 # step size is .01 ms
T = 50 # run simulation for 50 ms
alpha = .3
beta = .1

numerical_P, t = prob_open_euler_integration(init_Po, dt, T, alpha, beta)

plot_trajectories(t, numerical_P, init_Po, alpha, beta)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[4], line 39
     36 alpha = .3
     37 beta = .1
---> 39 numerical_P, t = prob_open_euler_integration(init_Po, dt, T, alpha, beta)
     41 plot_trajectories(t, numerical_P, init_Po, alpha, beta)

Cell In[4], line 28, in prob_open_euler_integration(init_Po, dt, T, alpha, beta)
     25   dP_dt = ...
     27   # Update P[n+1] based on P[n] and dP_dt
---> 28   P[n+1] = ...
     30 return P, t

TypeError: float() argument must be a string or a number, not 'ellipsis'

B) Checking our intuition#

Use this demo to look at the trajectory given different initial conditions. Try out all three types outlined in 2C and verify your answers there were correct.

You don’t need any answer here, just try out the demo.

Execute this cell for the demo

# @markdown Execute this cell for the demo

@widgets.interact
def plot_euler_integration(init_Po = (-1, 1, .1) ):
  dt = .01 # step size is .01 ms
  T = 50 # run simulation for 50 ms
  alpha = .3
  beta = .1

  numerical_P, t = prob_open_euler_integration(init_Po, dt, T, alpha, beta)


  plot_trajectories(t, numerical_P, init_Po, alpha, beta)

C) Euler integration accuracy#

Use the interactive demo below to assess the accuracy of the numerical approximation for different step sizes (values of dt), then answer the questions below.

Execute this cell for the demo

# @markdown Execute this cell for the demo

@widgets.interact
def plot_euler_integration(dt = widgets.SelectionSlider(options=[("%g"%i,i) for i in np.arange(0.001, 5.001, 0.01)])):
  # Have to do this clunky word around to show small values in slider accurately
  # (from https://github.com/jupyter-widgets/ipywidgets/issues/259)
  init_Po = .4
  T = 50 # run simulation for 50 ms
  alpha = .3
  beta = .1

  numerical_P, t = prob_open_euler_integration(init_Po, dt, T, alpha, beta)

  plot_trajectories(t, numerical_P, init_Po, alpha, beta)

i) Does the approximation get better or worse with larger steps?

ii) Why?

Answer#

Fill in your text answer here so it is green. Use

for line breaks

(Optional) D) Playing around with alpha/beta#

Change alpha and beta and see how the equilibrium solution and speed at which the system converges changes

Execute this cell for the demo

# @markdown Execute this cell for the demo

@widgets.interact
def plot_euler_integration(alpha = (0.01, 1, .01), beta = (0.01, 1, .01)):
  # Have to do this clunky word around to show small values in slider accurately
  # (from https://github.com/jupyter-widgets/ipywidgets/issues/259)
  init_Po = .5
  T = 50 # run simulation for 50 ms
  dt = 0.01

  numerical_P, t = prob_open_euler_integration(init_Po, dt, T, alpha, beta)

  plot_trajectories(t, numerical_P, init_Po, alpha, beta)

(Optional) Exercise 4: Analytical solution: integrating the differential equation#

This differential equation does have an analytical solution that we are able to compute without crazy integration.

Remember that the general form of a first order linear differential equation is:

\[ \frac{dx}{dt} + p(t)x = g(t)\]

In this case, we can compute the solution as:

\[ x(t) = \frac{\int \mu(t) g(t) dt + C}{\mu(t)} \]

where $\mu(t) = e^{\int p(t) dt}$ Go forth and try it out!

Answer#

Fill in your text answer here so it is green. Use

for line breaks

Mathematical Tools for Neuroscientists

Tutorial 1

Contents

Tutorial 1#

Key concept review & coding tips#

Dynamical systems#

Differential equations#

Exercise 1: Defining the differential equation#

A) Rate of change of the proportion of open channels#

Answer#

B) Proportion of channels opening#

Answer#

C) Proportion of channels closing#

Answer#

D) Getting dP_o/dt#

Answer#

Exercise 2: Graphical solution: assessing critical points and stability#

A) Plotting derivative vs variable#

B) Identifying equilibrium points and stability#

Answer#

C) Describing potential trajectories#

Answer#

(Optional) D) Getting the critical point for any value of alpha and beta#

Answer#

Exercise 3: Numerical solution: implementing Euler’s method and examining trajectories#

A) Complete the function below to integrate our differential equation numerically#

B) Checking our intuition#

C) Euler integration accuracy#

Answer#

(Optional) D) Playing around with alpha/beta#

(Optional) Exercise 4: Analytical solution: integrating the differential equation#

Answer#