Unable to perform assignment because the left and right sides have a different number of elements.

Hello! I am trying to model the calcium T-current using the hodgkin huxley model but am getting an error message "unable to perform assignment because the left and right sides have a different number of elements" for the line with the code: V(j+1) = V(j) +(k1v +k2v)*dt/2;

I am not sure how to fix this error so I was wondering if I could help on how to edit it. Thank you!

This is the equation that I am trying to model:

My code is:

% Biophysical parameters

C = 1;
El = -52;
Ena = 55;
Ek = -75;
ECa = 120;   % Calcium reversal potential (mV)
Gna = 120;
Gl = 0.3;
Gk = 36;
GCa = 0.3;
Iapp = 1.83;

D = 0;
phi=1;

minf=@(v) 1./(1+exp(-(v+40)/9));
hinf=@(v) 1./(1+exp((v+62)/10));
ninf=@(v) 1./(1+exp(-(v+53)/16));
cinf=@(V) 1./(1+exp(-(V+59)/6.2));
sinf=@(v) 1./(1+exp((v+83)/4));

taum=@(v) 0.3+0.000001*v;
tauh=@(v) 1+11./(1+exp((v+62)/10));
taun=@(v) 1+6./(1+exp((v+53)/16));
tauc=@(v) (22.7+0.27)./(exp((v + 48)/4)+exp(-(v + 407)/50));

% Voltage-dependent activation/inactivation curves & time constants

vv=-100:0.1:100;

figure
hold on
plot(vv,minf(vv),'b','linewidth',2);
plot(vv,hinf(vv),'r','linewidth',2);
plot(vv,ninf(vv),'g','linewidth',2);
plot(vv,cinf(vv),'k','linewidth',2);
plot(vv,sinf(vv),'y','linewidth',2);
axis([-100 120 -0.1 1.1]);
set(gca,'fontsize',20);
xlabel('V');
legend('m_{\infty}','h_{\infty}','n_{\infty}','c_{\infty}','s_{\infty}');

figure
hold on
plot(vv,taum(vv),'b','linewidth',2);
plot(vv,tauh(vv),'r','linewidth',2);
plot(vv,taun(vv),'g','linewidth',2);
plot(vv,tauc(vv),'y','linewidth',2);
axis([-100 120 -0.1 12]);
set(gca,'fontsize',20);
xlabel('V');
legend('\tau_{m}','\tau_{h}','\tau_{n}','\tau_{c}');

% Computation of the solution

Tmax = 3000;
dt = 0.05;
t = 0:dt:Tmax;

ti = dt;
tf = Tmax;
H = zeros(1,length(t));
H(floor(ti/dt):floor(tf/dt))=1;

% White noise

eta = randn(1,length(t));

% Numerical Computation

V = zeros(1,length(t));
m = zeros(1,length(t));
h = zeros(1,length(t));
n = zeros(1,length(t));
c = zeros(1,length(t));
s = zeros(1,length(t));

V(1) = 60;
m(1) = 0;
h(1) = 0;
n(1) = 0;
c(1) =1;
s(1)=1;

% V(1) = -61.8310;
% m(1) = 0.0816;
% h(1) = 0.4946;
% n(1) = 0.3661;

faraday=96520;
rgas=8.3134;
temp=293.15;
xi=V.*faraday*2/(rgas*1000*temp);
cao=2;
cai=1^-4;
I0=10;

cfedrive=0.002.*faraday.*xi.*(cai-cao*exp(-xi))/(1-exp(-xi));
i_cat=cinf(V).^2.*h*cfedrive;

