Extended Dynamical Causal Modelling for Phase Coupling...

(1) Overview

Introduction

An oscillatory process is a widespread phenomenon in nature that occurs in physical and biological systems [2, 3, 4]. Oscillatory and rhythmic activities play a prominent role in the interaction between biological systems, especially in the communication between brain areas [5]. Therefore it is of particular interest to understand how the interaction between the elements of these oscillating systems, such as coupled physical oscillators, or between oscillatory activities in distant brain regions, occurs. In this context phase reduction is a useful method to analyze a large oscillating network by representing every oscillating system with one variable – the phase.

In the theory of synchronization the interaction between oscillating systems is analyzed by the model of weakly coupled phase oscillators [3, 6, 7]:

1

{\dot{φ}}_{i} = ω_{i} + Q_{i} (φ_{1}, \dots, φ_{N}), i = 1, \dots, N,

where φ_i is the phase of an oscillating system (an oscillator) i, ω_i is the natural frequency of the oscillator i, Q_i is the interaction function (coupling function) with other oscillators. For weak coupling the coupling function Q_i in Eq. (1) can be simplified as the sum of the pairwise coupling functions q_i,j:

2

{\dot{φ}}_{i} = ω_{i} + \sum_{j = 1, j \neq i}^{N} q_{i, j} (φ_{i}, φ_{j}) .

The pairwise coupling functions q_i,j, on the other hand, can be represented as 2D surfaces (Figure 1).

The model of two weakly coupled phase oscillators. A coupling function can be represented as a surface.

One can obtain the coupling between the phase oscillators in the model Eq. (2) by finding the coupling functions q_i,j directly from the experimental data [8]. This method was actively developed in the past decades and used to find the directionality of the couplings, the causal relations and to build dynamical models [9, 10, 11, 12, 13, 14, 15, 16, 17, 18].

Kralemann and colleagues [9, 10] have shown that the phase extracted from the experimental signals (e.g. by using Hilbert or Wavelet transformation) behaves differently than the phase described in the theoretical model Eq. (1). The difference appears if the limit cycle of the oscillating system, the signal which we are measuring, is not circular, and the oscillatory signal has a non-sinusoidal form. This leads to a non-linear growth of the measured phase, even if the oscillating system has no input, as demonstrated in Figure 2. Moreover, it was shown in [10] that this non-linear growth of the measured phase can cause spurious couplings in the system under investigation. The problem was resolved by introducing a transformation function between the observable and the theoretical phases (referred to as the proto- and the true phases in [10]). Kralemann and colleagues implemented this approach into their DAMOCO toolbox¹ by approximating the transformation function from the observable phase to the theoretical phase using the distribution of the observable phase over a long period of time.

The difference between the theoretical and the observable (measured) phase when the limit cycle of the oscillating system is not circular (left panel). The signal of the theoretical phase is a sine wave, whereas the one of the observable phase is not (center panel). The observable phase grows non-linearly (right panel).

The software tool presented in this metapaper – extended Dynamical Causal Modelling for Phase Coupling (eDCM PC) – addresses the problem of a potential non-linear growth in the observable phase. It is implemented as an extension to Dynamical Causal Modelling for phase coupling (DCM PC) [1]. eDCM PC aims at finding the effective coupling in a network of oscillating systems by inferring the coupling function from the measured signals. In particular, the software tool is designed to reconstruct the coupling functions for the cases when the phase is not uniformly distributed, i.e. in the case when the transformation between the observable and the theoretical phases should be taken into account. Moreover, eDCM PC extends DCM PC by allowing to find the coupling between different frequency bands, thereby making it possible to analyze n:m synchronization cases. In contrast to the DAMOCO toolbox, eDCM PC uncovers the transformation functions from the theoretical to the observable phases together with the coupling functions. It uses the capability of DCM, based on Bayesian inference, to analyze and compare several possible network structures at the same time.

The theoretical results related to eDCM PC and the numerical testing on synthetic data sets were presented in our previous work [19]. Since its first presentation in [19] eDCM PC was supplemented with a documentation, additional tests and plotting functions for data sets and results, and was rearranged in a user-friendly structure. In this work we give a detailed description on the architecture and usage of eDCM PC, as well as provide testing examples.

Implementation and architecture

The architecture of eDCM PC is presented in Figure 3 as a data flow with different interface levels: data preprocessing level, user-interface level, eDCM PC level, and SPM/DCM level. We describe the levels in bottom up direction, namely starting from the SPM/DCM level up to the eDCM PC level by indicating specific implementation details. Thereafter, in the usage part, we describe the user-interface level.

