Gaussian-Linear Hidden Markov Model

This notebook shows an example of training and inspecting a Gaussian-Linear Hidden Markov Model (GLHMM). This model is fit to two sets of timeseries, such as a neuroimaging/electrophysiological recording and a corresponding behavioural or other physiological timeseries. If you want to model just one set of timeseries, you can e.g., use the Standard Gaussian HMM. If you are new to the HMM or the GLHMM toolbox, we recommend starting there for a thorough introduction.

Outline

1. Background

2. Example: Modelling time-varying interaction between brain and physiological data

   • Preparation

   • Load data

   • Initialise and train a GLHMM

   • Inspect model

     • Interaction between brain and physiological measures

     • State means and covariances

     • Dynamics: Transition probabilities and Viterbi path

Background

The GLHMM is a generalization of the HMM, introduced in Vidaurre et al., 2023. We here assume that the observations Y at time point t were generated by a Gaussian distribution with parameters \(\mu\) and \(\Sigma\) (similar to the standard Gaussian HMM), as well as regression coefficients \(\beta\) that relate the second set of variables X to Y. When state k is active at time point t, the GLHMM thus assumes that the timeseries Y follows the following distribution:

\[Y_t\sim N(\mu^k+X_t\beta^k,\Sigma^k)\]

Compared to the standard Gaussian HMM, the GLHMM thus adds the \(X_t\beta^k\) term, which allows modelling the relationship to a second set of variables.

The remaining HMM parameters are essentially the same as in the standard HMM, i.e. the transition probabilities \(\theta\):

\[P(s_t=k|s_{t-1}=l)=\theta_{k,l}\]

the initial state probabilities \(\pi\):

\[P(s_t=k)=\pi_k\]

as well as the posterior estimates for both X and Y:

\[\gamma_{t,k}:=P(s_t=k|s_{>t},s_{<t},X,Y)\]

\[\xi_{t,k,l}:=P(s_t=k,s_{t-1}=l|s_{>t},s_{<t-1},X,Y)\]

The GLHMM can be used to model, in addition to the patterns described by the standard HMM (such as time-varying amplitude or functional connectivity), temporal changes in the relationship between two timeseries. This could be, for instance, the interaction between one group of brain areas in the prefrontal cortex and another group of brain areas in the occipital cortex, the relationship between BOLD-signal across the whole brain and respiration, or the interaction between EEG recordings from two participants recorded simultaneously.

Example: Modelling time-varying interaction between brain and physiological data

We will now go through an example illustrating how to fit and inspect a GLHMM. The example uses simulated data that can be found in the example_data folder. The data were generated to resemble one set of fMRI timeseries and two corresponding non-brain physiological (e.g., heart rate and respiration) timeseries. Our goal is to estimate time-varying amplitude and functional connectivity (FC) within the fMRI recordings and temporal changes in the relationship between the fMRI and the physiological data.

Preparation

If you dont have the GLHMM-package installed, then run the following command in your terminal:

pip install glhmm

We then need to import the relevant modules:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sb
from glhmm import glhmm, preproc, utils, graphics

Load data

The GLHMM requires two inputs: a timeseries Yand a second timeseries X. When running the model on a concatenated timeseries, e.g., a group of subjects or several scanning sessions, you also need to provide the indices indicating where each subject/session in the concatenated timeseries Y starts and ends. Loading data and data formats are explained in more detail in the Standard Gaussian HMM tutorial.

Synthetic data for this example are provided in the glhmm/docs/notebooks/example_data folder. The file data.csv contains synthetic data for the first set of timeseries Y. This is the brain data, containing fMRI recordings from 20 subjects with 1,000 timepoints each, concatenated along the first dimension, and 50 brain areas. The file dataX.csv contains the corresponding physiological measures, heart rate and respiration. The sessions have the same duration as the brain data, i.e., there are 20 concatenated subjects with 1,000 timepoints each, but only 2 variables (heart rate and respiration). The file T.csv specifies the beginning and end of each subject’s session (same for X and Y).

brain_data = pd.read_csv('./example_data/data.csv', header=None).to_numpy()
phys_data = pd.read_csv('./example_data/dataX.csv', header=None).to_numpy()
T_t = pd.read_csv('./example_data/T.csv', header=None).to_numpy()

NOTE: It is important to standardise your timeseries and, if necessary, apply other kinds of preprocessing before fitting the model.

