Dynamical system analysis

We first performed the numerical analysis of the Reduced Wong Wang model (Wong & Wang, 2006) using PyDSTool (Clewley et al., 2012). We use this to model a single frontostriatal circuit.

Illustration of the single frontostriatal circuit.

Script overview

The main part of the dynamical system analysis code is in OCD_modeling.analysis.rww_dst_analysis.py. The script can be run with the following options:

usage: python rww_dst_analysis.py [-h] [--create_model] [--compute_nullclines]
                                  [--compute_epc] [--save_model]
                                  [--run_stability_analysis]
                                  [--load_stability_analysis] [--plot_figs]
                                  [--plot_phasespace_grid]
                                  [--plot_bifurcation_diagrams] [--save_figs]
                                  [--save_outputs] [--timeout TIMEOUT]
                                  [--n_jobs N_JOBS] [--n_trajs N_TRAJS]
                                  [--n_op N_OP] [--load_sample_rww]
                                  [--plot_timeseries_phasespace_bif]

Named Arguments

--create_model

create symbolic model

--compute_nullclines

compute nullclines numerically

--compute_epc

compute equilibrium point curves numerically

--save_model

save symbolic model with its computed attributes

--run_stability_analysis

run stability analysis: find fixed points (semi-analytically) and perform linear stability analysis around them

--load_stability_analysis

load previously completed stability analysis

--plot_figs

plot figures

--plot_phasespace_grid

plot grid of phase spaces using discretized order parameters

--plot_bifurcation_diagrams

plot grid of bifurcation diagrams using discretized order parameters

--save_figs

save figures

--save_outputs

save analysis outputs

--timeout

timeout of the stability analysis (per parameter combination invoked)

Default: 3000

--n_jobs

number of processes used in parallelization

Default: 20

--n_trajs

number of trajectories (traces) to compute for phase space projections

Default: 10

--n_op

number of values taken for each order parameters

Default: 5

--load_sample_rww

load a sample of ReducedWongWang model object previously ran for illustration

--plot_timeseries_phasespace_bif

plot neat figure of timeseries, phase space and bifurcations with trajectories (paper quality)

Reduced Wong-Wang Model

First, we instanciate the Reduced Wong-Wang model (Wong & Wang, 2006) in PyDSTool (Clewley et al., 2012) following the implementation provided in Deco et al. (Deco et al., 2013):

OCD_modeling.analysis.create_model(params, args=None)

Create the Dynamical System in PyDSTool

\[\dot{S_i} = - \cfrac{S_i}{\tau_S} + (1 - S_i) \gamma H(x_i) + \sigma v_i\]

with

\[H(x) = \cfrac{ax-b}{1-\exp{(-d(ax-b))}}\]

and

\[x_i = w J_N S_i + G J_N \sum_j{C_{ij} S_j} + I_0\]

Parameters:

• params (dict) – Parameters of the model

• args (Argparse.Namespace) – Optional arguments

Returns:

rww – PyDSTool object containing the dynamical system

Return type:

PyDSTool.Vode_ODEsystem

Default model’s paremeters

Parameter	Value	(unit) Interpretation
\(a\)	270	(n/C) slope of the transfer function
\(b\)	108	(Hz) offset of the transfer function
\(d\)	0.154	decay of the transfer function
\(G\)	2.5	(N/A) global coupling gain of the system
\(J_N\)	0.2609	(nA) synaptic scaling factor
\(I_0\)	0.3	(nA) external input
\(\tau_S\)	0.1	timescale of the local population activity
\(\gamma\)	0.000641	(Hz) timescale of the coupled activity
\(w\)	0.9	(N/A) recurrent excitation factor
\(\sigma\)	0.01	(nA) noise amplitude

Stability analysis

The stability analysis involves 3 major steps, which are executed for each single discretization of the order parameters:

1. Compute the nullclines (where each state variable is at equilibrium) and finding the fixed points (where nullclines cross).

2. Compute trajectories to illustrate the system behavior in the phase space (optional).

