Dynamic coupling
================

The nnsa package supports a general framework for computation of dynamic coupling (coupling as function of time)
between two signals. nnsa contains easy-to-use code to apply the same framework to any signals.
Examples of use cases are autoregulation, neurovascular coupling and multimodal graphs.

This dynamic coupling framework takes two or more signals, divides them into segments and
computes the coupling between the signals. This is done per segment, such that a time course of coupling values
is obtained (hence the term dynamic coupling).

Additionally, through Monte Carlo simulations using surrogate signals, p-values can be obtained for the coupling values.
It supports a range of different coupling function, such as correlation, coherence, RBF kernel, transfer
function gain, mutual information. Additionally, you can define your own coupling function.

This script demonstrates how the code can be used to compute dynamic coupling.

Author: Tim Hermans (tim-hermans@hotmail.com).

Link to script: `dynamic_coupling.py <https://gitlab.com/timhermans/nnsa_public/-/blob/main/examples/dynamic_coupling.py>`_



.. code-block:: python

    import numpy as np
    from matplotlib import pyplot as plt
    
    from nnsa import DynamicCoupling
    from nnsa.utils.dummy_data import generate_timeseries
    
    plt.close('all')


Generate two random signals that exhibit transient coupling 



.. code-block:: python

    
    # Generate random walks.
    fs = 100
    x = generate_timeseries(size=int(60*10*fs), demean=True, seed=43)  # 10 minutes.
    y = generate_timeseries(size=len(x), demean=True, seed=10)
    t = np.arange(len(x))/fs
    
    # Introduce coupling in the middle half of the time period by mixing the two random signals there.
    idx = len(x)//4
    gain = 4
    y[idx:-idx] = y[idx:-idx] + gain*x[idx:-idx]
    
    # Plot the signals.
    fig, ax = plt.subplots(tight_layout=True)
    ax.plot(t, x, label='x')
    ax.plot(t, y, label='y')
    ax.legend()
    ax.set_xlabel('Time (s)')
    ax.set_ylabel('Signal (a.u.)')
    ax.set_title('Test signals')
    

.. figure:: ../examples/figs/dynamic_coupling_1.png


The dynamic coupling can be computed using the DynamicCoupling class in nnsa. The options for the coupling analysis can be specified when initiating the class. The options are described in the code for default_parameters(). We can print the default parameters: 



.. code-block:: python

    print(DynamicCoupling.default_parameters())
    
    # Create an instance of the DynamicCoupling class with custom parameters, overruling some defaults:
    dc = DynamicCoupling(
        method='correlation',  # Choose from 'correlation', 'coherence', 'transfer_function'/'TF'
        method_kwargs={},
        segmentation={'segment_length': 30, 'overlap': 15},  # In seconds.
        surrogates=dict(n_surrogates=50, how='IAAFT', per_segment=True, seed=43)
    )


Run the algorithm on the test data. 



.. code-block:: python

    result = dc.process(x, y, fs=fs, verbose=1)
    
    # The result contains a coupling value for each segment:
    print('shape result.coupling :', result.coupling.shape)
    
    # If surrogates were computed, these are stored in the result.coupling_surrogates attribute:
    print('shape result.coupling_surrogates :', result.coupling_surrogates.shape)
    
    # Plot the result.
    plt.figure(tight_layout=True)
    result.plot()
    plt.legend(loc='upper left', bbox_to_anchor=(1, 1))
    

.. figure:: ../examples/figs/dynamic_coupling_2.png


Compare different coupling functions. 



.. code-block:: python

    fig, axes = plt.subplots(3, 1, tight_layout=True, sharex='all')
    for ax, method in zip(axes, ['correlation', 'coherence', 'TF']):
        print(f'Computing with {method}...')
        dc = DynamicCoupling(
            method=method,
            method_kwargs={},
            segmentation={'segment_length': 30, 'overlap': 15},
            surrogates=dict(n_surrogates=50, how='IAAFT', per_segment=False, seed=43)  # Do not compute surrogate coupling values per segment to speed computations up.
        )
        result = dc.process(x, y, fs=fs, verbose=1)
        result.plot(ax=ax)
        plt.legend(loc='upper left', bbox_to_anchor=(1, 1))
        ax.set_title(method)
    

.. figure:: ../examples/figs/dynamic_coupling_3.png


You can also create your own coupling function. The function should accept two inputs (x, y) and return a single coupling value (should be high for high coupling, low for low coupling). For example, to compute coupling with a Spearman correlation coefficient, we could do the following: 



.. code-block:: python

    
    def coupling_fun(x, y):
        # Spearman correlation.
        from scipy.stats import spearmanr
        r, pval = spearmanr(x, y)
        return abs(r)
    
    
    # Pass the custom function to DynamicCoupling().
    dc = DynamicCoupling(
        method=coupling_fun,
        method_kwargs={},
        segmentation={'segment_length': 30, 'overlap': 15},  # In seconds.
        surrogates=dict(n_surrogates=50, how='IAAFT', per_segment=False, seed=43)
    )
    
    # Process. Here we can also give a label to the result (which shows up in the plot below).
    result = dc.process(x, y, fs=fs, label='Spearman', verbose=1)
    
    # Plot the result.
    plt.figure(tight_layout=True)
    result.plot()
    plt.legend(loc='upper left', bbox_to_anchor=(1, 1))
    

.. figure:: ../examples/figs/dynamic_coupling_4.png