This will be done separately for each session/subject as specified in the indices. The data provided here are already close to standardised (so the code below will not do much), but see Prediction tutorial to see the effect on real data.

brain_data,_ = preproc.preprocess_data(brain_data, T_t)
phys_data,_ = preproc.preprocess_data(phys_data, T_t)

Initialise and train a GLHMM

We first initialise the glhmm object, which we here call brainphys_glhmm. By specifying the parameters of the glhmm object, we define which type of model we want to fit and how states should be defined. In the case of the GLHMM, we need to set the model_beta parameter to indicate that we wish to model an interaction between two sets of variables. There are two options for modelling the interaction: global, meaning we estimate regression coefficients that are static over the timeseries in which case model_beta='shared', or state-dependent, meaning we estimate regression coefficients that vary over time in which case model_beta='state'. In this example, we want to model the time-varying interaction between the brain data and the physiological data, so we set model_beta='state'. For the other parameters, we here want to use the same set-up as in the standard Gaussian HMM, i.e., each state will be defined in terms of mean and covariance. That means, we will also estimate the time-varying amplitude and functional connectivity in the brain timeseries. We model 4 states by setting the parameter K=4, which you can compare to the states from the standard Gaussian HMM.

brainphys_glhmm = glhmm.glhmm(model_beta='state', K=4, covtype='full')

We can check the hyperparameters of the object to make sure the model is defined as we planned:

print(brainphys_glhmm.hyperparameters)

{'K': 4, 'covtype': 'full', 'model_mean': 'state', 'model_beta': 'state', 'dirichlet_diag': 10, 'connectivity': None, 'Pstructure': array([[ True,  True,  True,  True],
       [ True,  True,  True,  True],
       [ True,  True,  True,  True],
       [ True,  True,  True,  True]]), 'Pistructure': array([ True,  True,  True,  True])}

We then train the model. To model the interaction, we need to train the GLHMM using both the brain data and the physiological data. The “main” timeseries we are looking to model, i.e., the data for which we also want to estimate mean and covariance, is Y, so in this case Y=brain_data. The secondary timeseries from which we model only the interaction with the main timeseries is called X so X=phys_data. When fitting the model to a group of subjects or sessions, we also need to specify the indices of the start and end of each session (here called T_t). Note that the indices are shared for X and Y, so the beginning and end of the sessions in X and Y need to correspond.

brainphys_glhmm.train(X=phys_data, Y=brain_data, indices=T_t)

