Research: Practice on Matlab Day 3 (MIT 6.094)

Wednesday was my third day of practice, I finished all the code but couldn’t get a solution on the problem regarding Hodgkin-Huxley differential equations using Matlab ode45 solver. It was so slow. Well I did get a solution for timescale 0 to 0.5 millisecond. But as I increased it to 20 ms, I left it overnight and when I woke up it still ran!  In fact it went up to about 70.5 million matrix element, that couldn’t be right. So I researched about Hodgkin-Huxley and found that a paper (Chance & Clore, 1981)  on it describes the HH equation partial differential equation as a stiff differential equations. Which means there are terms that can lead to rapid variation and operates on large and small timescales. With the usual numerical methods the solution set will become unstable. Or solved maybe very slowly.

I thought this was the problem, so maybe I should use ode23s or  ode15s solver, but because I don’t want to just brute force it and leave it overnight again I went to campus to ask Dr. Kyle Webb on my problem. So he did agree that if it was stiff I should try and use the ode15s stiff solver but since the question explicitly says I am recommended to use ode45 it is far more likely that I made a mistake in my implementation.

So he did what I should’ve done in the beginning. He checked all my equations and compared with the problem. Lo and behold, I forgot a negative sign on one of the equations.

Untitled.png

The lecture (link) exercises and homework (link) associated with Wednesday’s practice is posted below

Lectures

polyFittingLecture.m – this illustrates polynomial fitting, representation, and evaluation

MIT6094-L3-1polyfit.png

%Matlab Lecture 3
clc;
close all;

%y=x^2
p=[1,0,0];
x=-4:0.1:4;
%evaluate polynomial
y=polyval(p,x);

%add noise
yNoisy=y+rand(size(x));
%plot noisy data
plot(x,yNoisy,'.');
%curve fit the noisy data
pFit=polyfit(x,yNoisy,2);
%evaluate fitted polynomial
yFit=polyval(pFit,x);
hold on;
plot(x,yFit,'r-');

chemode.m – on setting up ordinary differential equations (ODE)

%Matlab Lecture 3 - ODE
function dydt=chemode(t,y)
    dydt=zeros(2,1);
    dydt(1)=-10*y(1)+50*y(2);
    dydt(2)=10*y(1)-50*y(2);
end

solveChemode.m – on solving ODEs

solveChemOde_01.png

%Matlab Lecture 3 - ODE Solver
clc;
close all;
%solve the chemode equations from t=0 to 0.5 and
%initial value y(1)=0 and y(2)=1
[t,y]=ode45('chemode',[0,0.5],[0,1]);
%plotting the amount of y(1) or substance A over time
plot (t,y(:,1),'k','LineWidth',1.5); hold on;
%plotting the amount of y(2) or substance B over time
plot (t,y(:,2),'r','LineWidth',1.5);
%formatting
legend('A','B');
xlabel('Time (s)');
ylabel('Amount of chemical (g)');
title('Chem reaction');

pendulum.m – tricks to deal with higher order ODEs

basically rephrase the higher order term as multiple variables and equations with only single order terms

%Matlab Lecture 3 - higher order ODE

%d(dtheta) + g/L sin(theta) = 0
%has to be broken down to a system of first order equations

function dxdt = pendulum(t,x)
    %0. Constant
    L=1;        %length of pendulum
    g=9.81;     %gravity

    %1. Input setup
    theta=x(1); %the angle
    gamma=x(2); %derivative of the angle, defined below

    %2. ODE
    dtheta=gamma;
    dgamma=-(g/L)*sin(theta);

    %3. The derivatives setup
    dxdt=zeros(2,1);    %derivatives initialization
    dxdt(1)=dtheta;
    dxdt(2)=dgamma;

end

Homework

Problem 1 – Linear System of Equations

% Homework 3 Problem 1
% Linear system of equations. Solve the following system of equations using \. Compute and
% display the error vector
% 3a + 6b + 4c = 1
% a + 5b = 2
% 7b + 7c = 3

%rephrase knowns
A=[3,6,4;1,5,0;0,7,7];
b=[1;2;3];

%check the rank
rankProblem=rank(A);
disp('rankProblem');
disp(rankProblem);

%Ax=b -> x=A\b
x=A\b;
disp('x');disp(x);

%Error vector abs(b-Ax)
Error = abs(b-A*x);
disp('Error');disp(Error);

