how to use multivariate hidden markov models

I'm trying to identify states in a high-dimensional binary time series, I want to use Hidden Markov Modelling specifically. The binary time series represent the output of a neural network with 500 neurons. The entire network is always in one of five possible states. The documentation for hmmtrain and hmmgenerate is sparse and makes no mention of multivariate emissions. However hmmtrain does allow a matrix to me used as the sequence, suggesting that it can use a multivariate emissions sequence. It doesn't say what dimension to put the extra emissions variable on so I guessed the 3rd dimension. It does work (see below).

Here is a failed example of using hmmgenerate:

T=150;
N_neurons=500;
N_states=5;
mean_time_in_state=7;

transition_mat=(1.1-0.2*rand(N_states)).* (ones(N_states)+(mean_time_in_state-2).*eye(N_states))/mean_time_in_state; 
transition_mat=transition_mat./sum(transition_mat,2);

emission_mat=0.1*rand(N_states,1,N_neurons); % chance that spike is fired at each timestep
emission_mat=cat(2,1-emission_mat,emission_mat); % emmissions are all binary, but are multivariate

[seq,states] = hmmgenerate(T,transition_mat,emission_mat);
sample_spike_train=seq-1; % make 1 or 0 (per convention in neural networks)

"sample_spike" should have size 500x150. Instead it is 1x150. Nonetheless I can use "states" to generate a spike series and test "hmmtrain".

Here is a failed example of using hmmtrain:

sample_spike_train=squeeze(emission_mat(states,2,:))'>rand(N_neurons,T); % each row is a sequence from each of the neurons (time synchronized)

spike_rates(1,1,:)=1-mean(spike_raster2,2);
spike_rates(1,2,:)=mean(spike_raster2,2);
guess_emit_mat=repmat(spike_rates,N_states,1,1);

guess_transition_mat=ones(N_states)/N_states;

[ESTTR,ESTEMIT] = hmmtrain(1+sample_spike_train,guess_transition_mat,guess_emit_mat) ;

This yeilds obviously wrong results. ESTTR is identical to guess_transition_mat, and (bizarrely) ESTEMIT has size N_states x 394 irrespective of how many neurons I use. I expect ESTTR~transition_mat and ESTEMIT~emission_mat.

I did find a multivariate HMM toolbox written a long time ago by one "Kevin Murphy" https://www.cs.ubc.ca/~murphyk/Software/HMM/hmm.html

It variously produces poor results (the transistion matrix is not even close to correct) or runs into infinity/nan issues when T is large.

Here is a failed example of using an external package.

%create a hidden markov model
  N_neurons;      %Number of features
  T;      %Number of observations 
  nex = 1;    %Number of sequences 
  M = 1;      %Number of mixtures 
  N_states;      %Number of states 
  cov_type = 'full';

  % initial guess at multivariate parameters
prior0 = normalise(rand(N_states,1));
transmat0 = ones(N_states)/N_states;
  [mu0, Sigma0] = mixgauss_init(N_states*M, 1+sample_spike_train', cov_type);
  mixmat0 = mk_stochastic(rand(N_states,M));

[LL, prior1, transmat1, mu1, Sigma1, mixmat1] = ...
    mhmm_em(1+sample_spike_train', prior0, transmat0, mu0, Sigma0, mixmat0, 'max_iter', 100);% Fitting category1 to object 1

Any help, or explanation is appreciated. I've seen HMMs used in neuroscience for multiple neurons at once, but it looks like custom software was used and it is not available. Why does hmmtrain allow matrix input if it isn't multivariate?

https://www.pnas.org/content/104/47/18772.short#sec-6

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2576226/

http://papers.nips.cc/paper/2886-efficient-estimation-of-hidden-state-dynamics-from-spike-trains.pdf