3. Compute equilibrium point curves. These are the branches of the bifurcation diagram indicating stability regimes and transition types.

OCD_modeling.analysis.stability_analysis(order_params, default_params, out_queue, args, pdomain={'C_12': [-0.5, 1.5]})

Create model and analyse dynamics using PyDSTool.

Parameters:

• order_params (dict) – Fixed parameters, for which to analyse the system using discretized values, for example {'C_21': np.linspace(0.2,0.8,4)}.

• default_params (dict) – Default model’s parameters.

• out_queue (Queue) – Queue to put results in (used for parallel computation) for each values of order parameter.

• args (Argparse.Namespace) – Structure of necessary options. For example, must include args.compute_epc = True to compute equilibrium point curves.

• pdomain (dict) – Free variable (or order parameter) to perform the bifurcation analyis from.

Returns:

A dict with model, nullclines (ncs), fixed points (fps), trajectories (trajs), and a pickled (dilled) continuation object (if equilibrium curves are asked in args), is appended to the queue.

Return type:

None

OCD_modeling.analysis.get_fixed_points(model, params, xdomain={'S1': [0, 1], 'S2': [0, 1]}, args=None)

Get model’s nullclines \(\frac{dS_1}{dt}=0\) and \(\frac{dS_2}{dt}=0\) and fixed points \(\frac{dS_1}{dt}=\frac{dS_2}{dt}=0\) for a given set of parameters.

Parameters:

• model (PyDSTool.Vode_ODEsystem) – Model object in PyDSTool.

• params (dict) – Model parameters.

• xdomain (dict) – Variable and lower/upper bounds.

• args (Argparse.Namespace) – Optional extra arguments.

Returns:

• model (PyDSTool.Vode_ODEsystem) – Updated model object with given parameters.

• fps (PyDSTool.Toolbox.phaseplane.fixedpoint_2D) – Fixed points of the system.

• (nulls_x, nulls_y) (tuple) – Tuple containing nullclines (arrays of paired xs and ys).

OCD_modeling.analysis.compute_trajectories(model, n, tdata=[0, 1000])

Compute n trajectories from model, each with different initial conditions.

Parameters:

• model (PyDSTool.Vode_ODEsystem) – PyDSTool model object.

• n (int) – Number of trajectories to compute.

• tdata (list) – Time interval of the saved trajectories.

OCD_modeling.analysis.compute_equilibrium_point_curve(model, fps, pdomain)

Find equilibrium point curve(s) of the system, starting from each fixed point (if exist).

Parameters:

• model (PyDSTool.Vode_ODEsystem) – Model object in PyDSTool.

• fps (PyDSTool.Toolbox.phaseplane.fixedpoint_2D) – Fixed points of the system.

• pdomain (dict) – Free variable (or order parameter) to perform the bifurcation analyis from.

Returns:

cont – PyDSTool Continuation Class object populated with equilibrium point curves.

Return type:

PyDSTool.ContClass

Parallel processing

Since we conducted the stability analysis while varying coupling parameters, we wrapped the stability analysis from the section above such that we can parallelize its execution with custom ranges of parameters.

OCD_modeling.analysis.run_stability_analysis(order_params, default_params, args)

Starts a pool of parallel processes to run the stability analysis.

Parameters:

• order_params (dict) – Fixed parameters, for which to analyse the system using discretized values, for example {'C_21': np.linspace(0.2,0.8,4)}.

• default_params (dict) – Default model’s parameters.

• args (Argparse.Namespace) – Structure of necessary options. For example, must include args.compute_epc = True to compute equilibrium point curves.

Returns:

• outputs (list of dict) – Stability analysis of each values (or combination of values) given in order_params.

