Model architecture and components design

A) Two-stage model architecture consisting of a static and a dynamic part. Inset: structure of the “drifting grating” visual stimulation trial. B) Schematics of spike trains recorded from neuron i and neuron j for three trials and corresponding output of the static part of the model Wst. C) Schematics of the dynamic part of the mode consisting of a Fourier features transformation, followed by recurrent long-short-term-memory model, multilayer perceptron, embedding, and positional encoding.

Model training and implementation

A) Loss optimization. Blue curve is actual spike rate for a representative unit taken from the VCN dataset (see below) as a function of trial time. The dotted yellow curve is a spike rate predicted by DyNetCP model during training. Inset: Loss function is designed to push the predicted spike rate to approach the actual spike rate during training. B) Model performance with and without Fourier features. Blue bars (left y-axis) represent a histogram of joint spiking probability (JPSTH anti-diagonal) for an example neuron pair with direct connection taken from the synthetic dataset (see below). Red curves (right y-axis) correspond to trial-averaged learned dynamic weights, while dotted magenta line represents static weight value. DyNetCP model without Fourier features (top plot) failed to properly localize the period where induced spiking is occurring, while use of Fourier features (bottom plot) enables to do so correctly.

Experimental dataset recorded in-vivo in 6 visual cortices.

Sessions used for model training and analysis taken from the publicly available Allen Brain Observatory-Visual Coding Neuropixels (VCN) database. We use a “functional connectivity” subset of these recordings corresponding to the “drifting grating 75 repeats” setting recorded for 26 animals. Filtering of recorded units (see statistical analysis criteria below) and classification of valid pairs results in a significant reduction of units (5.6% of raw units or 23% of clean units) used for dynamic connectivity inference. The last two columns represent the number and percentage of trials when the animal was running with average speed over 5cm/sec

Comparison of inferred static connectivity to alternative methods.

A) Comparison of static DyNetCP (red curve) with jitter-corrected CCG (black curve) for a representative pair of local units (session 771160300, spatially close units from Layer 5 of the VISam area). Residuals (bottom panel) and residual analysis (right panels). B) Comparison of static DyNetCP (red curve) with jitter-corrected CCG (black curve) for a representative pair of spatially distant units (session 771160300, units from VISam and VISp areas). C) Comparison of jitter-corrected CCG (black), GLMCC (magenta), and DyNetCP (red) weights for a pair of excitatory HH neurons from a synthetic dataset. D) MCC calculated for GLMCC (black) and for DyNetCP with αw =5SD (red) and 7SD (blue). E) Static connectivity Wst for a subset of excitatory HH neurons with GT (crossed circles) compared to connections recovered by GLMCC (red triangles) and DyNetCP (blue triangles).

TP and FP rates as a function of pair spike threshold

A) Recovery of True Positive (TP) connections in HH synthetic dataset. TP connections recovered by DyNetCP model at different pair spike thresholds αp. Yellow, green, and red curves correspond to peak threshold taken as 3SD, 5SD, and 7SD, respectively. TP connections recovered by GLMCC model are shown by blue curve. Since a predicted connection is counted as a TP only once (i.e., i → j and j→ i do not count as two negatives if i is not connected to j) and the number of GT connections decreases when αp increases, the MCC exhibits an increase beyond αp=100. B) Generation of False Positive (TP) connections in HH synthetic dataset. FP connections generated by the DyNetCP model at different pair spike rates. Yellow, green, and red curves correspond to peak threshold αp taken as 3SD, 5SD, and 7SD, respectively. FP connections recovered by GLMCC model are shown by blue curve.

Static connectivity inferred by the model is biologically interpretable.

A) Feedforward and feedback connections between 6 cortical areas ordered with respect to their anatomical hierarchy scores (after Ref.4). B) Histograms of the functional delays extracted from the jitter-corrected CCG for all pairs across all 26 sessions for each visual area combination. Dashed red line indicates the median of the distribution. C) Comparison of median functional delays vs anatomical hierarchy score difference extracted by DyNetCP (blue triangles) and by jitter-corrected CCG (red triangles) (VCN dataset, n=26 animals). D) Histograms of the functional delays inferred by static DyNetCP across all 26 sessions of VCN dataset for each visual area combination. Dashed red line indicates the median of the distribution.