Init repetition 1 free energy = 1382582.111953315
Init repetition 2 free energy = 1381255.915434804
Init repetition 3 free energy = 1381968.4204982717
Init repetition 4 free energy = 1381954.0127431145
Init repetition 5 free energy = 1382507.1131980629
Best repetition: 2
Cycle 1 free energy = 1381159.0951647148
Cycle 2 free energy = 1380383.240749712
Cycle 3, free energy = 1380006.3504278269, relative change = 0.3269503731612878
Cycle 4, free energy = 1379686.2278781703, relative change = 0.21734649997392563
Cycle 5, free energy = 1379381.6827700434, relative change = 0.1713418388664007
Cycle 6, free energy = 1379127.4059850879, relative change = 0.12515535718027884
Cycle 7, free energy = 1378913.8422795346, relative change = 0.09511788492197021
Cycle 8, free energy = 1378748.607515083, relative change = 0.06854827257748364
Cycle 9, free energy = 1378604.8125436073, relative change = 0.05629563866087043
Cycle 10, free energy = 1378461.980603623, relative change = 0.052957312991007324
Cycle 11, free energy = 1378334.721119588, relative change = 0.045057588691036106
Cycle 12, free energy = 1378227.4081101185, relative change = 0.03660452410883656
Cycle 13, free energy = 1378155.895424771, relative change = 0.023812164204874295
Cycle 14, free energy = 1378094.9364589371, relative change = 0.019894193378045797
Cycle 15, free energy = 1378025.3937514941, relative change = 0.022191874168231786
Cycle 16, free energy = 1377964.048344247, relative change = 0.019200159088152507
Cycle 17, free energy = 1377909.2980234837, relative change = 0.01684730411898447
Cycle 18, free energy = 1377866.3687983947, relative change = 0.01303759265517248
Cycle 19, free energy = 1377828.808854305, relative change = 0.011278292791922752
Cycle 20, free energy = 1377798.1014032778, relative change = 0.009136420120574806
Cycle 21, free energy = 1377768.2823060234, relative change = 0.008794085222955324
Cycle 22, free energy = 1377730.613546661, relative change = 0.010987009282461886
Cycle 23, free energy = 1377702.6380304338, relative change = 0.008093696852145612
Cycle 24, free energy = 1377679.0804313878, relative change = 0.006769396353535655
Cycle 25, free energy = 1377657.042579992, relative change = 0.006292838517594068
Cycle 26, free energy = 1377634.8069956757, relative change = 0.006309241256617101
Cycle 27, free energy = 1377615.7110052563, relative change = 0.005389195627692172
Cycle 28, free energy = 1377600.9551944989, relative change = 0.0041470574178971474
Cycle 29, free energy = 1377587.9228006254, relative change = 0.003649332080547571
Cycle 30, free energy = 1377573.0552210677, relative change = 0.004145960388433189
Cycle 31, free energy = 1377559.6928766614, relative change = 0.0037123787053929954
Cycle 32, free energy = 1377549.0240116878, relative change = 0.0029553060095899183
Cycle 33, free energy = 1377538.9476217818, relative change = 0.0027834196773739352
Cycle 34, free energy = 1377527.1775483335, relative change = 0.003240732497681216
Cycle 35, free energy = 1377516.807483455, relative change = 0.002847129547717743
Cycle 36, free energy = 1377508.8880270675, relative change = 0.0021695909543666915
Cycle 37, free energy = 1377498.3199301306, relative change = 0.002886846708608011
Cycle 38, free energy = 1377485.7716303766, relative change = 0.003416061677316688
Cycle 39, free energy = 1377473.251468426, relative change = 0.003396823897633399
Cycle 40, free energy = 1377458.7211633923, relative change = 0.003926712550836546
Cycle 41, free energy = 1377445.360765979, relative change = 0.003597564063256813
Cycle 42, free energy = 1377436.209226455, relative change = 0.002458184235461908
Cycle 43, free energy = 1377428.6383459598, relative change = 0.0020294781210659614
Cycle 44, free energy = 1377421.987087933, relative change = 0.001779787442605406
Cycle 45, free energy = 1377416.1551511774, relative change = 0.0015581165433029653
Cycle 46, free energy = 1377410.2474398022, relative change = 0.001575873924142783
Cycle 47, free energy = 1377402.8537476328, relative change = 0.0019683751251223175
Cycle 48, free energy = 1377394.4627820465, relative change = 0.002228893749340791
Cycle 49, free energy = 1377389.711051598, relative change = 0.0012606118946676386
Cycle 50, free energy = 1377385.0707803369, relative change = 0.0012295286909934388
Cycle 51, free energy = 1377379.5402138445, relative change = 0.0014632851127292688
Cycle 52, free energy = 1377373.7801129224, relative change = 0.0015216965677279878
Cycle 53, free energy = 1377368.1245746121, relative change = 0.0014918444171163018
Cycle 54, free energy = 1377363.2829751906, relative change = 0.001275510794484917
Cycle 55, free energy = 1377357.781955106, relative change = 0.00144713676071088
Cycle 56, free energy = 1377351.2342237076, relative change = 0.0017195300720829244
Cycle 57, free energy = 1377344.5150379678, relative change = 0.0017614483158033911
Cycle 58, free energy = 1377338.2740037923, relative change = 0.0016334274525414475
Cycle 59, free energy = 1377333.3748893393, relative change = 0.001280573094832465
Cycle 60, free energy = 1377329.752282783, relative change = 0.0009460125844278817
Cycle 61, free energy = 1377325.0276440016, relative change = 0.00123227844991071
Cycle 62, free energy = 1377320.6036971905, relative change = 0.0011525222469681224
Cycle 63, free energy = 1377315.622214722, relative change = 0.0012960888585295453
Cycle 64, free energy = 1377309.5813480178, relative change = 0.0015692544544170957
Cycle 65, free energy = 1377303.052055925, relative change = 0.0016932622143169826
Cycle 66, free energy = 1377297.357582863, relative change = 0.0014745883015729537
Cycle 67, free energy = 1377291.3708837442, relative change = 0.001547860882483115
Cycle 68, free energy = 1377284.819213344, relative change = 0.0016910696301523883
Cycle 69, free energy = 1377280.16629858, relative change = 0.0011995359865320297
Cycle 70, free energy = 1377276.875106017, relative change = 0.0008477604342421029
Cycle 71, free energy = 1377272.3309546842, relative change = 0.0011691348091558173
Cycle 72, free energy = 1377266.5474446379, relative change = 0.001485790403154828
Cycle 73, free energy = 1377260.0169396126, relative change = 0.001674884331186348
Cycle 74, free energy = 1377254.4773379925, relative change = 0.0014187308120508185
Cycle 75, free energy = 1377248.3056488875, relative change = 0.0015781184540896226
Cycle 76, free energy = 1377239.763802657, relative change = 0.00217941415041289
Cycle 77, free energy = 1377230.6284576755, relative change = 0.0023254225281021637
Cycle 78, free energy = 1377224.4334878002, relative change = 0.0015744606230363631
Cycle 79, free energy = 1377219.6282424927, relative change = 0.0012197704416303248
Cycle 80, free energy = 1377216.0356592361, relative change = 0.0009111156581707884
Cycle 81, free energy = 1377211.1164966065, relative change = 0.001245995240393448
Cycle 82, free energy = 1377207.3153851007, relative change = 0.0009618733122329076
Cycle 83, free energy = 1377204.966184519, relative change = 0.0005941132910514533
Cycle 84, free energy = 1377201.2942741648, relative change = 0.0009277653059648482
Cycle 85, free energy = 1377198.2779925617, relative change = 0.0007615301267484
Cycle 86, free energy = 1377195.6336285844, relative change = 0.0006671854774272866
Cycle 87, free energy = 1377191.3930790294, relative change = 0.0010687671260267413
Cycle 88, free energy = 1377185.188377184, relative change = 0.0015613606903841603
Cycle 89, free energy = 1377179.1073981924, relative change = 0.0015278888650981767
Cycle 90, free energy = 1377175.5107524558, relative change = 0.0009028667060582213
Cycle 91, free energy = 1377172.9599472224, relative change = 0.0006399193941603314
Cycle 92, free energy = 1377171.0211961465, relative change = 0.0004861371908250845
Cycle 93, free energy = 1377169.2835832324, relative change = 0.0004355125244825633
Cycle 94, free energy = 1377166.919232986, relative change = 0.0005922460049185309
Cycle 95, free energy = 1377164.4790553874, relative change = 0.0006108666093857384
Cycle 96, free energy = 1377163.215945759, relative change = 0.00031610305495521034
Cycle 97, free energy = 1377161.8636427387, relative change = 0.0003383099059908136
Cycle 98, free energy = 1377160.2806213445, relative change = 0.00039587267102307046
Cycle 99, free energy = 1377158.5352321828, relative change = 0.00043628621769622604
Cycle 100, free energy = 1377156.0877308918, relative change = 0.0006114156247563455
Finished training in 996.79s : active states = 4