Problem 2 – Numerical Integration Using Trapezoids and Adaptive Quadrature

MIT6094-H3-2integral.png

% Homework 3 Problem 2
clc;close all;
% what is the value of int_0^5 xe^(-x/3)dx
% using adaptive quadrature
quadInt=integral(@(x) x.*exp(-x./3),0,5);
% using trapezoidal rule
x=linspace(0,5,100);
y=x.*exp(-x./3);
trapzInt=trapz(x,y);
figure;plot(x,y,'.-');
% compute difference with -24e^(-5/3)+9
quadError=quadInt-(-24*exp(-5/3)+9);
trapzError=trapzInt-(-24*exp(-5/3)+9);
%display results
disp(['quadInt=' num2str(quadInt) ', trapzInt=' num2str(trapzInt)]);
disp(['quadError' num2str(quadError) ', trapzError=' num2str(trapzError)]);

Problem 3 – Inverse

% Homework 3 Problem 3
% Computing the inverse. Calculate the inverse of A and verify that when you multiply the
% original matrix by the inverse, you get the identity matrix (inv). Display the inverse matrix as well
% as the result of the multiplication of the original matrix by its inverse

A=[1,2;3,4];
invA=inv(A)
I=A*invA

Problem 4 – Polynomial Fit

MIT6094-H3-4polynomialfit.png

% Homework 3 Problem 4
% Fitting polynomials. Write a script to load the data file randomData.mat (which contains
    % variables x and y) and
clc; clear all;
load randomData.mat;
% fit first, second, third, fourth, and fifth degree polynomials to it.
[p1,s1,mu1]=polyfit(x,y,1);
[p2,s2,mu2]=polyfit(x,y,2);
[p3,s3,mu3]=polyfit(x,y,3);
[p4,s4,mu4]=polyfit(x,y,4);
[p5,s5,mu5]=polyfit(x,y,5);
%evaluate fitted polynomial
[y1,delta1] = polyval(p1,x,s1,mu1);
[y2,delta2] = polyval(p2,x,s2,mu2);
[y3,delta3] = polyval(p3,x,s3,mu3);
[y4,delta4] = polyval(p4,x,s4,mu4);
[y5,delta5] = polyval(p5,x,s5,mu5);

% Plot the data as blue dots on a figure, and plot all five polynomial fits using lines of different colors on
    % the same axes.
close all;
figure;
plot(x,y,'k.',x,y1,'c',x,y2,'b',x,y3,'m',x,y4,'r',x,y5,'y');
% Label the figure appropriately.
xlabel('x');
ylabel('y');
title('polynomial fit');
legend('data','p1','p2','p3','p4','p5');
% To get good fits, you’ll have to use the centering and scaling version of polyfit
    % (the one that returns three arguments, see help) and its counterpart in polyval (the one that accepts
    % the centering and scaling parameters).

Problem 5 – ODEs – Hodgkin Huxley Model of Neurons*

This is is my favorite problem, because I got to learn some neat stuff about how neurons behaves, about differential equations, numerical methods associated with them, and numerical problems.

This is the ODE setup

% Homework 3 Problem 5 - Hodgkin Huxley ODE
function dydt=HHode(t,x)
    % Constants
    C = 1 ;
    GK = 36 ;
    GNa = 120 ;
    GL = 0.3 ;
    EK = -72 ;
    ENa = 55 ;
    EL = -49.4 ;

    %Input Setup
    n=x(1);
    m=x(2);
    h=x(3);
    V=x(4);

    % ODE
    dndt=(1-n)*alphan(V)-n*betan(V);
    dmdt=(1-m)*alpham(V)-m*betam(V);
    dhdt=(1-h)*alphah(V)-h*betah(V);
    dVdt=-1/C*(GK*n^4*(V-EK)+GNa*m^3*h*(V-ENa)+GL*(V-EL)); %K.Webb: fix negative sign mistake

    % Derivative Setup
    dydt=[dndt;dmdt;dhdt;dVdt];
end

This is the solver

% Homework 3 Problem 5 - Main
% Hodgkin-Huxley model of the neuron
clc; close all;
% You will write an ODE file to describe the spiking of a neuron, based on the equations developed by
    %Hodgkin and Huxley in 1952 (they received a Nobel Prize for this work)