Comparison of inferred dynamic connectivity to descriptive statistical methods.

A) An example of JPSTH for a local pair of excitatory units (session 766640955 of the VCN dataset, spatially close units from Layer 5 of the VISp area). Panels on the left and at the bottom show corresponding spike rates for receive and send units, correspondingly. White dotted square corresponds to the area of interest. See also (SM, Movie 1) B) Construction of the conditional cJPSTH. Pixels from the white dotted square in A) are normalized by spiking rate of the send neuron. Row of pixels marked by red rectangle is used in C). C) cJPSTH corresponding to the white dotted square in A). Red rectangle corresponds to red rectangle pixels in B). D) Dynamic weight Wdyn for the same pair of excitatory units inferred for L2 of 1 (see also SM, Movie 1). E) cJPSTH for the same pair of units.

Tuning loss parameter.

A) Set of dynamic weights Wdyn for a local pair of excitatory VISp neurons (session 771160300, same as in Fig.4D,E) learnt with different L2 values. The corresponding cJPSTH is shown on the bottom. Left column shows a set of corresponding static weights Wst. B) A histogram of Pearson correlation coefficients between the dynamic weights and cJPSTH for each pair (comparing each pixel of the weight/cJPSTH matrix) calculated for all 41 pairs found to be significant for this session. C) The two-sided paired equivalence tests (TOST) for all 41 pairs in this session indicate that the average difference between normalized weights inferred from both models is smaller than 0.09.

Experimental dataset used for reconstruction of local circuit in Fig.5 and Fig.6

Clean units (passed unit quality metrics) that from valid pairs (passed through all three thresholds αs, αp, and αw) from recordings in VISam area in session 819701982. Units are classified as excitatory (waveform duration longer than 0.4ms) or inhibitory (less than 0.4ms) and are assigned to cortical depth and cortical layers based on their spatial coordinates. Local circuit reconstructed from static weights matrices (max_static weight) is shown in Fig.5B. Maximum dynamic weight (max_dynamic_weight) is taken at functional time lag. Distance between units (send_recv_distance) is estimated from vertical positions of corresponding principal electrodes (probe_vertical_position).

Model enables discovery of cellular-level dynamic connectivity patterns in visual cortex.

A) Jitter-corrected static weights Wst for a connected pairs in an example local circuit (session # 819701982, VISam area). Units are numbered by local indices (Methods, Table 2). B) Dynamic weights Wdyn for the same pairs. C) Schematic diagram of static directional connectivity in pairs in a local circuit in A). Circles and triangles correspond to excitatory and inhibitory units, respectively. D) An example of a dynamic weight Wdyn for a pair (3→4) classified as un-connected.

Temporal variation of dynamic weight and reconstruction of dynamic local circuit.

A) Normalized dynamic weights for all valid pairs in an example local circuit (session # 819701982, VISam area). Units are numbered by local indices. B) Schematic diagram of dynamic connectivity in a local circuit for different trial times. Circles and triangles correspond to excitatory and inhibitory units, respectively. The size of blue (red) arrows reflects the dynamic weight of the inhibitory (excitatory) connection at a given time.

Reconstruction of a dynamic connectivity in a local circuit.

Dynamic weights Wdyn for 6 valid pairs in an example local circuit in Fig.6 taken at different trial times. A) at 20ms. B) at 70ms. C) at 150ms. D) at 200ms. Reconstruction of dynamic connectivity is shown on top of each figure. Circles and triangles correspond to excitatory and inhibitory units, respectively. The size of blue (red) arrows reflects the dynamic weight of the inhibitory (excitatory) connection at a given time. See also SM, Movie 2.

Dynamic functional connectivity prediction along the hierarchy axis of visual cortices.

A) Woff at different trial times (VISp area, n=26 animals). Units ordered by cortical depth and assigned to layers. The circles diameter is proportional to the dynamic offset. Red, blue, and magenta circles correspond to excitatory-excitatory, inhibitory-inhibitory, and mixed connections, respectively. B) Woff matrices for all visual areas (n=26 animals) ordered by their hierarchical

Intrinsic dimensionality of dynamic connectivity is large