(array([[4.34527663e-01, 5.59899855e-01, 5.55166292e-03, 2.08185023e-05],
        [4.04311441e-01, 5.82935691e-01, 1.27359859e-02, 1.68821646e-05],
        [7.56607855e-02, 9.08121522e-01, 1.61925823e-02, 2.51097459e-05],
        ...,
        [2.08614002e-09, 9.99997065e-01, 2.93303106e-06, 6.92658514e-13],
        [2.15182464e-08, 9.99977038e-01, 2.29406217e-05, 1.68292840e-10],
        [7.20794365e-05, 9.99848696e-01, 7.91689412e-05, 5.54241780e-08]]),
 array([[[4.00801315e-01, 2.60648715e-02, 7.64790279e-03, 1.35737916e-05],
         [3.22505286e-03, 5.55410163e-01, 1.26160195e-03, 3.03721924e-06],
         [2.82035607e-04, 1.44548664e-03, 3.82398896e-03, 1.51701608e-07],
         [3.03751913e-06, 1.51693705e-05, 2.49216046e-06, 1.19452190e-07]],

        [[7.55797801e-02, 3.16063756e-01, 1.26434338e-02, 2.44710349e-05],
         [5.26171126e-05, 5.82702149e-01, 1.80450761e-04, 4.73741020e-07],
         [2.83371183e-05, 9.33917357e-03, 3.36832946e-03, 1.45719140e-07],
         [5.12036218e-08, 1.64434083e-05, 3.68301810e-07, 1.92508740e-08]],

        [[4.81836534e-04, 7.42403206e-02, 9.37867537e-04, 7.60822275e-07],
         [2.22540849e-06, 9.08030397e-01, 8.88022173e-05, 9.77148341e-08],
         [1.19705426e-06, 1.45357565e-02, 1.65559871e-03, 3.00200584e-08],
         [2.10676341e-09, 2.49274565e-05, 1.76319889e-07, 3.86280084e-09]],

        ...,

        [[1.39362508e-13, 2.52299504e-08, 9.53022255e-12, 1.66417350e-18],
         [2.08581896e-09, 9.99994447e-01, 2.92419272e-06, 6.92622387e-13],
         [1.81696128e-13, 2.59238683e-06, 8.82881116e-09, 3.44597820e-17],
         [1.61467942e-19, 2.24480613e-12, 4.74774249e-16, 2.23893773e-21]],

        [[1.18512040e-13, 2.07987834e-09, 6.14313096e-12, 3.33329812e-17],
         [2.15160623e-08, 9.99974179e-01, 2.28645356e-05, 1.68283580e-10],
         [2.06558632e-12, 2.85694905e-06, 7.60799316e-08, 9.22718511e-15],
         [5.13101847e-19, 6.91514232e-13, 1.14360129e-15, 1.67578507e-19]],

        [[3.41869759e-09, 1.79205328e-08, 1.78921432e-10, 9.45654351e-14],
         [7.20251858e-05, 9.99827679e-01, 7.72784422e-05, 5.54017385e-08],
         [5.08316221e-08, 2.09994487e-05, 1.89031907e-06, 2.23316037e-11],
         [4.14692637e-13, 1.66931634e-10, 9.33193622e-13, 1.33199048e-14]]]),
 array([1381159.09516471, 1380383.24074971, 1380006.35042783,
        1379686.22787817, 1379381.68277004, 1379127.40598509,
        1378913.84227953, 1378748.60751508, 1378604.81254361,
        1378461.98060362, 1378334.72111959, 1378227.40811012,
        1378155.89542477, 1378094.93645894, 1378025.39375149,
        1377964.04834425, 1377909.29802348, 1377866.36879839,
        1377828.8088543 , 1377798.10140328, 1377768.28230602,
        1377730.61354666, 1377702.63803043, 1377679.08043139,
        1377657.04257999, 1377634.80699568, 1377615.71100526,
        1377600.9551945 , 1377587.92280063, 1377573.05522107,
        1377559.69287666, 1377549.02401169, 1377538.94762178,
        1377527.17754833, 1377516.80748345, 1377508.88802707,
        1377498.31993013, 1377485.77163038, 1377473.25146843,
        1377458.72116339, 1377445.36076598, 1377436.20922646,
        1377428.63834596, 1377421.98708793, 1377416.15515118,
        1377410.2474398 , 1377402.85374763, 1377394.46278205,
        1377389.7110516 , 1377385.07078034, 1377379.54021384,
        1377373.78011292, 1377368.12457461, 1377363.28297519,
        1377357.78195511, 1377351.23422371, 1377344.51503797,
        1377338.27400379, 1377333.37488934, 1377329.75228278,
        1377325.027644  , 1377320.60369719, 1377315.62221472,
        1377309.58134802, 1377303.05205592, 1377297.35758286,
        1377291.37088374, 1377284.81921334, 1377280.16629858,
        1377276.87510602, 1377272.33095468, 1377266.54744464,
        1377260.01693961, 1377254.47733799, 1377248.30564889,
        1377239.76380266, 1377230.62845768, 1377224.4334878 ,
        1377219.62824249, 1377216.03565924, 1377211.11649661,
        1377207.3153851 , 1377204.96618452, 1377201.29427416,
        1377198.27799256, 1377195.63362858, 1377191.39307903,
        1377185.18837718, 1377179.10739819, 1377175.51075246,
        1377172.95994722, 1377171.02119615, 1377169.28358323,
        1377166.91923299, 1377164.47905539, 1377163.21594576,
        1377161.86364274, 1377160.28062134, 1377158.53523218,
        1377156.08773089]))