for j=1:length(t)-1
    k1v = (Iapp*H(j)-Gl*(V(j)-El)-Gna*m(j)^3*h(j)*(V(j)-Ena)-Gk*n(j)^4*(V(j)-Ek)-I0-Gl*(V(j)-El)-i_cat)/C
    k1m = phi*(minf(V(j))-m(j))/taum(V(j));
    k1h = phi*(hinf(V(j))-h(j))/tauh(V(j));
    k1n = phi*(ninf(V(j))-n(j))/taun(V(j));
    k1c= phi*(sinf(V(j))-h(j))/tauc(V(j));
    av = V(j)+k1v*dt;
    av = av + sqrt(2*D*dt)*eta(j);
    am = m(j)+k1m*dt;
    ah = h(j)+k1h*dt;
    an = n(j)+k1n*dt;
    ac = s(j)+k1c*dt;
    k2v = (Iapp*H(j+1)-Gl*(av-El)-Gna*am^3*ah*(av-Ena)-Gk*an^4*(av-Ek)-I0-Gl*(ac-El)-I_cat)/C
    k2m = phi*(minf(av)-am)/taum(av);
    k2h = phi*(hinf(av)-ah)/tauh(av);
    k2n = phi*(ninf(av)-an)/taun(av);
    k2c = phi*(sinf(av)-an)/tauc(av);
    V(j+1) = V(j) + (k1v + k2v)*dt/2;
    V(j+1) = V(j+1) + sqrt(2*D*dt)*eta(j);
    m(j+1) = m(j) + (k1m + k2m)*dt/2;
    h(j+1) = h(j) + (k1h + k2h)*dt/2;
    n(j+1) = n(j) + (k1n + k2n)*dt/2;
    c(j+1) = s(j) + (k1c + k2c)*dt/2;
end

k1v = 1×60001

-75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700 -75.3700

Unrecognized function or variable 'I_cat'.

figure
plot(t,V,'b','linewidth',2);
axis([0 Tmax -80 80]);
set(gca,'fontsize',20);
xlabel('t');
ylabel('V');

aux = 0;
spkcnt = 0;
tspk = zeros(1);
if D>0
    for j=1:length(t)
        if V(j)>0 && aux==0
            spkcnt = spkcnt+1;
            tspk(spkcnt) = t(j);
            aux = 1;
        elseif V(j)<0 && aux==1
            aux = 0;
        end        
    end
    ISI = diff(tspk);
    figure
    histogram(ISI,20);
    set(gca,'fontsize',20);
    xlabel('ISI')
    title('Histogram');
end

if D > 0

    % "Manual Histogram"

    binsize = 10;
    Edges = 0:binsize:1000;
    [Hisi,Edgescnt] = histcounts(ISI,Edges);
    Hmax = max(Hisi);

    figure
    hold on
    histogram(ISI,Edges,'Facecolor','b');
    set(gca,'fontsize',20);
    xlabel('ISI')
    title('Histogram');
end

[Rik]

I recovered the removed content from the Bing cache (something which anyone can do). Editing away your question is very rude. Someone spent time reading your question, understanding your issue, figuring out the solution, and writing an answer. Now you repay that kindness by ensuring that the next person with a similar question can't benefit from this answer.

Answer 1

[Matt J]

Undoubtedly, it is because (k1v +k2v)*dt/2 is a vector, whereas your code pretends that it is a scalar. You should stop the code at that line and inspect k1v and k2v.

Answer 2

[Torsten]

k1v and k2v both are vectors of size 1x60001, but they must be scalars in the assignment.

Maybe you mean

V(j+1) = V(j) + (k1v(j) + k2v(j))*dt/2;

instead of

V(j+1) = V(j) + (k1v + k2v)*dt/2;

Or make i_cat a scalar (at the moment, it is a vector of size 1x60001).

[Michelle Park]

I tried adding (j) but the code still showed an error

[Torsten]

If you change I_cat to i_cat in

k2v = (Iapp*H(j+1)-Gl*(av-El)-Gna*am^3*ah*(av-Ena)-Gk*an^4*(av-Ek)-I0-Gl*(ac-El)-I_cat)/C

, it seems to work (at least for the MATLAB online runtime).

Or which error do you encounter ?

[Michelle Park]

When I changed I_cat to i_cat the same error showed up. When I tried adding (j) k1v and k2v, the code would not run. Sorry to bother you, how would I make i_cat a scalar?

[Torsten]

This code worked on MATLAB online. I changed Tmax from 3000 to 300.

