Matlab to R conversion

Question

I am very close to converting this Matlab script to R, and successfully converted several similar scripts. However, this one I just cannot find the bug. Can anyone spot it?

A way to test is to run both scripts and see if the plot is the same. The Matlab plot is correct.

Note that I use a Matlab plugin for R so that some of the code can stay in tact.

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])

The graph should look like this:

One issue is that the firings variable results in a 412369 x 2 matrix in R, but a 7383x2 matrix in Matlab. I simply cannot find where this takes place.

Answer 1

I went through you codes both in R and MATLAB, finding two places that you did not correctly translated, which should be

• c <- matrix(c(-65+15*re^2, 65*ones(Ni,1))) to c <- matrix(c(-65+15*re^2, -65*ones(Ni,1)))
• sum(S[,fired]) to rowSums(S[,fired,drop = FALSE])

Below is the plot I got from R