% The main idea behind the model is that ion channels in the neuron’s membrane have voltage-sensitive gates
    %that open or close as the transmembrane voltage changes
% Once the gates are open, charged ions can flow through them, affecting the transmembrane voltage
% The equations are nonlinear and coupled, so they must be solved numerically

% a. Download the HH.zip file from the class website and unzip its contents into your homework folder
% This zip folder contains 6 m-files: alphah.m, alpham.m, alphan.m, betah.m, betam.m, betan.m
% These functions return the voltagedependent opening ? (V ) ) and closing ( ? (V ) ) rate constants
    % for the h, m, and n ?( V ) gates: ? ( V )

% b. Write an ODE file to return the following derivatives
% 	and the following constants ( C is membrane capacitance, G are the conductances and E are the reversal potentials of the
% 	potassium ( K ), sodium ( Na ), and leak ( L ) channels): C = 1 GK = 36 GNa = 120 GL = 0.3 EK = ?72 ENa = 55 EL = ?49.4

% c. Write a script called HH.m which will solve this system of ODEs and plot the transmembrane voltage
% First, we’ll run the system to steady state
% Run the simulation for 20ms (the timescale of the equations is ms) with initial values:
    %n = 0.5; m = 0.5;h = 0.5;V = ?60 (ode45)
[t,y]=ode45('HHode',[0,20],[0.5,0.5,0.5,-60]);

% Store the steady-state value of all 4 parameters in a vector called ySS
n=y(:,1);
m=y(:,2);
h=y(:,3);
V=y(:,4);
ySS=y(end,:);

% Make a new figure and plot the timecourse of V (t) to verify that it reaches steady state by the end of the simulation
figure(1);
plot(t,V);
xlabel('t [ms]');ylabel('Transmembrane Voltage [mV]');title('Approaching Steady State');
% d. Next, we’ll explore the trademark feature of the system: the all-or-none action potential
% Neurons are known to ‘fire’ only when their membrane surpasses a certain voltage threshold
% To find the threshold of the system, solve the system 10 times, each time using ySS as ...
    %the initial condition but increasing the initial value of V by 1, 2, … 10 mV from its steady state value
figure(2);
clear V;
Vthres=10;
for i=1:10
   y0=ySS+[0,0,0,i];
   [t,y]=ode45('HHode',[0,20],y0);

   % After each simulation, check whether the peak voltage surpassed 0mV, and if it did, plot the voltage with a red line,
    %but if it didn’t, plot it with a black line
   if max(y(:,4))>0
       plot(t,y(:,4),'r');
       hold on;
       Vthres=min(Vthres,i);    %record the voltage threshold
   else
       plot(t,y(:,4),'k');
       hold on;
   end
end

%Figure formatting
xlabel('t [ms]');ylabel('Transmembrane Voltage [mV]');title('Threshold Behavior');

% Plot all the membrane voltage trajectories on the same figure, like below
% We see that if the voltage threshold is surpassed, then the neuron ‘fires’ an action potential;
    % otherwise it just returns to the steady state value
% You can zoom in on your figure to see the threshold value (the voltage that separates the red lines from the black lines).

MIT6094-H3-5steadystate.png

MIT6094-H3-5threshold behavior.png

Problem 6 Optional – Recursion – Julia Set*

My second favorite problem. Just because it was made using a very simple process over and over again to produce intricate patterns. This particular problem uses a complex quadratic polynomial with complex parameter c, the simplest nontrivial case is:

f_{c}(z)=z^{2}+c

The Julia set of f(z) for c =-0.297491 + i0.641051 is shown below

MIT6094-H3-6juliasetpng.png

The escape velocity function

% a. It has been shown that if the modulus of zn becomes larger than 2 for some n then it is guaranteed
    % that the orbit will tend to infinity.
% The value of n for which this becomes true is called the ‘escape velocity’ of a particular z0 .
% Write a function that returns the escape velocity of a given z0 and c .
% The function declaration should be: n=escapeVelocity(z0,c,N) where N is the maximum allowed
% escape velocity (basically, if the modulus of zn does not exceed 2 for n<N return N as the escape velocity.
% This will prevent infinite loops). Use abs to calculate the modulus of a complex number
function n=escapeVelocity(z0,c,N)
    clear i;        %remove potential problems with variables vs imaginary number
    z=z0;
    n=0;            %invariant n : subcript of z
    while (abs(z)<2) && (n ~= N)
        zLast = z;      %store previous number, for speed
        z = z^2+c;      %the dynamic recursion
        n = n + 1;
        if z==zLast     %premature end if it becomes constant, for speed
            n = N;
        end
    end
