# How to use interpolation instead of integrating in chunks with ODE45?

1 view (last 30 days)

Show older comments

Samuele Bolotta
on 21 Mar 2021

Answered: Bjorn Gustavsson
on 22 Mar 2021

I have written this code to integrate voltage V and gating conductances m, n and h in 30ms with injection of current from ms 10 to 20. The code works well and it shows the behaviour I want to see. I know it's suboptimal, but just out of curiosity I'm trying to implement the same code with interpolation instead of integration in chunks. However, I'm failing to do so.

This is the integration with interpolation, which is giving the incorrect result:

myode = @(T,y0) [((1/Cm)*interp1(t, I, T, 'linear','extrap')-(INa+IK+Il)); % Normal case

alpham(y0(1,1))*(1-y0(2,1))-betam(y0(1,1))*y0(2,1);

alphan(y0(1,1))*(1-y0(3,1))-betan(y0(1,1))*y0(3,1);

alphah(y0(1,1))*(1-y0(4,1))-betah(y0(1,1))*y0(4,1)];

[time,V] = ode45(myode, t, y0, I);

This is the correct one, in chunks:

[time,V] = ode45(@ODEMAT, [t(chunk), t(chunk+1)], y);

%%% Calls ODEMAT, in which theres' the integration (extensive code posted below):

dydt = [((1/Cm)*(I(chunk)-(INa+IK+Il))); % Normal case

alpham(V)*(1-m)-betam(V)*m;

alphan(V)*(1-n)-betan(V)*n;

alphah(V)*(1-h)-betah(V)*h];

This is the code with the integration in chunks:

function Z2_chunks ()

%% Initial values

V=-60; % Initial Membrane voltage

m1=alpham(V)/(alpham(V)+betam(V)); % Initial m-value

n1=alphan(V)/(alphan(V)+betan(V)); % Initial n-value

h1=alphah(V)/(alphah(V)+betah(V)); % Initial h-value

y0=[V;m1;n1;h1];

t = [1:30];

I = [zeros(1,10),ones(1,10),zeros(1,10)];

% Plotting purposes (set I(idx) equal to last value of I)

idx = numel(t);

I(idx) = 0.1;

chunks = numel(t) - 1;

for chunk = 1:chunks

if chunk == 1

V=-60; % Initial Membrane voltage

m=alpham(V)/(alpham(V)+betam(V)); % Initial m-value

n=alphan(V)/(alphan(V)+betan(V)); % Initial n-value

h=alphah(V)/(alphah(V)+betah(V)); % Initial h-value

y=[V;m;n;h];

else

y = V(end, :); % Final position is initial value for next interval

end

[time,V] = ode45(@ODEMAT, [t(chunk), t(chunk+1)], y);

if chunk == 1

def_time = time;

def_v = V;

else

def_time = [def_time; time];

def_v = [def_v; V];

end

end

OD = def_v(:,1);

ODm = def_v(:,2);

ODn = def_v(:,3);

ODh = def_v(:,4);

time = def_time;

%% Plots

%% Voltage

figure

subplot(3,1,1)

plot(time,OD);

legend('ODE45 solver');

xlabel('Time (ms)');

ylabel('Voltage (mV)');

title('Voltage Change for Hodgkin-Huxley Model');

%% Current

subplot(3,1,2)

stairs(t,I)

ylim([0 5*max(I)])

legend('Current injected')

xlabel('Time (ms)')

ylabel('Ampere')

title('Current')

%% Gating variables

subplot(3,1,3)

plot(time,[ODm,ODn,ODh]);

legend('ODm','ODn','ODh');

xlabel('Time (ms)')

ylabel('Value')

title('Gating variables')

function [dydt] = ODEMAT(t,y)

%% Constants

ENa=55; % mv Na reversal potential

EK=-72; % mv K reversal potential

El=-49; % mv Leakage reversal potential

%% Values of conductances

gbarl=0.003; % mS/cm^2 Leakage conductance

gbarNa=1.2; % mS/cm^2 Na conductance

gbarK=0.36; % mS/cm^2 K conductancence

Cm = 0.01; % Capacitance

% Values set to equal input values

V = y(1);

m = y(2);

n = y(3);

h = y(4);