• futures (list of concurrent.futures.Future) – (deprecated) if using futures.concurrent library for parallel process, return the list of Future objects (https://docs.python.org/3/library/concurrent.futures.html#concurrent.futures.Future).

OCD_modeling.analysis.launch_stability_analysis(order_params, default_params, out_queue, args)

Ghost process that launches the stability analysis for a set of defined order parameter, creating a child process with a set timeout per child process, such that the stbaility analysis does not hang waiting for the continuation to terminate if it does not converge.

Parameters:

• order_params (dict) – Fixed parameters, for which to analyse the system using discretized values, for example {'C_21': np.linspace(0.2,0.8,4)}.

• default_params (dict) – Default model’s parameters.

• out_queue (Queue) – Output queue on which to append the results.

• args (Argparse.Namespace) – Structure of necessary options. For example, must include args.compute_epc = True to compute equilibrium point curves.

Plotting

We visualized phase space topology via a wrapper of PyDSTool functions to plot fixed points and vector fields, and customed plotting of nullclines:

OCD_modeling.analysis.plot_phasespace(model, fps, nullclines, trajs, ax=None, args=None)

Plot vector field, nullclines and fixed points of a model previously set with parameters.

Parameters:

• model (PyDSTool.Vode_ODEsystem) – PyDSTool model object.

• fps (PyDSTool.Toolbox.phaseplane.fixedpoint_2D) – Fixed points of he system

• nullclines (list) – list of PyDSTool nullcline objects, containing arrays of xs and ys of nullclines.

• trajs (list of dict) – Simulated trajectories.

• ax (matplotlib.Axis) – Axis in which to plot the phasespace.

It can then be called multiple times to plot grids of phase spaces as in Supplementary Section I.

OCD_modeling.analysis.plot_phasespace_grid(outputs, order_params, args=None)

Plot a grid of phasespaces from stability analysis outputs.

Parameters:

• outputs (list) – Outputs from stability analysis.

• order_params (dict) – Fixed parameters, for which to analyse the system using discretized values, for example {'C_21': np.linspace(0.2,0.8,4)}.

• args (Argparse.Namespace) – Optional extra arguments in argprase Namespace, such as args.save_figs=True to save figure and args.plot_figs=True to plot figures.

Grid of phase space with fixed point, nullclines and trajectories for different values of couplings.

OCD_modeling.analysis.plot_bifurcation_row(outputs, order_params, rww=None, t_range=None, args=None)

Plot a row of bifurcation diagrams (ie. a 1 by n grid).

Parameters:

• outputs (list) – Outputs from stability analysis.

• order_params (dict) – Order parameters of the analysis in {‘param_name’: np.array} format where np.array is the list of order parameters param_name.

• rww (OCD_modeling.models.ReducedWongWangOU) – Model instance that ran.

• t_range (list) – [start, stop] timestamp values of the RWW model traces to plot in the diagrams.

• args (Argparse.Namespace) – Extra options.

Bifurcation diagrams of the model variables \(S1\), \(S2\) with free parameter \(C_{12}\) and order parameter \(C_{21}\). Lines denote equilibrium point curves (plain: stable, dash: unstable); LP: saddle-node points (i.e. LimitPoint), B: Boundary point (bound values of the analysis).

OCD_modeling.analysis.plot_timeseries_phasespace_bif(outputs, rww, df_eta_sigma, args)

Show \(S1\), \(S2\), \(C_{12}\) timeseries, \(S1-S2\) phase space with trajectories and \(C_{12} - S1|S2\) phase space. with bifurcation diagram.

Parameters:

• outputs (list) – Outputs from stability analysis.

• rww (OCD_modeling.models.ReducedWongWangOU) – Model instance that ran.

• df_eta_sigma (pandas.DataFrame) – Data from the simulations varying eta and sigma parameters.

• args (Argparse.Namespace) – Extra options.

Timeseries of the model variables \(S1\), \(S2\) and stochastic coupling \(C_{12}\), with corresponding projections in \(S1-S2\) state space, bifurcation diagram with free parameter \(C_{12}\), transition rates and frontostriatal functional connectivity (FC) according to \(\eta_{12}\) and \(\sigma_{12}\) parameters. Note that transition rates and FC were precomputed separately and loaded for visualization only.