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

DBS(; name, namespace=nothing, frequency=130.0, amplitude=2.5, pulse_width=0.066, 
    offset=0.0, start_time=0.0, smooth=1e-4)

Create a continuous deep brain stimulation (DBS) stimulus with regular pulses.

Arguments:

• name: Name given to ODESystem object within the blox
• namespace: Additional namespace above name if needed for inheritance
• frequency: Pulse frequency in Hz
• amplitude: Pulse amplitude in arbitrary units
• pulse_width: Duration of each pulse in ms
• offset: Baseline value of the signal between pulses
• start_time: Time delay before stimulation begins in ms
• smooth: Smoothing parameter for pulse transitions, set to 0 for sharp transitions

Returns a DBS stimulus blox that outputs square pulses with specified parameters.

Generic2dOscillator(name, namespace, ...)

The Generic2dOscillator model is a generic dynamic system with two state
variables. The dynamic equations of this model are composed of two ordinary
differential equations comprising two nullclines. The first nullcline is a
cubic function as it is found in most neuron and population models; the
second nullcline is arbitrarily configurable as a polynomial function up to
second order. The manipulation of the latter nullcline's parameters allows
to generate a wide range of different behaviours.

Equations:

```math
        \begin{align}
        \dot{V} &= d \, \tau (-f V^3 + e V^2 + g V + \alpha W + \gamma I) \\
        \dot{W} &= \dfrac{d}{	au}\,\,(c V^2 + b V - \beta W + a)
        \end{align}
```

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: FitzHugh, R., Impulses and physiological states in theoretical models of nerve membrane, Biophysical Journal 1: 445, 1961.

Nagumo et.al, An Active Pulse Transmission Line Simulating Nerve Axon, Proceedings of the IRE 50: 2061, 1962.

Stefanescu, R., Jirsa, V.K. Reduced representations of heterogeneous mixed neural networks with synaptic coupling. Physical Review E, 83, 2011.

Jirsa VK, Stefanescu R. Neural population modes capture biologically realistic large-scale network dynamics. Bulletin of Mathematical Biology, 2010.

Stefanescu, R., Jirsa, V.K. A low dimensional description of globally coupled heterogeneous neural networks of excitatory and inhibitory neurons. PLoS Computational Biology, 4(11), 2008).

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. Defaults to 1 for cortical regions, 14 for subcortical.
• H: See equation for use. Defaults to 0.02 for both cortical and subcortical regions.
• λ: See equation for use. Defaults to 5 for cortical regions, 400 for subcortical.
• r: See equation for use. Defaults to 0.15 for cortical regions, 0.1 for subcortical.
• 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

KuramotoOscillator(name, namespace, ...)

Simple implementation of the Kuramoto oscillator as described in the original paper [1].
Useful for general models of synchronization and oscillatory behavior.
The general form of the Kuramoto oscillator is given by:
Equations:

```math
        \begin{equation}
        \dot{\theta_i} = \omega_i + \frac{1}{N}\sum_{j=1}^N{K_{i, j}\text{sin}(\theta_j - \theta_i)}
        \end{equation}
```

Where this describes the connection between regions $i$ and $j$. An alternative form
which includes a noise term for each region is also provided, taking the form:

```math
        \begin{equation}
        \dot{\theta_i} = \omega_i + \zeta dW_i \frac{1}{N}\sum_{j=1}^N{K_{i, j}\text{sin}(\theta_j - \theta_i)}
        \end{equation}
```

where $W_i$ is a Wiener process and $\zeta_i$ is the noise strength.

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. Default parameter values are taken from [2].

Citations:

1. Kuramoto, Y. (1975). Self-entrainment of a population of coupled non-linear oscillators. In: Araki, H. (eds) International Symposium on Mathematical Problems in Theoretical Physics. Lecture Notes in Physics, vol 39. Springer, Berlin, Heidelberg. https://doi.org/10.1007/BFb0013365

2. Sermon JJ, Wiest C, Tan H, Denison T, Duchet B. Evoked resonant neural activity long-term dynamics can be reproduced by a computational model with vesicle depletion. Neurobiol Dis. 2024 Jun 14;199:106565. doi: 10.1016/j.nbd.2024.106565. Epub ahead of print. PMID: 38880431.

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)

PINGNeuronExci(name, namespace, C, g_Na, V_Na, g_K, V_K, g_L, V_L, I_ext, τ_R, τ_D)

Create an excitatory neuron from Borgers et al. (2008).
The formal definition of this blox is:

\[\frac{dV}{dt} = \frac{1}{C}(-g_{Na}*m_{\infty}^3*h*(V - V_{Na}) - g_K*n^4*(V - V_K) - g_L*(V - V_L) + I_{ext} + jcn) \m_{\infty} = \frac{a_m(V)}{a_m(V) + b_m(V)} \frac{dn}{dt} = a_n(V)*(1 - n) - b_n(V)*n \frac{dh}{dt} = a_h(V)*(1 - h) - b_h(V)*h \frac{ds}{dt} = \frac{1}{2}*(1 + \tanh(V/10))*(\frac{1 - s}{\tau_R} - \frac{s}{\tau_D})\]

where $jcn$ is any input to the blox. Note that this is a modified Hodgkin-Huxley formalism with an additional synaptic accumulation term. Synapses are added into the $jcn$ term by connecting the postsynaptic neuron's voltage to the presynaptic neuron's output:

\[jcn = w*s*(V_E - V)\]