The data flow and the architecture scheme of eDCM PC. The green background area denotes the elements introduced in eDCM PC.

Dynamic causal modelling

The structure of the eDCM PC is imposed by the architecture of the Dynamic Causal Modelling (DCM) [20]. DCM is an opensource toolbox within the Statistical Parametric Mapping (SPM) software [21], developed to analyze connectivity in the brain. DCM can work with different modalities (functional magnetic resonance imaging (fMRI), EEG, MEG, local field potentials (LFP), functional near-infrared spectroscopy (fNIRS)) [20, 22, 23, 24, 25, 26, 27, 28, 29, 30]. The usage of a wide variety of modalities is provided by a particular feature of the DCM architecture that allows the modification of its different parts without changing the common, basic structure of the software [30].

DCM uses Bayesian inference to find the parameters of the system and the coupling between the brain regions. The elements of the common structure of DCM are the data preprocessing component (Figure 3a), the modelling component (Figure 3f), and the statistical component (Figure 3i). These elements can be modified with respect to the modality and the modelling, however the interaction between them remains the same for all versions of DCM.

Initial raw data are transformed into time courses of observables, e.g. into observable phases $\hat{θ}$ (Figure 3a). In the modelling component (Figure 3f) synthetic signals are generated by numerically integrating the model with an initial guess for the system parameters (priors). The modelling component consists of the evolution equations (Figure 3g), a system of ODEs that represents the hidden state of the system and the connections between parts of the system, and the observation equations (Figure 3h), functions that express the observables using the hidden states. In the statistical component of DCM (Figure 3i) the observables $\hat{θ}$ obtained from the data and the observables θ from the synthetic signals are “compared” by means of the Bayesian model inversion, namely Variational Laplace (VL) (Figure 3j). In a simplified manner, the result of this comparison gives the correction for the system parameters. Using these corrected parameter values the modelling component generates new synthetic signals. This procedure repeats until the required convergence between the measured and the synthetically generated observables is reached (Figure 3k).

Modelling component

eDCM PC introduces a new modelling component. As mentioned before the modelling component contains the evolution equation and the observation equation. The extended evolution equation is defined as

3

{\dot{φ}}_{i} = ω_{i} + \sum_{j = 1, j \neq i}^{N_{q}} q_{ij} (φ_{i}, φ_{j}),

where the pairwise coupling functions are constructed in the way to detect ‘n:m‘ synchronization

4

q_{ij} (φ_{i}, φ_{j}) = \sum_{n, m = - N_{q}, n, m \neq 0}^{N_{q}} Q_{ij} e^{i (n φ_{i} + m φ_{j})},

with complex Fourier coefficients Q_ij truncated at a given number of terms N_q. In the code, however, the coupling functions q_ij are represented with real-valued Fourier coefficients as

5

\begin{matrix} q_{ij} (φ_{i}, φ_{j}) & = \sum_{n = 1}^{N_{q}} \sum_{m = 1}^{N_{q}} [a_{ij}^{(n, m)} \cos (n φ_{i}) \cos (m φ_{j}) \\ + b_{ij}^{(n, m)} \cos (n φ_{i}) \sin (m φ_{j}) \\ + c_{ij}^{(n, m)} \sin (n φ_{i}) \cos (m φ_{j}) \\ + d_{ij}^{(n, m)} \sin (n φ_{i}) \sin (m φ_{j})] . \end{matrix}

To introduce the transformation function between the theoretical and observable phases (the forward transformation) we use the observation equation. Following [10] we define the transformation as follows

6

θ_{i} = Θ (φ_{i}) = φ_{i} + \sum_{k = - N_{ρ}, k \neq 0}^{N_{ρ}} (e^{i k φ_{i}} - 1),

also truncated at a given number of terms N_ρ. The forward transformation function is represented in the code with real-valued coefficients:

7

θ_{i} = Θ (φ_{i}) = φ_{i} + \sum_{k = 1}^{N_{ρ}} \frac{1}{k} [α_{i}^{(k)} \sin (n φ_{i}) - β_{i}^{(k)} \cos (n φ_{i}) + β_{i}^{(k)}] .

Thus, in eDCM PC the hidden states are the theoretical phases φ_i and the observables are the observable phases θ_i defined for every region and every frequency band.

The resulting reconstructed coupling and transformation functions are (1) the real-valued coefficients (matrices) $a_{ij}^{(n, m)}$ , $b_{ij}^{(n, m)}$ , $c_{ij}^{(n, m)}$ , $d_{ij}^{(n, m)}$ defined for every connected pair of regions and (2) the real-valued coefficients (vectors) $α_{i}^{(k)}$ and $β_{i}^{(k)}$ defined for every region, respectively.