Inspect model

Interaction between brain and physiological measures

Let’s start by retrieving the parameters describing the interaction, i.e., the \(\beta\) values, between the two sets of timeseries: the brain data and the physiological data (heart rate & respiration). The beta-values can be obtained from the trained model using the get_betas() function (or alternatively get_beta(k) to obtain only the beta-values for state k):

K = brainphys_glhmm.hyperparameters["K"] # the number of states
q = brain_data.shape[1] # the number of parcels/channels
state_betas = np.zeros(shape=(2,q,K))
state_betas = brainphys_glhmm.get_betas()

Since we here defined \(\beta\) to be time-varying, i.e., state-dependent, we have a matrix describing the interaction between each of the 50 brain regions and the 2 physiological measures for each of the 4 states:

cmap = "coolwarm"
ytick =["Heart rate", "Respiration"]
for k in range(K):
    plt.subplot(2,2,k+1)
    plt.imshow(state_betas[:,:,k], cmap=cmap,aspect='auto', interpolation='none')
    plt.colorbar()
    plt.ylabel('Physiological data')
    plt.yticks(np.arange(2), ytick)
    plt.xlabel('Brain area')
    plt.title(f"Betas for state #{k+1}")
plt.subplots_adjust(hspace=0.5, wspace=1)
plt.show()