where $w$ is the weight of the synapse and $V_E$ is the reversal potential of the excitatory synapse.

Inputs:

• name: Name given to ODESystem object within the blox.
• namespace: Additional namespace above name if needed for inheritance.
• C: Membrane capacitance (defaults to 1.0).
• g_Na: Sodium conductance (defaults to 100.0).
• V_Na: Sodium reversal potential (defaults to 50.0).
• g_K: Potassium conductance (defaults to 80.0).
• V_K: Potassium reversal potential (defaults to -100.0).
• g_L: Leak conductance (defaults to 0.1).
• V_L: Leak reversal potential (defaults to -67.0).
• I_ext: External current (defaults to 0.0).
• τ_R: Rise time of synaptic conductance (defaults to 0.2).
• τ_D: Decay time of synaptic conductance (defaults to 2.0).

PINGNeuronInhib(name, namespace, C, g_Na, V_Na, g_K, V_K, g_L, V_L, I_ext, τ_R, τ_D)

Create an inhibitory neuron from Borgers et al. (2008).
The formal definition of this blox is:

\[\frac{dV}{dt} = \frac{1}{C}(-g_{Na}*m_{\infty}^3*h*(V - V_{Na}) - g_K*n^4*(V - V_K) - g_L*(V - V_L) + I_{ext} + jcn) \m_{\infty} = \frac{a_m(V)}{a_m(V) + b_m(V)} \frac{dn}{dt} = a_n(V)*(1 - n) - b_n(V)*n \frac{dh}{dt} = a_h(V)*(1 - h) - b_h(V)*h \frac{ds}{dt} = \frac{1}{2}*(1 + \tanh(V/10))*(\frac{1 - s}{\tau_R} - \frac{s}{\tau_D})\]

where $jcn$ is any input to the blox. Note that this is a modified Hodgkin-Huxley formalism with an additional synaptic accumulation term. Synapses are added into the $jcn$ term by connecting the postsynaptic neuron's voltage to the presynaptic neuron's output:

\[jcn = w*s*(V_I - V)\]

where $w$ is the weight of the synapse and $V_I$ is the reversal potential of the inhibitory synapse.

Inputs:

• name: Name given to ODESystem object within the blox.
• namespace: Additional namespace above name if needed for inheritance.
• C: Membrane capacitance (defaults to 1.0).
• g_Na: Sodium conductance (defaults to 35.0).
• V_Na: Sodium reversal potential (defaults to 55.0).
• g_K: Potassium conductance (defaults to 9.0).
• V_K: Potassium reversal potential (defaults to -90.0).
• g_L: Leak conductance (defaults to 0.1).
• V_L: Leak reversal potential (defaults to -65.0).
• I_ext: External current (defaults to 0.0).
• τ_R: Rise time of synaptic conductance (defaults to 0.5).
• τ_D: Decay time of synaptic conductance (defaults to 10.0).

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.

ProtocolDBS(; name, namespace=nothing, frequency=130.0, amplitude=2.5,
              pulse_width=0.066, offset=0.0, start_time=0.0, smooth=1e-4,
              pulses_per_burst=10, bursts_per_block=12, 
              pre_block_time=200.0, inter_burst_time=200.0)

Create a deep brain stimulation (DBS) stimulus consisting of a block of pulse bursts.

Arguments:

• name: Name given to ODESystem object within the blox
• namespace: Additional namespace above name if needed for inheritance
• frequency: Pulse frequency in Hz
• amplitude: Pulse amplitude in arbitrary units
• pulse_width: Duration of each pulse in ms
• offset: Baseline value of the signal between pulses
• start_time: Time delay before stimulation begins in ms
• smooth: Smoothing parameter for pulse transitions, set to 0 for sharp transitions
• pulsesperburst: Number of pulses in each burst
• burstsperblock: Number of bursts in the block
• preblocktime: Time before the block starts in ms
• interbursttime: Time between bursts in ms

Returns a DBS stimulus blox that outputs a block of pulse bursts.

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/unknowns of the system that are neither external inputs nor measurements, i.e. these are the dynamic states
- `idx`: indices of these states

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 `Connector` object with all relevant connections 
from lower levels and this level.

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 Connector),
that will later be composed and will automatically get the outermost namespace.

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 delay and weight into a single vector

function paramscoping(;tunable=true, 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, indices)

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)
- `indices`  : 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

system_from_graph(g::MetaDiGraph, p=Num[]; name, simplify=true, graphdynamics=false, kwargs...)

Take in a MetaDiGraph g describing a network of neural structures (and optionally a vector of extra parameters p) and construct a System which can be used to construct various Problem types (i.e. ODEProblem) for use with DifferentialEquations.jl solvers.

If simplify is set to true (the default), then the resulting system will have structural_simplify called on it with any remaining keyword arguments forwared to structural_simplify. That is,

@named sys = system_from_graph(g; kwarg1=x, kwarg2=y)

is equivalent to

@named sys = system_from_graph(g; simplify=false)
sys = structural_simplify(sys; kwarg1=x, kwarg2=y)

See the docstring for structural_simplify for information on which options it supports.

If graphdynamics=true (defaults to false), the output will be a GraphSystem from GraphDynamics.jl, and the kwargs will be sent to the GraphDynamics constructor instead of using ModelingToolkit.jl. The GraphDynamics.jl backend is typically significantly faster for large neural systems than the default backend, but is experimental and does not yet support all Neuroblox.jl features.

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)