API Documentation

Arguments:

• name: Name given to ODESystem object within the blox.
• namespace: Additional namespace above name if needed for inheritance.
• lnκ: logarithmic prefactor to signal decay H[1], set to 0 for standard parameter value.
• lnτ: logarithmic prefactor to transit time H[3], set to 0 for standard parameter value.
• lnϵ: logarithm of ratio of intra- to extra-vascular signal

NB: the prefix ln of the variables u, ν, q as well as the parameters κ, τ denotes their transformation into logarithmic space to enforce their positivity. This transformation is considered in the derivates of the model equations below.

Citations:

1. Stephan K E, Weiskopf N, Drysdale P M, Robinson P A, and Friston K J. Comparing Hemodynamic Models with DCM. NeuroImage 38, no. 3 (2007): 387–401. doi: 10.1016/j.neuroimage.2007.07.040
2. Hofmann D, Chesebro A G, Rackauckas C, Mujica-Parodi L R, Friston K J, Edelman A, and Strey H H. Leveraging Julia's Automated Differentiation and Symbolic Computation to Increase Spectral DCM Flexibility and Speed, 2023. doi: 10.1101/2023.10.27.564407

HarmonicOscillator(name, namespace, ω, ζ, k, h)

Create a harmonic oscillator blox with the specified parameters.
The formal definition of this blox is:

\[\frac{dx}{dt} = y-(2*\omega*\zeta*x)+ k*(2/\pi)*(atan((\sum{jcn})/h) \frac{dy}{dt} = -(\omega^2)*x\]

where ``jcn`` is any input to the blox.

Arguments:

• name: Name given to ODESystem object within the blox.
• namespace: Additional namespace above name if needed for inheritance.
• ω: Base frequency. Note the default value is scaled to give oscillations in milliseconds to match other blocks.
• ζ: Damping ratio.
• k: Gain.
• h: Threshold.

JansenRit(name, namespace, τ, H, λ, r, cortical, delayed)

Create a Jansen Rit blox as described in Liu et al.
The formal definition of this blox is:

\[\frac{dx}{dt} = y-\frac{2}{\tau}x \frac{dy}{dt} = -\frac{x}{\tau^2} + \frac{H}{\tau} [\frac{2\lambda}{1+\text{exp}(-r*\sum{jcn})} - \lambda]\]

where $jcn$ is any input to the blox.

Arguments:

• name: Name given to ODESystem object within the blox.
• namespace: Additional namespace above name if needed for inheritance.
• τ: Time constant. This is changed from the original source as the time constant was in seconds, while all our blocks are in milliseconds.
• H: See equation for use.
• λ: See equation for use.
• r: See equation for use.
• cortical: Boolean to determine whether to use cortical or subcortical parameters. Specifying any of the parameters above will override this.
• delayed: Boolean to indicate whether states are delayed

Citations:

1. Liu C, Zhou C, Wang J, Fietkiewicz C, Loparo KA. The role of coupling connections in a model of the cortico-basal ganglia-thalamocortical neural loop for the generation of beta oscillations. Neural Netw. 2020 Mar;123:381-392. doi: 10.1016/j.neunet.2019.12.021.

Jansen-Rit model block for canonical micro circuit, analogous to the implementation in SPM12

LarterBreakspear(name, namespace, ...)

Create a Larter Breakspear blox described in Endo et al. For a full list of the parameters used see the reference.
If you need to modify the parameters, see Chesebro et al. and van Nieuwenhuizen et al. for physiological ranges.

Arguments:

• name: Name given to ODESystem object within the blox.
• namespace: Additional namespace above name if needed for inheritance.
• Other parameters: See reference for full list. Note that parameters are scaled so that units of time are in milliseconds.

Citations:

1. Endo H, Hiroe N, Yamashita O. Evaluation of Resting Spatio-Temporal Dynamics of a Neural Mass Model Using Resting fMRI Connectivity and EEG Microstates. Front Comput Neurosci. 2020 Jan 17;13:91. doi: 10.3389/fncom.2019.00091.
2. Chesebro AG, Mujica-Parodi LR, Weistuch C. Ion gradient-driven bifurcations of a multi-scale neuronal model. Chaos Solitons Fractals. 2023 Feb;167:113120. doi: 10.1016/j.chaos.2023.113120.
3. van Nieuwenhuizen, H, Chesebro, AG, Polis, C, Clarke, K, Strey, HH, Weistuch, C, Mujica-Parodi, LR. Ketosis regulates K+ ion channels, strengthening brain-wide signaling disrupted by age. Preprint. bioRxiv 2023.05.10.540257; doi: https://doi.org/10.1101/2023.05.10.540257.

Ornstein-Uhlenbeck process Blox

variables: x(t): value jcn: input parameters: τ: relaxation time μ: average value σ: random noise (variance of OU process is τ*σ^2/2) returns: an ODE System (but with brownian parameters)

Ornstein-Uhlenbeck Coupling Blox This blox takes an input and multiplies that input with a OU process of mean μ and variance τ*σ^2/2

This blox allows to create edges that have fluctuating weights

variables: x(t): value jcn: input parameters: τ: relaxation time μ: average value σ: random noise (variance of OU process is τ*σ^2/2) returns: an ODE System (but with brownian parameters)

Subcortical blox
all subcprtical blox used in cortico-striatal model are defined here

WilsonCowan(name, namespace, τ_E, τ_I, a_E, a_I, c_EE, c_IE, c_EI, c_II, θ_E, θ_I, η)

Create a standard Wilson Cowan blox.
The formal definition of this blox is:

\[\frac{dE}{dt} = \frac{-E}{\tau_E} + \frac{1}{1 + \text{exp}(-a_E*(c_{EE}*E - c_{IE}*I - \theta_E + \eta*(\sum{jcn}))} \frac{dI}{dt} = \frac{-I}{\tau_I} + \frac{1}{1 + exp(-a_I*(c_{EI}*E - c_{II}*I - \theta_I)}\]

where $jcn$ is any input to the blox.

Arguments:

• name: Name given to ODESystem object within the blox.
• namespace: Additional namespace above name if needed for inheritance.
• Others: See equation for use.

WinnerTakeAllBlox

Creates a winner-take-all local circuit found in neocortex, typically 5 pyramidal (excitatory) neurons send synapses to a single interneuron (inhibitory) and receive feedback inhibition from that interneuron.

function LinearAlgebra.eigen(M::Matrix{Dual{T, P, np}}) where {T, P, np}

Dispatch of LinearAlgebra.eigen for dual matrices with complex numbers. Make the eigenvalue decomposition 
amenable to automatic differentiation. To do so compute the analytical derivative of eigenvalues
and eigenvectors. 

Arguments:
- `M`: matrix of type Dual of which to compute the eigenvalue decomposition. 

Returns:
- `Eigen(evals, evecs)`: eigenvalue decomposition returned as type LinearAlgebra.Eigen

ARVTarget Time series data is bandpass filtered and then the power spectrum is computed for a given time interval (control bin), returned as the average value of the power spectral density within a certain frequency band ([lb, ub]).

CDVTarget Time series data is bandpass filtered and hilbert-transformed. Phase angle is computed in radians. Circular difference is quantified as the angle of circular_location.

ControlError Returns the control error (deviation of the actual value from the target value).

PDVTarget Time series data is bandpass filtered and hilbert-transformed. Phase angle is computed in radians. Phase deviation is quantified as the angle difference between a given set of signals.

PLVTarget Time series data is bandpass filtered and hilbert-transformed. Phase angle is computed in radians.

function addnontunableparams(param, model)

Function adds parameters of a model that were not marked as tunable to a list of tunable parameters
and respects the MTK ordering of parameters.

Arguments:
- `paramlist`: parameters of an MTK system that were tagged as tunable
- `sys`: MTK system

Returns:
- `completeparamlist`: complete parameter list of a system, including those that were not tagged as tunable

bandpassfilter takes in time series data and bandpass filters it. It has the following inputs: data: time series data lb: minimum cut-off frequency ub: maximum cut-off frequency fs: sampling frequency order: filter order

Arguments:

• name: Name given to ODESystem object within the blox.
• lnϵ : logarithm of ratio of intra- to extra-vascular signal

NB: the prefix ln of the variables ν, q as well as the parameters ϵ denotes their transformation into logarithmic space to enforce their positivity.

Citations:

1. Stephan K E, Weiskopf N, Drysdale P M, Robinson P A, and Friston K J. Comparing Hemodynamic Models with DCM. NeuroImage 38, no. 3 (2007): 387–401. doi: 10.1016/j.neuroimage.2007.07.040
2. Hofmann D, Chesebro A G, Rackauckas C, Mujica-Parodi L R, Friston K J, Edelman A, and Strey H H. Leveraging Julia's Automated Differentiation and Symbolic Computation to Increase Spectral DCM Flexibility and Speed, 2023. doi: 10.1101/2023.10.27.564407

complexwavelet creates a complex morlet wavelet by windowing a complex sine wave with a Gaussian taper. The morlet wavelet is a special case of a bandpass filter in which the frequency response is Gaussian-shaped. Convolution with a complex wavelet is equivalent to performing a Hilbert transform of a bandpass filtered signal.

It has the following inputs: data: time series data dt : data sampling rate lb : lower bound wavelet frequency (in Hz) ub : upper bound wavelet frequency (in Hz) a : amplitude of the Gaussian taper, default is 1 n : number of wavelet cycles of the Gaussian taper, defines the trade-off between temporal precision and frequency precision larger n gives better frequency precision at the cost of temporal precision default is 6 Hz m : x-axis offset, default is 0 num_wavelets : number of wavelets to create, default is 5

And outputs: complex_wavelet : a family of complex morlet wavelets

This function converts a cross-spectral density (CSD) into a multivariate auto-regression (MAR) model. It first transforms the CSD into its cross-correlation function (Wiener-Kinchine theorem) and then computes the MAR model coefficients. csd : cross-spectral density matrix of size MxN; M: number of samples, N: number of cross-spectral dimensions (number of variables squared) w : frequencies dt : time step size p : number of time steps of auto-regressive model

This function returns coeff : array of length p of coefficient matrices of size sqrt(N)xsqrt(N) noise_cov : noise covariance matrix

This function implements equation 2 of the spectral DCM paper, Friston et al. 2014 "A DCM for resting state fMRI".
Note that nomenclature is taken from SPM12 code and it does not seem to coincide with the spectral DCM paper's nomenclature. 
For instance, Gu should represent the spectral component due to external input according to the paper. However, in the code this represents
the hidden state fluctuations (which are called Gν in the paper).
Gn in the code corresponds to Ge in the paper, i.e. the observation noise. In the code global and local components are defined, no such distinction
is discussed in the paper. In fact the parameter γ, corresponding to local component is not present in the paper.

function get_dynamic_states(sys)

Function extracts states from the system that are dynamic variables, 
get also indices of external inputs (u(t)) and measurements (like bold(t))
Arguments:
- `sys`: MTK system

Returns:
- `sts`  : states of the system that are neither external inputs nor measurements, i.e. these are the dynamic states
- `idx_u`: indices of states that represent external inputs
- `idx_m`: indices of states that represent measurements

Plain implementation of idft because AD dispatch versions for ifft don't work still!

Returns the complete namespace EXCLUDING the outermost (highest) level.
This is useful for manually preparing equations (e.g. connections, see BloxConnector),
that will later be composed and will automatically get the outermost namespace.

Returns the equations for all input variables of a system, 
assuming they have a form like : `sys.input_variable ~ ...`
so only the input appears on the LHS.

Input equations are namespaced by the inner namespace of blox
and then they are returned. This way during system `compose` downstream,
the higher-level namespaces will be added to them.

If blox isa AbstractComponent, it is assumed that it contains a `connector` field,
which holds a `BloxConnector` object with all relevant connections 
from lower levels and this level.

This function computes learning rate. It has the following inputs: outcomes: vector of 1's and 0's for behavioral outcomes windows: number of windows to split the outcome data into And the following outputs: rate: the learning rate across each window

This function converts multivariate auto-regression (MAR) model parameters to a cross-spectral density (CSD). A : coefficients of MAR model, array of length p, each element contains the regression coefficients for that particular time-lag. Σ : noise covariance matrix of MAR p : number of time lags freqs : frequencies at which to evaluate the CSD sf : sampling frequency

This function returns: csd : cross-spectral density matrix of size MxN; M: number of samples, N: number of cross-spectral dimensions (number of variables squared)

Maximum likelihood estimator of a multivariate, or vector auto-regressive model. y : MxN Data matrix where M is number of samples and N is number of dimensions p : time lag parameter, also called order of MAR model return values mar["A"] : model parameters is a NxNxP tensor, i.e. one NxN parameter matrix for each time bin k ∈ {1,...,p} mar["Σ"] : noise covariance matrix

function matlab_norm(A, p)

Simple helper function to implement the norm of a matrix that is equivalent to the one given in MATLAB for order=1, 2, Inf. 
This is needed for the reproduction of the exact same results of SPM12.

Arguments:
- `A`: matrix
- `p`: order of norm

Helper to merge delays and weights into a single vector

function paramscoping(;kwargs...)

Scope arguments that are already a symbolic model parameter thereby keep the correct namespace 
and make those that are not yet symbolic a symbol.
Keyword arguments are used, because parameter definition require names, not just values.

phasecosblox is creating a cos with angular frequency ω and variable phase phaseinter has the following parameters: ω: angular frequency t: time phaseinter: a function that returns phase as a function of time and returns: the resulting value

Usage: phaseint = phaseinter(0:0.1:50,phasedata) phaseout(t) = phasecosblox(0.1,t,phaseint) which is now a function of time and can be used in an input blox you can also use the dot operator to calculate time-series signal = phaseout.(collect(0:0.01:50))

phaseinter is creating a function that interpolates the phase data for any time given phaseinter has the following parameters: phaserange: a range, e.g. 0:0.1:50 which should reflect the time points of the data phasedata: phase at equidistant time points and returns: an function that returns an interpolated phase for t in range

phasesinblox is creating a sin with angular frequency ω and variable phase phaseinter has the following parameters: ω: angular frequency t: time phaseinter: a function that returns phase as a function of time and returns: the resulting value

Usage: phaseint = phaseinter(0:0.1:50,phasedata) phaseout(t) = phasesinblox(0.1,t,phaseint) which is now a function of time and can be used in an input blox you can also use the dot operator to calculate time-series signal = phaseout.(collect(0:0.01:50))

phaseangle takes in time series data, hilbert transforms it, and estimates the phase angle.

random_initials creates a vector of random initial conditions for an ODESystem that is composed of a list of blox. The function finds the initial conditions in the blox and then sets a random value in between range tuple given for that state.

It has the following inputs: odesys: ODESystem blox : list of blox

And outputs: u0 : Float64 vector of initial conditions

Non-symbolic, time-block-based way of `@register_symbolic sample_poisson(λ)`.

function setup_sDCM(data, stateevolutionmodel, initcond, csdsetup, priors, hyperpriors, params_idx)

Interface function to performs variational inference to fit model parameters to empirical cross spectral density.
The current implementation provides a Variational Laplace fit (see function above `variationalbayes`).

Arguments:
- `data`        : dataframe with column names corresponding to the regions of measurement.
- `model`       : MTK model, including state evolution and measurement.
- `initcond`    : dictionary of initial conditions, numerical values for all states
- `csdsetup`    : dictionary of parameters required for the computation of the cross spectral density
-- `dt`         : sampling interval
-- `freq`       : frequencies at which to evaluate the CSD
-- `p`          : order parameter of the multivariate autoregression model
- `priors`      : dataframe of parameters with the following columns:
-- `name`       : corresponds to MTK model name
-- `mean`       : corresponds to prior mean value
-- `variance`   : corresponds to the prior variances
- `hyperpriors` : dataframe of parameters with the following columns:
-- `Πλ_pr`      : prior precision matrix for λ hyperparameter(s)
-- `μλ_pr`      : prior mean(s) for λ hyperparameter(s)
- `params_idx`  : indices to separate model parameters from other parameters. Needed for the computation of AD gradient.

function spm_logdet(M)

SPM12 style implementation of the logarithm of the determinant of a matrix.

Arguments:
- `M`: matrix

vecparam(param::OrderedDict)

Function to flatten an ordered dictionary of model parameters and return a simple list of parameter values.

Arguments:
- `param`: dictionary of model parameters (may contain numbers and lists of numbers)