gNa = gbarNa*m^3*h;

gK = gbarK*n^4;

gL = gbarl;

INa=gNa*(V-ENa);

IK=gK*(V-EK);

Il=gL*(V-El);

dydt = [((1/Cm)*(I(chunk)-(INa+IK+Il))); % Normal case

alpham(V)*(1-m)-betam(V)*m;

alphan(V)*(1-n)-betan(V)*n;

alphah(V)*(1-h)-betah(V)*h];

end

end

This is the code with the integration with interpolation. As you can see, the plots are really different:

function Z1_interpol ()

%% Initial values

V=-60; % Initial Membrane voltage

m1=alpham(V)/(alpham(V)+betam(V)); % Initial m-value

n1=alphan(V)/(alphan(V)+betan(V)); % Initial n-value

h1=alphah(V)/(alphah(V)+betah(V)); % Initial h-value

y0=[V;m1;n1;h1];

t = [1:30];

I = [zeros(1,10),ones(1,10),zeros(1,10)];

% Plotting purposes (set I(idx) equal to last value of I)

idx = numel(t);

I(idx) = 0.1;

V=-60; % Initial Membrane voltage

m=alpham(V)/(alpham(V)+betam(V)); % Initial m-value

n=alphan(V)/(alphan(V)+betan(V)); % Initial n-value

h=alphah(V)/(alphah(V)+betah(V)); % Initial h-value

y=[V;m;n;h];

%% Constants

ENa=55; % mv Na reversal potential

EK=-72; % mv K reversal potential

El=-49; % mv Leakage reversal potential

%% Values of conductances

gbarl=0.003; % mS/cm^2 Leakage conductance

gbarNa=1.2; % mS/cm^2 Na conductance

gbarK=0.36; % mS/cm^2 K conductancence

Cm = 0.01; % Capacitance

%% Initial values

V=-60; % Initial Membrane voltage

m=alpham(V)/(alpham(V)+betam(V)); % Initial m-value

n=alphan(V)/(alphan(V)+betan(V)); % Initial n-value

h=alphah(V)/(alphah(V)+betah(V)); % Initial h-value

y0=[V;m;n;h];

gNa = gbarNa*m^3*h;

gK = gbarK*n^4;

gL = gbarl;

INa=gNa*(V-ENa);

IK=gK*(V-EK);

Il=gL*(V-El);

myode = @(T,y0) [((1/Cm)*interp1(t, I, T, 'linear','extrap')-(INa+IK+Il)); % Normal case

alpham(y0(1,1))*(1-y0(2,1))-betam(y0(1,1))*y0(2,1);

alphan(y0(1,1))*(1-y0(3,1))-betan(y0(1,1))*y0(3,1);

alphah(y0(1,1))*(1-y0(4,1))-betah(y0(1,1))*y0(4,1)];

[time,V] = ode45(myode, t, y0, I);

OD=V(:,1);

ODm=V(:,2);

ODn=V(:,3);

ODh=V(:,4);

time = time;

%% Plots

%% Voltage

figure

subplot(3,1,1)

plot(time,OD);

legend('ODE45 solver');

xlabel('Time (ms)');

ylabel('Voltage (mV)');

title('Voltage Change for Hodgkin-Huxley Model');

%% Current

subplot(3,1,2)

stairs(t,I)

ylim([0 5*max(I)])

legend('Current injected')

xlabel('Time (ms)')

ylabel('Ampere')

title('Current')

%% Gating variables

subplot(3,1,3)

plot(time,[ODm,ODn,ODh]);

legend('ODm','ODn','ODh');

xlabel('Time (ms)')

ylabel('Value')

title('Gating variables')

end

Thanks!

##### 4 Comments

darova
on 21 Mar 2021

### Accepted Answer

Bjorn Gustavsson
on 22 Mar 2021

Too long code. But this is what I use for cases where I have to integrate ODEs with time-varying data for one parameter:

function drdt = myODE(t,r,data,t4data)

g = -9.82;

D = interp1(t4data,data,t,'pchip');

dydt = [r(3);

-D*r(3)*abs(r(3));

r(4);

g-D*r(4)*abs(r(4))];

end

Perhaps this is the type of solution you're looking for.

HTH

##### 0 Comments

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!