Я очень близок к преобразованию этого сценария Matlab в R и успешно преобразовал несколько похожих сценариев. Тем не менее, этот я просто не могу найти ошибку. Кто-нибудь может это заметить?
Один из способов проверить - запустить оба сценария и посмотреть, совпадает ли сюжет. Сюжет Matlab правильный.
Обратите внимание, что я использую плагин Matlab для R, так что часть кода может оставаться в такте.
Matlab:
% Excitatory neurons Inhibitory neurons
Ne=800; Ni=200;
re=rand(Ne,1); ri=rand(Ni,1);
a=[0.02*ones(Ne,1); 0.02+0.08*ri];
b=[0.2*ones(Ne,1); 0.25-0.05*ri];
c=[-65+15*re.^2; -65*ones(Ni,1)];
d=[8-6*re.^2; 2*ones(Ni,1)];
S=[0.5*rand(Ne+Ni,Ne),-rand(Ne+Ni,Ni)];
v=-65*ones(Ne+Ni,1); % Initial values of v
u=b.*v; % Initial values of u
firings=[]; % spike timings
for t=1:1000 % simulation of 1000 ms
I=[5*randn(Ne,1);2*randn(Ni,1)]; % thalamic input
fired=find(v>=30); % indices of spikes
ding = fired;
if ~isempty(fired)
firings=[firings; t+0*fired, fired];
v(fired)=c(fired);
u(fired)=u(fired)+d(fired);
jaja = S(:,fired);
I=I+sum(S(:,fired),2);
end;
v=v+0.5*(0.04*v.^2+5*v+140-u+I);
v=v+0.5*(0.04*v.^2+5*v+140-u+I);
u=u+a.*(b.*v-u);
end;
plot(firings(:,1),firings(:,2),'.');
xlabel('Time (ms)'); ylabel('Neurons');
R :
require(matlab)
require(ramify)
require(ggplot2)
Ne <- 800
Ni <- 200
re <- rand(Ne,1)
ri <- rand(Ni,1)
# Excitatory neurons # Inhibitory neurons (i.e different column vectors)
a <- matrix(c(0.02*ones(Ne,1), 0.02+0.08*ri))
b <- matrix(c(0.2*ones(Ne,1) , 0.25-0.05*ri))
c <- matrix(c(-65+15*re^2, 65*ones(Ni,1)))
d <- matrix(c(8-6*re^2, 2*ones(Ni,1)))
S <- cbind(0.5*rand(Ne+Ni,Ne),-rand(Ne+Ni,Ni))
v <- -65*ones(Ne+Ni,1) # Initial values of v
u <- b*v # Initial values of u
firings <- matrix(ncol=2) # spike timings
for (t in 1:1000) { # simulation of 1000 ms
I <- matrix(c(5*randn(Ne,1), 2*randn(Ni,1))) # thalamic input
fired <- matlab::find(v>=30) # indices of spikes
if (!isempty(fired)) {
firings <- rbind(firings,cbind(t+0*fired, fired))
v[fired] <- c[fired]
u[fired] <- u[fired]+d[fired]
I <- I + sum(S[,fired])
}
v <- v + 0.5 * (0.04*v^2+5*v+140-u+I)
v <- v + 0.5 * (0.04*v^2+5*v+140-u+I)
u <- u + a * (b*v-u)
}
plot(firings[,1], firings[,2])
График должен выглядеть следующим образом:
Одна из проблем заключается в том, что переменная срабатывания приводит к 412369 x 2 матрицы в R, но матрица 7383x2 в Matlab. Я просто не могу найти, где это происходит.