% Biophysical parameters
C = 1;
El = -52;
Ena = 55;
Ek = -75;
ECa = 120;   % Calcium reversal potential (mV)
Gna = 120;
Gl = 0.3;
Gk = 36;
GCa = 0.3;
Iapp = 1.83;
D = 0;
phi=1;
minf=@(v) 1./(1+exp(-(v+40)/9));
hinf=@(v) 1./(1+exp((v+62)/10));
ninf=@(v) 1./(1+exp(-(v+53)/16));
cinf=@(V) 1./(1+exp(-(V+59)/6.2));
sinf=@(v) 1./(1+exp((v+83)/4));
taum=@(v) 0.3+0.000001*v;
tauh=@(v) 1+11./(1+exp((v+62)/10));
taun=@(v) 1+6./(1+exp((v+53)/16));
tauc=@(v) (22.7+0.27)./(exp((v + 48)/4)+exp(-(v + 407)/50));
% Voltage-dependent activation/inactivation curves & time constants
vv=-100:0.1:100;
figure
hold on
plot(vv,minf(vv),'b','linewidth',2);
plot(vv,hinf(vv),'r','linewidth',2);
plot(vv,ninf(vv),'g','linewidth',2);
plot(vv,cinf(vv),'k','linewidth',2);
plot(vv,sinf(vv),'y','linewidth',2);
axis([-100 120 -0.1 1.1]);
set(gca,'fontsize',20);
xlabel('V');
legend('m_{\infty}','h_{\infty}','n_{\infty}','c_{\infty}','s_{\infty}');

figure
hold on
plot(vv,taum(vv),'b','linewidth',2);
plot(vv,tauh(vv),'r','linewidth',2);
plot(vv,taun(vv),'g','linewidth',2);
plot(vv,tauc(vv),'y','linewidth',2);
axis([-100 120 -0.1 12]);
set(gca,'fontsize',20);
xlabel('V');
legend('\tau_{m}','\tau_{h}','\tau_{n}','\tau_{c}');

% Computation of the solution
Tmax = 300;
dt = 0.05;
t = 0:dt:Tmax;
ti = dt;
tf = Tmax;
H = zeros(1,length(t));
H(floor(ti/dt):floor(tf/dt))=1;
% White noise
eta = randn(1,length(t));
% Numerical Computation
V = zeros(1,length(t));
m = zeros(1,length(t));
h = zeros(1,length(t));
n = zeros(1,length(t));
c = zeros(1,length(t));
s = zeros(1,length(t));
V(1) = 60;
m(1) = 0;
h(1) = 0;
n(1) = 0;
c(1) =1;
s(1)=1;
% V(1) = -61.8310;
% m(1) = 0.0816;
% h(1) = 0.4946;
% n(1) = 0.3661;
faraday=96520;
rgas=8.3134;
temp=293.15;
xi=V.*faraday*2/(rgas*1000*temp);
cao=2;
cai=1^-4;
I0=10;
cfedrive=0.002.*faraday.*xi.*(cai-cao*exp(-xi))/(1-exp(-xi));
i_cat=cinf(V).^2.*h*cfedrive;
for j=1:length(t)-1
    k1v = (Iapp*H(j)-Gl*(V(j)-El)-Gna*m(j)^3*h(j)*(V(j)-Ena)-Gk*n(j)^4*(V(j)-Ek)-I0-Gl*(V(j)-El)-i_cat)/C;
    k1m = phi*(minf(V(j))-m(j))/taum(V(j));
    k1h = phi*(hinf(V(j))-h(j))/tauh(V(j));
    k1n = phi*(ninf(V(j))-n(j))/taun(V(j));
    k1c= phi*(sinf(V(j))-h(j))/tauc(V(j));
    av = V(j)+k1v*dt;
    av = av + sqrt(2*D*dt)*eta(j);
    am = m(j)+k1m*dt;
    ah = h(j)+k1h*dt;
    an = n(j)+k1n*dt;
    ac = s(j)+k1c*dt;
    k2v = (Iapp*H(j+1)-Gl*(av-El)-Gna*am^3*ah*(av-Ena)-Gk*an^4*(av-Ek)-I0-Gl*(ac-El)-i_cat)/C;
    k2m = phi*(minf(av)-am)/taum(av);
    k2h = phi*(hinf(av)-ah)/tauh(av);
    k2n = phi*(ninf(av)-an)/taun(av);
    k2c = phi*(sinf(av)-an)/tauc(av);
    V(j+1) = V(j) + (k1v(j) + k2v(j))*dt/2;
    V(j+1) = V(j+1) + sqrt(2*D*dt)*eta(j);
    m(j+1) = m(j) + (k1m + k2m)*dt/2;
    h(j+1) = h(j) + (k1h + k2h)*dt/2;
    n(j+1) = n(j) + (k1n + k2n)*dt/2;
    c(j+1) = s(j) + (k1c + k2c)*dt/2;