end

The Julia Set

% Homework 3 Optional Problem 6 - Julia Sets
% Optional, but highly recommended: Julia Sets
% In this problem you will generate quadratic Julia Sets. The following description is adapted from Wikipedia at http://en.wikipedia.org/wiki/Julia
% For more information about Julia Sets please read the entire article there.
% Given two complex numbers, c and z0 , we define the following recursion:
% z_n = (z_(n?1))^2 + c
% This is a dynamical system known as a quadratic map
% Given a specific choice of c and z0 , the above recursion leads to a sequence of complex numbers z1, z2 , z3 … called the orbit of z0
% Depending on the exact choice of c and z0 , a large range of orbit patterns are possible
% For a given fixed c , most choices of z0 yield orbits that tend towards infinity
% (That is, the modulus z grows without limit as n increases.)
% For some values of c certain choices of z0 yield orbits n that eventually go into a periodic loop
% Finally, some starting values yield orbits that appear to dance around the complex plane, apparently at random
% (This is an example of chaos.) These starting values, z0 , make up the Julia set of the map, denoted Jc
% In this problem, you will write a MATLAB script that visualizes a slightly different set,
% called the filled-in Julia set (or Prisoner Set), denoted Kc , which is the set of all z0 with orbits which do not tend towards infinity
% The "normal" Julia set Jc is the edge of the filled-in Julia set
% The figure below illustrates a Julia Set for one particular value of c
% You will write MATLAB code that can generate such fractals in this problem.

% a. It has been shown that if the modulus of zn becomes larger than 2 for some n then it is guaranteed that the orbit will tend to infinity.
% The value of n for which this becomes true is called the ‘escape velocity’ of a particular z0 . Write a function that returns the escape velocity of a given z0 and c .
% The function declaration should be: n=escapeVelocity(z0,c,N) where N is the maximum allowed escape velocity (basically, if the modulus of zn does not exceed 2 for n<N return N as the escape velocity.
% This will prevent infinite loops). Use abs to calculate the modulus of a complex number

% b. To generate the filled in Julia Set, write the following function M=julia(zMax,c,N).
% zMax will be the maximum of the imaginary and complex parts of the various z0 ’s for which we will compute escape velocities.
% c and N are the same as defined above, and M is the matrix that contains the escape velocity of various z0 ‘s.
function M=hw3p6o(zMax,c,iteration,resolution)
    clear i;clc;
    % i. In this function, you first want to make a 500x500 matrix that contains complex numbers with real part between –zMax and zMax, and imaginary part between –zMax and zMax.
        % Call this matrix Z. Make the imaginary part vary along the y axis of this matrix. You can most easily do this by using linspace and meshgrid, but you can also do it with a loop.
    if nargin
    re=linspace(-zMax,zMax,resolution);
    im=1i*re;
    [Re,Im]=meshgrid(re,im);
    Z=Re+Im;
    % ii. For each element of Z, compute the escape velocity (by calling your escapeVelocity) and store it in the same location in a matrix M.
        % When done, the matrix M should be the same size as Z and contain escape velocities with values between 1 and N.
    M=zeros(size(Z)); %allocate M
    for x=1:size(Z,1)
        for y=1:size(Z,2)
            M(x,y)=escapeVelocity(Z(x,y),c,iteration);
        end
    end

    % iii. Run your julia function with various zMax, c, and N values to generate various fractals. To display the fractal nicely, use imagesc to visualize atan(0.1*M),
        % (taking the arctangent of M makes the image look nicer; you may also want to use axis xy so the y values aren’t flipped). WARNING: this function may take a while to run.
    close all;
    figure;
    imagesc(atan(0.1*M));
    xlabel('Re(z)');
    ylabel('Im(z)');
    axis xy;
    colormap(jet);
end

References

Chance, E. M., Clore, G. M., Curtis, A. R., & Shephard, E. P. (1981). Numerical solution of the Hodgkin-Huxley equations in a moving coordinate system. Simulation of nerve impulse transmission over long distances. Journal of Computational Physics, 40(2), 318-326.

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s