Initial approximation of the observation equation

Good convergence is provided by a good initial guess of the parameters of the system. Therefore, the initial approximation of the transformation function, i.e. the observation equation, is essential for a successful reconstruction of the system parameters and the coupling functions. eDCM PC uses the approach presented in [10] to approximate the inverse transformation function φ = Φ(θ), as shown in Figure 4a–c.

Initial approximation of the inverse transformation function and updating the initial conditions.

The average distribution of the observable phases ${\hat{θ}}_{i}$ over a long period of time for a region i can be approximated by a 2π periodic function σ_i(θ) with the mean equal to 1. (Figure 4a,b), and be written as the following Fourier series [10]

8

σ_{i} (θ) = 1 + \sum_{k = 1}^{N_{σ}} {\hat{α}}_{i}^{(k)} \cos (k φ_{i}) + {\hat{β}}_{i}^{(k)} \sin (k φ_{i}),

where the sum is truncated at N_σ. The integral of σ_i(θ) over one period gives us the inverse transformation function

9

Φ_{i} (θ) = θ_{i} + \sum_{k = 1}^{N_{σ}} \frac{1}{k} [{\hat{α}}_{i}^{(k)} \sin (k θ_{i}) - {\hat{β}}_{i}^{(k)} \cos (k φ_{i}) + {\hat{β}}_{i}^{(k)}],

which is the inverse function of the forward transformation function Θ(φ) defined in Eq. (6) (Figure 4c). Note, that the truncation orders N_σ and N_ρ are, in general, different.

After obtaining the coefficients $\hat{α}$ and $\hat{β}$ by approximating the average distribution of the observable phases σ_i(θ), eDCM PC finds the initial approximation of the coefficients α and β of the forward transformation Θ(φ).

Updating the initial conditions

In every step of the Bayesian estimation the updated parameters α and β change the transformation function Θ(φ). Therefore the initial conditions ${\hat{θ}}_{0}$ of the observations should also be recalculated in every step to obtain new initial conditions for the evolution equation φ₀ (Figure 4d). For this the updated inverse transformation function Φ(θ) is needed. Therefore, eDCM PC recalculates the coefficients $\hat{α}$ and $\hat{β}$ of Eq. (9) using the updated coefficients α and β after every step of the Bayesian estimation.

Usage

The data flow of the eDCM PC is presented in Figure 3. The eDCM PC package works with time courses of the observable phases obtained from oscillatory signals. The observable phase can be extracted using the Hilbert or Wavelet transformation applied to the raw signals filtered within a specific frequency band (Figure 3a).

In the user-interface level the usage of eDCM PC comes down to the initialization of the elements of the DCM structure (Figure 3b) and calling the modcomp.spm_dcm_po(DCM) function (Figure 3d).

The code of eDCM PC contains an example script test.m that runs the reconstruction procedure using a pre-simulated data set. For details of the DCM structure and calling the modcomp.spm_dcm_po(DCM) function we refer to the documentation included in the repository of eDCM PC. Here, we discuss only essential parts of the DCM structure.

The observable phases $\hat{θ}$ should be stored in the DCM structure in the fields DCM.xY.yn, where n is the trial index (Figure 3b). All trials should have the same length with an equally distributed time step.

The supposed network structure of the system is defined in the binary adjacency matrix DCM.A (Figure 3c). DCM allows the definition of several possible network structures and finds the structure best fitting the data using Bayesian model comparison. In that case DCM.A is an array of matrix cells, and different structures should be defined in different cells. For details of the definition of the structure, please refer to the documentation of DCM in Chapters 41,42, and 43 of [21].

The elements of the modelling component are defined in the substructure DCM.M, such as the mean frequency ω_i (DCM.M.freq) and the frequency band (DCM.M.fb) for every region, the truncation orders of the Fourier series for the coupling functions N_q (DCM.M.Nq), for the forward transformation functions N_ρ (DCM.M.Nrho), and for the inverse transformation functions N_σ (DCM.M.Nsig).

The initialized DCM structure is used as an argument to call the modcomp.spm_dcm_po(DCM) function (Figure 3d). This function performs the initial approximation of the inverse transformation function Φ(θ), defines the priors for all parameters of the system (modcomp.spm_dcm_po_priors), and assigns the functions of the model in the substructure DCM.M such as:

state activities generation function modcomp.spm_gen_po;
evolution equation function modcomp.spm_fx_po;
observable equation function modcomp.spm_gx_po;
linear observation function modcomp.spm_lx_po, which is a dummy function for compatibility with other versions of DCM. The function returns the values without change.

