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.
polyFittingLecture.m – this illustrates polynomial fitting, representation, and evaluation
%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
%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
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
% 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
% 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).
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:
The Julia set of f(z) for c =-0.297491 + i0.641051 is shown below
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
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.