end

figure
plot(t,V,'b','linewidth',2);
axis([0 Tmax -80 80]);
set(gca,'fontsize',20);
xlabel('t');
ylabel('V');

aux = 0;
spkcnt = 0;
tspk = zeros(1);
if D>0
    for j=1:length(t)
        if V(j)>0 && aux==0
            spkcnt = spkcnt+1;
            tspk(spkcnt) = t(j);
            aux = 1;
        elseif V(j)<0 && aux==1
            aux = 0;
        end        
    end
    ISI = diff(tspk);
    figure
    histogram(ISI,20);
    set(gca,'fontsize',20);
    xlabel('ISI')
    title('Histogram');
end
if D > 0
    % "Manual Histogram"
    binsize = 10;
    Edges = 0:binsize:1000;
    [Hisi,Edgescnt] = histcounts(ISI,Edges);
    Hmax = max(Hisi);
    figure
    hold on
    histogram(ISI,Edges,'Facecolor','b');
    set(gca,'fontsize',20);
    xlabel('ISI')
    title('Histogram');
end

[Torsten]

I just saw that you tried to define i_cat as a function of V. This doesn't work as you've done it in your code.

Try

% Biophysical parameters
C = 1;
El = -52;
Ena = 55;
Ek = -75;
ECa = 120;   % Calcium reversal potential (mV)
Gna = 120;
Gl = 0.3;
Gk = 36;
GCa = 0.3;
Iapp = 1.83;
D = 0;
phi=1;
minf=@(v) 1./(1+exp(-(v+40)/9));
hinf=@(v) 1./(1+exp((v+62)/10));
ninf=@(v) 1./(1+exp(-(v+53)/16));
cinf=@(V) 1./(1+exp(-(V+59)/6.2));
sinf=@(v) 1./(1+exp((v+83)/4));
taum=@(v) 0.3+0.000001*v;
tauh=@(v) 1+11./(1+exp((v+62)/10));
taun=@(v) 1+6./(1+exp((v+53)/16));
tauc=@(v) (22.7+0.27)./(exp((v + 48)/4)+exp(-(v + 407)/50));
% Voltage-dependent activation/inactivation curves & time constants
vv=-100:0.1:100;
figure
hold on
plot(vv,minf(vv),'b','linewidth',2);
plot(vv,hinf(vv),'r','linewidth',2);
plot(vv,ninf(vv),'g','linewidth',2);
plot(vv,cinf(vv),'k','linewidth',2);
plot(vv,sinf(vv),'y','linewidth',2);
axis([-100 120 -0.1 1.1]);
set(gca,'fontsize',20);
xlabel('V');
legend('m_{\infty}','h_{\infty}','n_{\infty}','c_{\infty}','s_{\infty}');

figure
hold on
plot(vv,taum(vv),'b','linewidth',2);
plot(vv,tauh(vv),'r','linewidth',2);
plot(vv,taun(vv),'g','linewidth',2);
plot(vv,tauc(vv),'y','linewidth',2);
axis([-100 120 -0.1 12]);
set(gca,'fontsize',20);
xlabel('V');
legend('\tau_{m}','\tau_{h}','\tau_{n}','\tau_{c}');