Thereafter, modcomp.spm_dcm_po(DCM) calls DCM routines with the DCM structure as the input argument to start the Bayesian estimation of the parameters of the system (Figure 3e).

The result of the modcomp.spm_dcm_po(DCM) function is an extended DCM structure that contains the estimated values of the system parameters in the substructure DCM.Ep (Figure 3k), namely the matrices of the Fourier coefficients of the coupling functions $a_{ij}^{(n, m)}$ , $b_{ij}^{(n, m)}$ , $c_{ij}^{(n, m)}$ , $d_{ij}^{(n, m)}$ in DCM.Ep.anm, DCM.Ep.bnm, DCM.Ep.cnm, DCM.Ep.dnm respectively, of the forward transformations in DCM.Ep.rc and DCM.Ep.rs, which correspond to $α_{i}^{(k)}$ and $β_{i}^{(k)}$ respectively.

Installing and testing

eDCM PC is written in MATLAB and uses Dynamic Causal Modelling (DCM), which is included in the Statistical Parametric Mapping (SPM) package [31]. The SPM package can be downloaded here https://www.fil.ion.ucl.ac.uk/spm/software/download/. The minimum version of SPM used by eDCM PC is SPM12. The detailed SPM installation information is provided here https://www.fil.ion.ucl.ac.uk/spm/docs/installation/.

eDCM PC can either be downloaded or cloned from our GitLab repository https://gitlab.com/azayeld/edcmpc. Add the directory containing eDCM PC package into MATLAB using Set Path. Use Add Folder, but not Add with Subfolders…. At restarting MATLAB will run startup.m script that allow to work with eDCM PC as a MATLAB-package. In the command line the message: “eDCM PC package ist initialized.” will occur as confirmation.

For installation details follow the instruction given in the documentation of the eDCM PC package https://gitlab.com/azayeld/edcmpc/-/blob/master/Docs/Doc.md#installation.

After installation, eDCM PC can be tested using the test.m script. Go to the folder, where eDCM PC is copied or cloned, and run the testing script by running test in the command line or if eDCM PC folder is just added to MATLAB run:

startup
test

Quality control

The users can easily test whether the installation was correct by calling a test.m script.

The package includes three other example scripts, which are also intended to test the quality of the reconstruction.

In the first example, eDCM PC successfully reconstructs the coupling and the transformation functions from the signals of the synthetically simulated system of two weakly coupled phase oscillators. In this example the parameters of the coupling are chosen such that one coupling is zero, i.e. making the coupling unidirectional. Moreover, an artificial distortion (transformation) is added to the theoretical phases in order to simulate the effect of non-linear growing of the observable phase. Thus, the parameter values reconstructed by eDCM PC can be compared with the original parameter values used by the simulation.

The second example analyzes two uni-directionally coupled neural mass models, namely, Jansen and Rit models, which simulate EEG-like signals. In this case neither the coupling function nor the non-linear relation between the observable and the theoretical phases are known, which simulates the real case scenario. In this example, eDCM PC demonstrates successful reconstruction of the coupling functions and shows the presence of a non-linear relation between the observable and theoretical phases in the neural mass models.

The third example analyses the openly accessible Steady-State Visual Evoked Potential (SSVEP) dataset from [32], specifically examining connectivity in single-frequency SSVEPs. The provided scripts download the necessary files, extract the time courses corresponding to a given stimulation frequency and prepare the data set in a form suitable for working with the eDCM PC. The results of this example demonstrate the presence of non-linear coupling between periodic visual stimulation and the EEG signal from the occipital cortex. This example serves as a minimal script demonstrating the analysis pipeline for a dataset using eDCM PC.

The eDCM PC package includes several routines which allow testing of the data set before running the reconstruction procedure and examining the quality of the reconstruction afterwards:

Routines that allow the analysis of the spectral relation between two regions by plotting an approximated 2D Fourier spectrum of the observable phases. These routines can be used to analyze the data before the reconstruction procedure to set the truncation orders of the Fourier series (DCM.M.Nq and DCM.M.Nrho).
Routines to plot the distribution of the observable phases together with the initial approximation of the inverse transformation function and the final reconstructed inverse transformation function. These routines can be used to set the value of the parameter DCM.M.Nsig, and test the final reconstructed inverse function.
Routines to plot the reconstructed pairwise coupling functions together with the measured data points. These functions are useful to examine the convergence after the reconstruction procedure. Since the coupling functions q_ij(φ_i,φ_j) and the measured data points (the measured observable phases ${\hat{θ}}_{i}$ and ${\hat{θ}}_{j}$ ) are in different domains, there are two possibilities to compare them: (1) by representing the data points on the domain of the theoretical phases φ using the inverse transformation functions ${\hat{φ}}_{i} = Φ_{i} ({\hat{θ}}_{i})$ and ${\hat{φ}}_{j} = Φ_{j} ({\hat{θ}}_{j})$ , (2) or by projecting the coupling functions onto the domain of the observable phases using the forward transformation functions $q_{ij} (θ_{i}, θ_{j}) = q_{ij} (θ_{i} = Θ_{i} (φ_{i}), θ_{i} = Θ_{i} (φ_{i}))$ .

