Thank you @lmazzuca for all your suggestions.
We are already fitting a different HMM to each single area (the results we attached follow this too). However, with all the areas (including secondary motor cortex, striatum, and others in the midbrain) we have fitted until now we have obtained results similar to the ones shown in the previous figure.
For the fitting we are using your matlab code, and we looked directly into the generated plots with the posterior probabilities for each trial. Then, we used the states (hmm_postfit.sequences) with higher confidence to make our plot (in which zero implies no state with confidence at all). For example, in the figure we showed, we have a lot of uncertainty in trials 40-100.
Considering that we were already following your suggestions (perhaps we are still misunderstanding parameters or any other concept), we were worried about the fact of the dynamics related to the task are not captured by the model. Also, due to the computational cost, exploration of different datasets becomes difficult.
We attach one example of the last simulation we run with the MRN area (only including correct trials). At least it seems we are capturing some pre and post-stimulus dynamics, but we still have a lot of uncertainty and variability. If you can see a clear problem that we are missing in the transition matrix (sparseness or whatever) or in any other thing, please let us know. Also, note that we have changed the color coding of the states to show uncertain states in white.
Thank you very much for your help and sorry for replying during the weekend,