A) Intrinsic dimensionality of neural activity manifold (session 831882777, no temporal smoothing) using spiking data for all VIS cortical units (red), units in selected pairs (blue), and for the same unit pairs using dynamic weights (black). Inset: Intrinsic dimensionality as a function of temporal gaussian kernel width. B) Intrinsic dimensionality (n=26 animals) in spike rate space normalized by the number of all VIS units (violet) and units in selected pairs (green). Intrinsic dimensionality in a weight space for units in selected pairs (orange). Box: SD, whiskers: 5%-95% range, box: mean, horizontal line: median, distribution: kernel Scott, **** is p-value<0.0001.

Selection of dynamic weight values and weights averaging

A)Selection of weights which contribute to averaging. Single trial spiking times for a given unit (black bars). Inferred dynamic weight values (orange curve) are selected at times when associated unit spikes (black dots), as the remaining weights do not contribute to the predicted spiking probability and therefore are not supervised to represent a meaningful value. B) Averaging of weights. Each blue point corresponds to a dynamic weight from a single trial. The orange curve represents a trial-average of every weight time-aligned with spikes. Green curve represents the Gaussian smoothing with a kernel including the 10 most recent bins.

Performance of the model on synthetic dataset with FP connections

A)Identification of FP connections for a “common input” triplet. Inset: schematics of a triplet with common input from the unit 0 to both units 1 and 2. Top row: connection is active during the first and last 125 bin period of the trial. Bottom row: Connection is active only during the middle 375 bin period. Blue bars (left y-axis) represent a histogram of joint spiking probability (JPSTH anti-diagonal). Red curves correspond to inferred dynamic weights (right y-axis), while dotted magenta line represents static weight value. Comparison of dynamic and static weights might help to distinguish TP connections 0→1 and 0→2 from an artifactual FP connection 1→2. B) Identification of FP connections for a “serial input” triplet. Inset: schematics of a triplet with serial input from the unit 0 to 1 followed by 1 to 2. Top row: connection is active during the first and last 125 bin period of the trial. Bottom row: Blue bars (left y-axis) represent a histogram of joint spiking probability (JPSTH anti-diagonal). Red curves correspond to inferred dynamic weights (right y-axis), while dotted magenta line represents static weight value. Comparison of dynamic and static weights may help to distinguish TP connections 0→1 and 1→2 from an artifactual FP connection 0→2.

Thresholds for units and pairs filtering

A) Number of valid units (session #771160300) that form pairs with distinct static peak (black curve) for different peak weight thresholds. Violet, green, blue, and red curves correspond to units that participate in formation of pairs, triplets, quintuplets, and connections above 10-tuples, correspondingly. B) Number of valid units that form pairs with static peak above 4SD (black curve) for different pair spike rate. Violet, green, blue, and red curves correspond to units that participate in formation of pairs, triplets, quintuplets, and connections above 10-tuples, correspondingly. C) 3D plot of a number of valid units that form pairs for different peak thresholds and pair spike rates.

Example of pairs and dynamic model performance

A) Example of a pair of units (session #771160300, excitatory units from VISam and VISpm areas with spike rates 1.92Hz and 0.35Hz) that pass a peak weight threshold of 7SD (left static weight matrix plot), however, has pair spike rate as small as 0.73 s-2. Corresponding cJPSTH (bottom right) and dynamic weight plots (top right) are sparse. B) Example of a pair of units (session #771160300, excitatory VISp unit with 7.2Hz spike rate and inhibitory VISp unit with 6.8Hz rate) that pass a peak weight threshold of 7SD (left static weight plot), however, has moderate pair spike rate of 50 s-2. Corresponding cJPSTH and dynamic weight plots are shown on the right. C) Example of a pair of units that pass peak weight threshold of 7SD (left static weight plot), that has a relatively high pair spike rate of 130 s-2. Corresponding cJPSTH and dynamic weight plots are shown on the right.

Filtered pairs and potential biases.

A) Histogram of pairs recorded in all VIS cortices for all 26 animals that pass 7SD peak weight threshold and are formed from units with spike rate exceeding 1Hz. Red arrow shows the pair spike threshold of 200 s-2. B) Histograms of inferred static and maximum dynamic weights for filtered pairs. C) Histograms of unit spike rates for clean units and for valid units belonging to filtered pairs. D) Histograms of waveform duration for clean units and valid units. Red arrow indicates a threshold to classify units as inhibitory or excitatory.