% Computation of the solution
Tmax = 300;
dt = 0.05;
t = 0:dt:Tmax;
ti = dt;
tf = Tmax;
H = zeros(1,length(t));
H(floor(ti/dt):floor(tf/dt))=1;
% White noise
eta = randn(1,length(t));
% Numerical Computation
V = zeros(1,length(t));
m = zeros(1,length(t));
h = zeros(1,length(t));
n = zeros(1,length(t));
c = zeros(1,length(t));
s = zeros(1,length(t));
V(1) = 60;
m(1) = 0;
h(1) = 0;
n(1) = 0;
c(1) =1;
s(1)=1;
% V(1) = -61.8310;
% m(1) = 0.0816;
% h(1) = 0.4946;
% n(1) = 0.3661;
faraday=96520;
rgas=8.3134;
temp=293.15;
xi=V.*faraday*2/(rgas*1000*temp);
cao=2;
cai=1^-4;
I0=10;
cfedrive=0.002.*faraday.*xi.*(cai-cao*exp(-xi))/(1-exp(-xi));
i_cat=@(v,h)cinf(v).^2.*h*cfedrive;
for j=1:length(t)-1
    k1v = (Iapp*H(j)-Gl*(V(j)-El)-Gna*m(j)^3*h(j)*(V(j)-Ena)-Gk*n(j)^4*(V(j)-Ek)-I0-Gl*(V(j)-El)-i_cat(V(j),h(j)))/C;
    k1m = phi*(minf(V(j))-m(j))/taum(V(j));
    k1h = phi*(hinf(V(j))-h(j))/tauh(V(j));
    k1n = phi*(ninf(V(j))-n(j))/taun(V(j));
    k1c= phi*(sinf(V(j))-h(j))/tauc(V(j));
    av = V(j)+k1v*dt;
    av = av + sqrt(2*D*dt)*eta(j);
    am = m(j)+k1m*dt;
    ah = h(j)+k1h*dt;
    an = n(j)+k1n*dt;
    ac = s(j)+k1c*dt;
    k2v = (Iapp*H(j+1)-Gl*(av-El)-Gna*am^3*ah*(av-Ena)-Gk*an^4*(av-Ek)-I0-Gl*(ac-El)-i_cat(V(j),h(j)))/C;
    k2m = phi*(minf(av)-am)/taum(av);
    k2h = phi*(hinf(av)-ah)/tauh(av);
    k2n = phi*(ninf(av)-an)/taun(av);
    k2c = phi*(sinf(av)-an)/tauc(av);
    V(j+1) = V(j) + (k1v + k2v)*dt/2;
    V(j+1) = V(j+1) + sqrt(2*D*dt)*eta(j);
    m(j+1) = m(j) + (k1m + k2m)*dt/2;
    h(j+1) = h(j) + (k1h + k2h)*dt/2;
    n(j+1) = n(j) + (k1n + k2n)*dt/2;
    c(j+1) = s(j) + (k1c + k2c)*dt/2;
end
figure
plot(t,V,'b','linewidth',2);
axis([0 Tmax -80 80]);
set(gca,'fontsize',20);
xlabel('t');
ylabel('V');

aux = 0;
spkcnt = 0;
tspk = zeros(1);
if D>0
    for j=1:length(t)
        if V(j)>0 && aux==0
            spkcnt = spkcnt+1;
            tspk(spkcnt) = t(j);
            aux = 1;
        elseif V(j)<0 && aux==1
            aux = 0;
        end        
    end
    ISI = diff(tspk);
    figure
    histogram(ISI,20);
    set(gca,'fontsize',20);
    xlabel('ISI')
    title('Histogram');
end
if D > 0
    % "Manual Histogram"
    binsize = 10;
    Edges = 0:binsize:1000;
    [Hisi,Edgescnt] = histcounts(ISI,Edges);
    Hmax = max(Hisi);
    figure
    hold on
    histogram(ISI,Edges,'Facecolor','b');
    set(gca,'fontsize',20);
    xlabel('ISI')
    title('Histogram');
end