How to code and solve ODE equations many variables?

Hi there,

I am new to MATLAB so i'm doing a few tutorials but they don't have memos, so any help would be appreciated! heres the required and what I've done so far:

Write an ODE file to return the following derivatives

dndt= 1− n* αn(V) − n*βn(V)

dmdt=(1 − m)*αm(V) − m*βm(V)

dhdt =(1 − h)*αh(V) − h*βh(V)

dvdt=−(1/C)* (Gk*n^4(V− Ek)+ Gna*m^3*h*(V − Ena )+ Gl(V − El ))

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

and αn(V),βn(V),αm(V),βm(V),αh(V)and βh(V) are all functions whose code have been given.

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)

seems simple enough.... my code is as follows:

function rk=chem(t,V)

C=1 Gk=36 Gna=120 Gl=0.3 Ek=-72 Ena=55 El=-49.4 V=zeros(4,1)

V(1)= (1-n).*'alphan' -n.*'betan'

V(2)= (1-m).*'alpham'-m.*'betam' V(3)= (1-h).*'alphah'-h.*'betah'

V(4)= -(1/C)*(Gk*n^4*(V-Ek)+Gna*m^3*h*(V-Ena)+Gl*(V-El)) rk=[V(1);V(2);V(3);V(4)] end

timerange=[0 20]

initialvalues=[0.5 0.5 0.5 -60]

[t,y]=ode45(@chem,timerange,initialvalues)

plot(t,y(:,1))

but the damn thing keeps telling me that t is undefined! please help, I know it looks long but any help would be appreciated.

Answer 1

[Star Strider]

‘but the damn thing keeps telling me that t is undefined! please help, I know it looks long but any help would be appreciated.’

Please don’t be angry with it! The ODE solvers require the ODE function to return [t,y]. Change your function line from:

function rk=chem(t,V)

to:

function [t,V]=chem(t,V)

and if there are no other problems, your code should run. You might want to define alphan, betan, and n before you run it.

Also, put semicolons (;) at the end of your lines in chem. Otherwise, it will output the result of every line in every iteration to the Command Window.

[Star Strider]

First please format your code. Use the [{}Code] radio button.

Second, write your equations this way:

function [t V]=chem(t,V)

C=1; Gk=36; Gna=120; Gl=0.3; Ek=-72; Ena=55; El=-49.4;

V(1)=  (1-n).*'alphan' -n.*'betan' ;
V(2)=  (1-m).*'alpham'-m.*'betam' ; 
V(3)= (1-h).*'alphah'-h.*'betah';
V(4)= -(1/C)*(Gk*n^4*(V-Ek)+Gna*m^3*h*(V-Ena)+Gl*(V-El));  
end

Don’t write them as anonymous functions because they aren’t.

Once you define alphan, betan, n, m, and h, in chem, use this code (with those revisions) as your entire ODE function in ode45. It should work, providing there aren’t any other errors.

[Matthew]

 I made those changes and put all the alpha and beta functions in the chem file, and defined h,m and n. but it now says my matrix dimensions don't agree when im defining the derivatives. sorry about the messy code before. thanks for all your help so far. here it is:

     function [t V]=chem(t,V)

C=1;
Gk=36;
Gna=120;
Gl=0.3;
Ek=-72;
Ena=55;
El=-49.4;

n=0.5;
m=0.5;
h=0.5;
function a = alphah(V)

a = 0.07*exp(-0.05*(V+60));
end

function b = betah(V)

b = 1./(exp(-(V+30)/10)+1);

end

function b = betam(V)

b = 4*exp(-(V+60)/18);
end

function b = betan(V)

b = 0.125*exp(-(V+60)/80);

end

function a = alpham(V)

badIndices = find(abs(V+35)<1e-4);
goodIndices = find(abs(V+35)>=1e-4);

a(goodIndices) = -0.1*(V(goodIndices)+35)./(exp(-(V(goodIndices)+35)/10)-1);
if(max(size(badIndices))>0)
    a(badIndices) = 1./(1-(V(badIndices)+35)/20);
end
end

function a = alphan(V)

badIndices = find(abs(V+50)<1e-4);
goodIndices = find(abs(V+50)>=1e-4);

a(goodIndices) = -0.01*(V(goodIndices)+50)./(exp(-(V(goodIndices)+50)/10)-1);
if(max(size(badIndices))>0)
    a(badIndices) = .1./(1-(V(badIndices)+50)/20);
end

end

V(1)= (1-n).*'alphan' -n.*'betan' ;

V(2)= (1-m).*'alpham'-m.*'betam' ;
V(3)=(1-h).*'alphah'-h.*'betah';

V(4)= -(1/C)*(Gk*n^4*(V-Ek)+Gna*m^3*h*(V-Ena)+Gl*(V-El));

end

and this is the error message:

Error using - Matrix dimensions must agree.

Error in chem (line 60) V(1)= (1-n).*'alphan' -n.*'betan' ;

[Star Strider]

I am now totally lost! I have no idea what you are doing. I can tell by inspection though that MATLAB is going to have serious problems with your code.

Please see the documentation on ode45 and the other solvers for details on how to write your ODE function. Then see ‘The Hodgkin Huxley Model’. They also implement it in MATLAB.

Answer 2

[Matthew]

thanks for the help, I changed what you said but its now giving me an error in the first line where I define one of the derivatives, heres my altered code:

function [t V]=chem(t,V)

C=1; Gk=36; Gna=120; Gl=0.3; Ek=-72; Ena=55; El=-49.4;

V(1)= @(n,t) (1-n).*'alphan' -n.*'betan' ;

V(2)= @(m,t) (1-m).*'alpham'-m.*'betam' ; V(3)=@(h,t) (1-h).*'alphah'-h.*'betah';

V(4)=@(V,t) -(1/C)*(Gk*n^4*(V-Ek)+Gna*m^3*h*(V-Ena)+Gl*(V-El)); [t V]=[V(1);V(2);V(3);V(4)] end

timerange=[0 20]

initialvalues=[0.5 0.5 0.5 -60]

[t,V]=ode45(@chem,timerange,initialvalues)

plot(t,V(:,1))