State means and covariances

We have defined the model so that each state also has a mean (amplitude) and covariance (functional connectivity), as in the Standard Gaussian HMM. We can retrieve them by using the get_mean and get_covariance_matrix functions:

state_means = np.zeros(shape=(q, K))
for k in range(K):
    state_means[:,k] = brainphys_glhmm.get_mean(k) # the state means in the shape (no. features, no. states)
state_FC = np.zeros(shape=(q, q, K))
for k in range(K):
    state_FC[:,:,k] = brainphys_glhmm.get_covariance_matrix(k=k) # the state covariance matrices in the shape (no. features, no. features, no. states)

And plot them:

plt.imshow(state_means,cmap=cmap, interpolation="none")
plt.colorbar(label='Activation Level') # Label for color bar
plt.title("State mean activation")
plt.xticks(np.arange(K), np.arange(1,K+1))
plt.gca().set_xlabel('State')
plt.gca().set_ylabel('Brain region')
plt.tight_layout()  # Adjust layout for better spacing
plt.show()

for k in range(K):
    plt.subplot(2,2,k+1)
    plt.imshow(state_FC[:,:,k], cmap=cmap, interpolation="none")
    plt.xlabel('Brain region')
    plt.ylabel('Brain region')
    plt.colorbar()
    plt.title("State covariance\nstate #%s" % (k+1))
plt.subplots_adjust(hspace=0.7, wspace=0.8)
plt.show()

Dynamics: Transition probabilities and Viterbi path

We can also look at the transition probabilities and the Viterbi path to understand the temporal sequence in which the states occur. See Standard Gaussian HMM for a detailed explanation. The transition probabilities with and without self-transitions:

TP = brainphys_glhmm.P.copy() # the transition probability matrix

# Plot Transition Probabilities
plt.figure(figsize=(7, 4))

# Plot 1: Original Transition Probabilities
plt.subplot(1, 2, 1)
plt.imshow(TP, cmap=cmap, interpolation='nearest')  # Improved color mapping
plt.title('Transition Probabilities')
plt.xlabel('To State')
plt.ylabel('From State')
plt.colorbar(fraction=0.046, pad=0.04)

# Plot 2: Transition Probabilities without Self-Transitions
TP_noself = TP - np.diag(np.diag(TP))  # Remove self-transitions
TP_noself2 = TP_noself / TP_noself.sum(axis=1, keepdims=True)  # Normalize probabilities
plt.subplot(1, 2, 2)
plt.imshow(TP_noself2, cmap=cmap, interpolation='nearest')  # Improved color mapping
plt.title('Transition Probabilities\nwithout Self-Transitions')
plt.xlabel('To State')
plt.ylabel('From State')
plt.colorbar(fraction=0.046, pad=0.04)

plt.tight_layout()  # Adjust layout for better spacing
plt.show()

And the Viterbi path:

vpath = brainphys_glhmm.decode(X=phys_data, Y=brain_data, indices=T_t, viterbi=True)

And plot the Viterbi path (see also graphics module):

graphics.plot_vpath(vpath, title="Viterbi path")

We can also visualize the Viterbi path for the first subject from timepoints 0-1000

num_subject = 0
graphics.plot_vpath(vpath[T_t[num_subject,0]:T_t[num_subject,1],:], title="Viterbi path")

As for the standard HMM, there is a range of useful summary metrics that you can compute to describe the obtained patterns. Have a look at the Standard Gaussian HMM tutorial to see how to obtain them from your trained model. These can be used, e.g., for statistical testing (see Statistical testing tutorial) or prediction/machine learning (see Prediction tutorial).