These functions are also included in the code of the examples. For more details we refer to the documentation of the eDCM PC package.

(2) Availability

Operating system

Windows (>=7 SP1, 32-bit or 64-bit),
Linux (>=Kernel 2.4.x or 2.6.x and >= glibc (glibc6) 2.3.4), Ubuntu 14.04 or 14.10
MacOS (>=X 10.9.5 (Mavericks))

Please refer to the system requirements of MATLAB R2015b, which can be found here: https://de.mathworks.com/support/requirements/previous-releases.html.

Programming language

The numerical computing language MATLAB.

Additional system requirements

Minimal requirements for MATLAB R2015b are any Intel or AMD x86 processor supporting the SSE2 instruction set, i.e processors Intel Core-2-Duo or AMD Athlon 64 (SSE2 required) and above, Disk space >=1024MB (500MB MATLAB and 350 MB SPM), RAM >=1024MB.

For better performance use processors Intel Core i5/i7 from 2012 or newer, Disk space >=3GB, RAM >=4GB.

Dependencies

eDCM PC needs MATLAB and SPM:

MATLAB (>=8.6, R2015b)
Statistical Parametric Mapping (SPM) (SPM12, >= Update 6225)

List of contributors

Please refer to the list of authors.

Software location

Archive

Name: extended Dynamic Causal Modelling for Phase Coupling (eDCM PC)
Persistent identifier: https://doi.org/10.5281/zenodo.17406441
Licence: GNU General Public License v3.0
Publisher: Azamat Yeldesbay
Version published: Version 1.0.0
Date published: 21/10/2025

Code repository

Name: GitLab
Persistent identifier: https://gitlab.com/azayeld/edcmpc/-/tags/v1.0.0
Licence: GNU General Public License v3.0 https://gitlab.com/azayeld/edcmpc/-/blob/master/LICENSE
Date published: 21/10/2025

Language

English

(3) Reuse potential

Despite being originally developed to analyze the connectivity between oscillating brain regions, eDCM PC makes no assumption about the origin of the oscillating system. Therefore, the reuse potential of eDCM PC is not restricted to brain signals, but it can also be used in other biological as well as physical and mechanical systems. In general, this software can be used to reconstruct the effective connectivity in any network of oscillating systems. The only requirement is the presence of the phase information extracted from the rhythmic signals for every element of the network.

Although the eDCM PC offers a framework for estimating effective connectivity between oscillating systems of any kind using their phase information, its applicability is limited by the weakly coupled oscillators model. Stronger interactions between oscillating systems may require the consideration of phase as well as amplitude information. Moreover, the computational performance of the eDCM PC is limited by the Bayesian inference algorithm underlying DCM. Despite its efficient algorithmic design, the software can lead to prolonged runtimes when processing large-scale computations or high-dimensional datasets. Further development could address this by implementing the algorithm provided in our recent work [33] on the reconstruction of amplitude-phase dynamics. This could also be implemented outside of the SPM (DCM) package.

The collection of codes of eDCM PC and the testing and example scripts, together with the documentation are uploaded in the GitLab repository https://gitlab.com/azayeld/edcmpc, where users can clone and modify the code. The modifications can also be included into the package by sending a merge request. Moreover, the GitLab repository has a build-in support mechanism in the form of an issue tracking system. In the issue tracking system the users of eDCM PC can post their suggestions and problems with the code by opening an issue, which is monitored and will be solved by the authors.

Notes

[1] The Data Analysis with Models Of Coupled Oscillators (DAMOCO) toolbox can be found here http://www.stat.physik.uni-potsdam.de/~mros/damoco.html.

Acknowledgements

We would like to thank Gereon Fink for useful discussions on this project.

Competing Interests

The authors have no competing interests to declare.

Author Contributions

Azamat Yeldesbay: Designing, implementing and testing the program code, writing the article and the documentation, creation of the figures.

Silvia Daun: Supervision, writing the article and the documentation, revising the text and the figures.

Extended Dynamical Causal Modelling for Phase Coupling (eDCM PC)

Full Article