OSCR

Serotonergic modulation of motor subspace dynamics drives a sleep-independent quiescent state.

Code ↔ Paper

15 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 15 matches · 2 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted
  1. [1] § Methods › Eigenvector-randomization surrogate ↔ analysis_from_mat.py, lines 545–633 · score 0.97 · random orthonormal basis, surrogate covariance matrix, eigenvalue spectrum, QR decomposition, Gaussian matrix, geometric structure
  2. [2] § Methods › Hyperbolic and Euclidean multidimensional scaling ↔ analysis_from_mat.py, lines 636–667 · score 0.78 · triangle inequality, proper metric, chord distance, hypersphere, vectors, correlated
  3. [3] § Results › DRN 5-HT neuron activation produces bout type–dependent, graded suppression of the motor subspace ↔ analysis_from_mat.py, lines 545–633 · score 0.77 · random orthonormal basis, eigenvalue spectrum, geometric structure, covariance matrix, randomized, symmetric
  4. [4] § Methods › Hyperbolic and Euclidean multidimensional scaling ↔ analysis_from_mat.py, lines 830–868 · score 0.76 · scikit learn, Euclidean MDS, precomputed, classical, metric, dissimilarity
  5. [5] § Methods › Subspace angle analysis ↔ code/figure4/dPCA_example.m, lines 457–482 · score 0.74 · sound related subspace, activation subspace, principal angle, motor related, orth, alignment
  6. [6] § Results › DRN 5-HT activation modulates motor circuits to suppress sound-evoked responses ↔ code/figure4/dPCA_example.m, lines 457–482 · score 0.71 · sound related subspace, geometric relationship, principal angle, motor related, alignment, activation
  7. [7] § Results › DRN 5-HT activation modulates motor circuits to suppress sound-evoked responses ↔ code/figure3/dpca_example.m, lines 358–373 · score 0.67 · geometric relationship, related subspace, principal angle, motor related, alignment, activation
  8. [8] § Results › DRN 5-HT neuron activation produces bout type–dependent, graded suppression of the motor subspace ↔ code/tools/computePearsonDistanceMatrix.m, the whole file · a weak match · score 0.66 · perfectly anti correlated, Pearson correlation, uncorrelated, distances, matrix, neuron
  9. [9] § Methods › Hyperbolic and Euclidean multidimensional scaling ↔ analysis_from_mat.py, lines 37–100 · score 0.64 · Shepard diagrams, embedded distances, original distances, diagonal, Euclidean, Hyperbolic
  10. [10] § Results › DRN 5-HT neuron activation produces bout type–dependent, graded suppression of the motor subspace ↔ analysis_from_mat.py, lines 37–100 · score 0.64 · original pairwise distances, Shepard diagram, Pairwise distance matrices, Original distance matrices, Euclidean, Embedding
  11. [11] § Methods › Subspace angle analysis ↔ code/figure3/dpca_example.m, lines 358–373 · score 0.62 · activation subspace, Principal angles, motor related, orth, alignment
  12. [12] § Methods › Subspace angle analysis ↔ code/figure4/subspace_angle_analysis.m, the whole file · a weak match · score 0.62 · random subspace, principal angle, orth, sound, tailed, DRN
  13. [13] § Methods › Hyperbolic and Euclidean multidimensional scaling ↔ metric_HMDS.py, lines 162–184 · score 0.59 · metric hmds, hyperbolic distances, Poincar, optimized, embedded
  14. [14] § Results › DRN 5-HT neuron activation produces bout type–dependent, graded suppression of the motor subspace ↔ analysis_from_mat.py, lines 724–827 · score 0.59 · hyperbolic fit, pairwise distances, Hyperbolic embeddings, baseline, space, curvature
  15. [15] § Methods › Hyperbolic and Euclidean multidimensional scaling ↔ analysis_from_mat.py, lines 830–868 · score 0.56 · Euclidean MDS, pairwise distance, distance matrix, uncertainty, multidimensional, variance

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 923 lines · 39 KB · MIT · 8 matches

  1. import metric_HMDS as HMDS
  2. import numpy as np
  3. from scipy.io import loadmat
  4. from scipy.io.matlab import MatlabOpaque
  5. import sys
  6. import argparse
  7. import contextlib
  8. import pickle
  9. import os
  10. import io
  13. # --- Code Quality & Clarity Suggestion: Use h5py for modern .mat files ---
  14. # MATLAB v7.3+ files are HDF5 format and require h5py to be read.
  15. # We import it here and will use it in the loading function.
  16. try:
  17. import h5py
  18. except ImportError:
  19. h5py = None
  20. # Add matplotlib and scikit-learn for plotting and MDS/PCA
  21. try:
  22. import matplotlib.pyplot as plt
  23. from sklearn.manifold import MDS
  24. from sklearn.decomposition import PCA
  25. from scipy.interpolate import interp1d
  26. from sklearn.metrics import pairwise_distances
  27. except ImportError:
  28. plt = None
  29. PCA = None
  30. MDS = None
  33. class MatReadError(Exception):
  34. """Custom exception for MATLAB file reading errors."""
  35. pass
  37. def plot_shepard_diagram(original_dmat, embedded_dmat, lambda_val=None, title="Shepard Diagram", output_file=None, sample_fraction=1.0):
  38. """
  39. Creates and displays a Shepard diagram to compare original and embedded distances.
  41. Args:
  42. original_dmat (np.ndarray): The original, normalized distance matrix.
  43. embedded_dmat (np.ndarray): The distance matrix from the embedding.
  44. lambda_val (float): The fitted curvature scale parameter.
  45. title (str): The title for the plot.
  46. output_file (str, optional): If provided, saves the plot to this file path as a PDF.
  47. sample_fraction (float, optional): Fraction of points to randomly sample for plotting (e.g., 0.1 for 10%). Default is 1.0 (all points).
  48. """
  49. if plt is None:
  50. print("\nWarning: matplotlib is not installed. Skipping plot. Please run 'pip install matplotlib'.")
  51. return
  53. print("Generating Shepard diagram...")
  55. # We use the upper triangle to avoid plotting each pair twice and the diagonal.
  56. N = original_dmat.shape[0]
  57. triu_indices_row, triu_indices_col = np.triu_indices(N, k=1)
  59. if sample_fraction < 1.0:
  60. num_pairs = len(triu_indices_row)
  61. sample_size = int(num_pairs * sample_fraction)
  62. if sample_size < 1:
  63. print(f"Warning: Sample fraction {sample_fraction} is too small for the number of pairs. No points will be plotted.")
  64. return
  65. print(f"Sampling {sample_size} of {num_pairs} pairs for Shepard diagram ({sample_fraction:.1%}).")
  66. sampled_idx = np.random.choice(num_pairs, size=sample_size, replace=False)
  67. indices = (triu_indices_row[sampled_idx], triu_indices_col[sampled_idx])
  68. else:
  69. indices = (triu_indices_row, triu_indices_col)
  71. original_distances = original_dmat[indices]
  72. if lambda_val is not None:
  73. embedded_distances = embedded_dmat[indices] / lambda_val
  74. ylabel = "Embedded Hyperbolic Distances / λ"
  75. else:
  76. embedded_distances = embedded_dmat[indices]
  77. ylabel = "Embedded Euclidean Distances"
  79. # Calculate R-squared (coefficient of determination)
  80. if len(original_distances) > 1 and len(embedded_distances) > 1:
  81. r_value = np.corrcoef(original_distances, embedded_distances)[0, 1]
  82. r_squared = r_value**2
  83. else:
  84. r_squared = np.nan
  86. fig, ax = plt.subplots(figsize=(8, 8))
  87. ax.scatter(original_distances, embedded_distances, alpha=0.5, s=15, edgecolors='k', linewidths=0.5)
  88. ax.plot([0, max(original_distances.max(), embedded_distances.max())], [0, max(original_distances.max(), embedded_distances.max())], 'r--', label=f'Perfect Match (y=x)\n$R^2 = {r_squared:.2f}$')
  89. ax.set_xlabel("Original Pairwise Distances (Normalized)")
  90. ax.set_ylabel(ylabel)
  91. ax.set_title(title)
  92. ax.grid(True)
  93. ax.set_aspect('equal', 'box')
  94. ax.legend()
  95. plt.show()
  97. if output_file:
  98. print(f"Saving Shepard diagram to {output_file}...")
  99. fig.savefig(output_file, format='pdf', bbox_inches='tight')
  100. plt.close(fig)
  102. def plot_poincare_2d(poincare_coords, title="2D Hyperbolic Embedding (Poincare Disk)", colors=None, output_file=None):
  103. """
  104. Creates a 2D visualization of points in the Poincare disk.
  106. Args:
  107. poincare_coords (np.ndarray): An N x 2 array of Poincare coordinates.
  108. title (str): The title for the plot.
  109. colors (optional): An array of colors for the points.
  110. output_file (str, optional): If provided, saves the plot to this file path.
  111. """
  112. if plt is None:
  113. print("\nWarning: matplotlib is not installed. Skipping 2D plot. Please run 'pip install matplotlib'.")
  114. return
  116. print("Generating 2D Poincare disk visualization...")
  117. fig, ax = plt.subplots(figsize=(8, 8))
  119. # Draw the boundary circle
  120. circle = plt.Circle((0, 0), 1.0, color='gray', fill=False, alpha=0.5)
  121. ax.add_artist(circle)
  123. # Use provided colors (and a colormap) or default to blue
  124. point_colors = colors if colors is not None else 'b'
  125. scatter = ax.scatter(poincare_coords[:, 0], poincare_coords[:, 1], c=point_colors, s=20, cmap='Reds')
  127. ax.set_title(title)
  128. ax.set_aspect('equal', 'box')
  130. # Add a colorbar if a color sequence is provided
  131. if colors is not None:
  132. fig.colorbar(scatter, ax=ax, label="Point Index")
  134. if output_file:
  135. print(f"Saving 2D plot to {output_file}...")
  136. fig.savefig(output_file, format='pdf', bbox_inches='tight')
  137. plt.close(fig)
  139. plt.show()
  141. def plot_poincare_3d(poincare_coords, colors=None, output_file=None):
  142. """Creates an interactive 3D visualization of points in the Poincare ball.
  144. Args:
  145. poincare_coords (np.ndarray): An N x 3 array of Poincare coordinates.
  146. colors (optional): An array of colors for the points.
  147. output_file (str, optional): If provided, saves the plot to this file path.
  148. """
  149. if plt is None:
  150. print("\nWarning: matplotlib is not installed. Skipping 3D plot. Please run 'pip install matplotlib'.")
  151. return
  153. if poincare_coords.shape[1] != 3:
  154. print(f"\nWarning: 3D plot is only available for 3D embeddings. Found {poincare_coords.shape[1]} dimensions. Skipping plot.")
  155. return
  157. print("Generating interactive 3D Poincare ball visualization...")
  159. fig = plt.figure(figsize=(9, 9))
  160. ax = fig.add_subplot(111, projection='3d')
  162. # Draw the boundary sphere (wireframe)
  163. u, v = np.mgrid[0:2*np.pi:20j, 0:np.pi:10j]
  164. x = np.cos(u)*np.sin(v)
  165. y = np.sin(u)*np.sin(v)
  166. z = np.cos(v)
  167. ax.plot_wireframe(x, y, z, color="gray", alpha=0.3)
  169. # Use provided colors (and a colormap) or default to blue
  170. point_colors = colors if colors is not None else 'b'
  171. scatter = ax.scatter(poincare_coords[:, 0], poincare_coords[:, 1], poincare_coords[:, 2], c=point_colors, s=20, cmap='Reds')
  173. # Add a colorbar if a color sequence is provided
  174. if colors is not None:
  175. fig.colorbar(scatter, ax=ax, shrink=0.6, aspect=20, label="Point Index")
  177. ax.set_xlabel("X coordinate")
  178. ax.set_ylabel("Y coordinate")
  179. ax.set_zlabel("Z coordinate")
  180. ax.set_title("3D Hyperbolic Embedding (Poincare Ball)")
  181. # Set aspect ratio to be equal
  182. ax.set_box_aspect([1,1,1])
  184. if output_file:
  185. print(f"Saving 3D plot to {output_file}...")
  186. fig.savefig(output_file, format='pdf', bbox_inches='tight')
  187. plt.close(fig)
  189. plt.show()
  191. def plot_poincare_3d_projections(poincare_coords, colors=None, output_file=None):
  192. """Creates 2D projections (XY, XZ, YZ) of a 3D Poincare embedding.
  194. Args:
  195. poincare_coords (np.ndarray): An N x 3 array of Poincare coordinates.
  196. colors (optional): An array of colors for the points.
  197. output_file (str, optional): If provided, saves the plot to this file path.
  198. """
  199. if plt is None:
  200. print("\nWarning: matplotlib is not installed. Skipping 2D projections plot.")
  201. return
  203. if poincare_coords.shape[1] != 3:
  204. return # Silently fail as the calling function should handle this.
  206. print("Generating 2D projections of the 3D embedding...")
  207. fig, axes = plt.subplots(1, 3, figsize=(24, 7))
  208. fig.suptitle("2D Projections of 3D Embedding", fontsize=16)
  210. projections = [
  211. (0, 1, 'X', 'Y'),
  212. (0, 2, 'X', 'Z'),
  213. (1, 2, 'Y', 'Z')
  214. ]
  216. for ax, (d1, d2, d1_name, d2_name) in zip(axes, projections):
  217. circle = plt.Circle((0, 0), 1.0, color='gray', fill=False, alpha=0.5)
  218. ax.add_artist(circle)
  219. scatter = ax.scatter(poincare_coords[:, d1], poincare_coords[:, d2], c=colors, s=20, cmap='Reds')
  220. ax.set_title(f"{d1_name} vs {d2_name} Projection")
  221. ax.set_xlabel(f"{d1_name} coordinate")
  222. ax.set_ylabel(f"{d2_name} coordinate")
  223. ax.set_aspect('equal', 'box')
  224. ax.grid(True)
  226. if colors is not None:
  227. fig.colorbar(scatter, ax=axes.ravel().tolist(), shrink=0.8, label="Point Index")
  229. if output_file:
  230. print(f"Saving 3D projections plot to {output_file}...")
  231. fig.savefig(output_file, format='pdf', bbox_inches='tight')
  232. plt.close(fig)
  234. plt.show()
  236. def visualize_embedding(fit_results, output_prefix=None):
  237. """
  238. Visualizes the embedding. For d=2 or d=3, it plots directly.
  239. For d > 3, it uses PCA to project the data to 2D for visualization.
  241. Args:
  242. fit_results (dict): The dictionary returned by run_embedding.
  243. output_prefix (str, optional): If provided, saves plots to files with this prefix.
  244. """
  245. if plt is None:
  246. print("\nWarning: matplotlib is not installed. Skipping visualization.")
  247. return
  249. coords = fit_results['cp']
  250. dim = coords.shape[1]
  252. # Generate a color array based on the index of each point
  253. num_points = coords.shape[0]
  254. colors = np.arange(num_points)
  256. def get_path(suffix):
  257. return f"{output_prefix}_{suffix}.pdf" if output_prefix else None
  259. if dim == 2:
  260. plot_poincare_2d(coords, colors=colors, output_file=get_path("2d"))
  261. elif dim == 3:
  262. plot_poincare_3d(coords, colors=colors, output_file=get_path("3d"))
  263. plot_poincare_3d_projections(coords, colors=colors, output_file=get_path("3d_projections"))
  264. elif dim > 3:
  265. print(f"\nEmbedding dimension is {dim} (>3). Projecting to 2D using PCA for visualization.")
  266. # Use PCA to find the three principal components
  267. pca = PCA(n_components=3)
  268. projected_coords = pca.fit_transform(coords)
  269. explained_variance = sum(pca.explained_variance_ratio_) * 100
  270. print(f"The 3D PCA projection explains {explained_variance:.2f}% of the variance.")
  272. title = f"PCA Projection of {dim}D Embedding to 3D Poincare Ball"
  273. plot_poincare_3d(projected_coords, colors=colors, output_file=get_path("pca_3d"))
  274. plot_poincare_3d_projections(projected_coords, colors=colors, output_file=get_path("pca_3d_projections"))
  276. def sample_submatrix(dmat, n, seed=None):
  277. """
  278. Randomly samples N neurons from an MxM distance matrix and returns an NxN submatrix.
  280. Args:
  281. dmat (np.ndarray): The MxM distance matrix.
  282. n (int): Number of neurons to sample.
  283. seed (int, optional): Random seed for reproducibility.
  285. Returns:
  286. tuple: (submatrix, indices) where submatrix is the NxN distance matrix
  287. and indices are the sampled row/column indices.
  288. """
  289. m = dmat.shape[0]
  290. if n > m:
  291. raise ValueError(f"Cannot sample {n} neurons from a {m}x{m} matrix.")
  292. rng = np.random.default_rng(seed)
  293. indices = np.sort(rng.choice(m, size=n, replace=False))
  294. submatrix = dmat[np.ix_(indices, indices)]
  295. print(f"Sampled {n} neurons from {m}x{m} matrix -> {n}x{n} submatrix. Indices: {indices}")
  296. return submatrix, indices
  299. def load_dmat_from_mat(mat_file_path, matrix_variable_name):
  300. """
  301. Loads a distance matrix from a .mat file, handling both old and v7.3 formats.
  303. Args:
  304. mat_file_path (str): Path to the .mat file.
  305. matrix_variable_name (str): The name of the distance matrix variable inside the .mat file.
  307. Returns:
  308. np.ndarray: The loaded distance matrix.
  309. """
  311. try:
  312. # First, try reading with scipy's loadmat, which handles older formats.
  313. mat_data = loadmat(mat_file_path, simplify_cells=True)
  314. except NotImplementedError:
  315. # This error is raised for v7.3 .mat files. We switch to h5py.
  316. print("Detected MATLAB v7.3 file, switching to h5py reader...")
  317. if h5py is None:
  318. raise MatReadError("MATLAB v7.3 file detected, but 'h5py' is not installed. Please run 'pip install h5py'.")
  320. try:
  321. with h5py.File(mat_file_path, 'r') as f:
  322. # h5py doesn't load the whole file into a dict, so we access the variable directly.
  323. if matrix_variable_name not in f:
  324. available_vars = list(f.keys())
  325. raise MatReadError(f"Variable '{matrix_variable_name}' not found in HDF5 file. Available variables: {available_vars}")
  327. # Extract the data and convert to a NumPy array.
  328. # Note: h5py may read data transposed compared to MATLAB. For a symmetric
  329. # distance matrix, this is not an issue.
  330. dmat = np.array(f[matrix_variable_name], dtype=float)
  331. except Exception as e:
  332. raise MatReadError(f"Failed to read HDF5 file with h5py: {e}")
  334. except FileNotFoundError:
  335. raise MatReadError(f"The file '{mat_file_path}' was not found.")
  337. except Exception as e:
  338. raise MatReadError(f"An unexpected error occurred while reading the .mat file: {e}")
  340. else:
  341. # This block runs if scipy.io.loadmat succeeded.
  342. if matrix_variable_name not in mat_data:
  343. available_vars = [k for k in mat_data.keys() if not k.startswith('__')]
  344. raise MatReadError(f"Variable '{matrix_variable_name}' not found in .mat file. Available variables: {available_vars}")
  346. dmat = mat_data[matrix_variable_name]
  348. # Handle complex data structures that can arise from MATLAB structs/cells
  349. if isinstance(dmat, MatlabOpaque):
  350. raise MatReadError("The variable is a complex MATLAB object. Please simplify it in MATLAB before saving.")
  352. # Ensure the matrix is a numpy array of the correct type (float)
  353. dmat = np.array(dmat, dtype=float)
  355. print(f"Successfully loaded matrix '{matrix_variable_name}' with shape: {dmat.shape}")
  357. return dmat
  359. def calculate_bic(dmat, emb_mat, n_params, is_hyperbolic=False, lambda_val=None, sig_vals=None, dmat_unc=None):
  360. """
  361. Calculates the Bayesian Information Criterion (BIC) for an embedding.
  363. Args:
  364. dmat (np.ndarray): The original distance matrix.
  365. dmat_unc (np.ndarray, optional): The matrix of uncertainties for each distance.
  366. emb_mat (np.ndarray): The distance matrix from the embedding.
  367. n_params (int): The number of parameters in the model.
  368. is_hyperbolic (bool): Flag to indicate if the model is hyperbolic.
  369. lambda_val (float, optional): The lambda scale parameter for hyperbolic models.
  370. sig_vals (np.ndarray, optional): The uncertainty parameters for hyperbolic models.
  372. Returns:
  373. float: The calculated BIC value.
  374. """
  375. N = dmat.shape[0]
  376. n_pairs = N * (N - 1) / 2
  377. indices = np.triu_indices(N, k=1)
  379. original_distances = dmat[indices]
  380. embedded_distances = emb_mat[indices]
  382. if is_hyperbolic:
  383. # Combine inferred and data uncertainties, matching the Stan model
  384. inferred_seff_sq = np.array([sig_vals[i]**2 + sig_vals[j]**2 for i in range(N) for j in range(i + 1, N)])
  385. data_unc_sq = dmat_unc[indices] if dmat_unc is not None else 0
  386. seff_sq = inferred_seff_sq + data_unc_sq
  388. residuals = original_distances - (embedded_distances / lambda_val)
  389. log_likelihood = -0.5 * np.sum(np.log(2 * np.pi * seff_sq) + (residuals**2 / seff_sq))
  390. else: # Euclidean
  391. # For the Euclidean case, we assume a single global model variance (sigma_model^2)
  392. # plus the known data variance for each pair.
  393. residuals = original_distances - embedded_distances
  394. rss = np.sum(residuals**2)
  396. # Estimate the single model variance parameter via MLE.
  397. # This is equivalent to the average residual sum of squares.
  398. sigma2_model_mle = rss / n_pairs
  400. data_unc_sq = dmat_unc[indices] if dmat_unc is not None else 0
  401. total_variance = sigma2_model_mle + data_unc_sq
  402. log_likelihood = -0.5 * np.sum(np.log(2 * np.pi * total_variance) + (residuals**2 / total_variance))
  404. return n_params * np.log(n_pairs) - 2 * log_likelihood
  406. def run_embedding(dmat, embedding_dim, dmat_unc=None, verbose=False, output_path=None):
  407. """
  408. Takes a distance matrix, normalizes it, and runs the HMDS embedding.
  410. Args:
  411. dmat (np.ndarray): The input distance matrix.
  412. dmat_unc (np.ndarray, optional): The matrix of uncertainties for each distance.
  413. embedding_dim (int): The target dimension for the hyperbolic embedding.
  414. verbose (bool): If True, prints Stan's optimization progress. Defaults to False.
  415. output_path (str, optional): If provided, saves the fit_results dictionary
  416. to this path using pickle.
  418. Returns:
  419. dict: The 'fit' dictionary containing all embedding results.
  420. """
  421. # Ensure the matrix is a numpy array of the correct type (float)
  422. if not isinstance(dmat, np.ndarray) or dmat.dtype != float:
  423. dmat = np.array(dmat, dtype=float)
  425. if dmat.ndim != 2 or dmat.shape[0] != dmat.shape[1]:
  426. raise ValueError(f"Input matrix must be square. Got shape: {dmat.shape}")
  428. # --- Pre-processing (as seen in tst.py) ---
  429. # It's often a good idea to normalize the distance matrix.
  430. # The original paper and tst.py normalize it to have a max value of 2.0.
  431. print("Normalizing distance matrix...")
  432. dmat = 2.0 * dmat / np.max(dmat)
  434. # --- Run the Embedding ---
  435. print(f"Starting HMDS embedding into {embedding_dim} dimensions...")
  436. # This step can take several minutes, especially the first time.
  437. if verbose:
  438. fit = HMDS.embed(embedding_dim, dmat, dij_unc=dmat_unc)
  439. else:
  440. # Suppress the C++ output from Stan by redirecting stdout and stderr
  441. print("Stan optimization output is suppressed. This may take a while...")
  442. f = io.StringIO()
  443. with contextlib.redirect_stdout(f), contextlib.redirect_stderr(f):
  444. fit = HMDS.embed(embedding_dim, dmat, dij_unc=dmat_unc)
  446. # --- Print Results ---
  447. print("\n--- Embedding Complete ---")
  448. print(f"Fitted curvature scale parameter (lambda): {fit['lambda']:.4f}")
  449. print(f"Embedding contains {len(fit['euc'])} points.")
  450. print("Results are available in the 'fit' dictionary.")
  452. # --- Calculate BIC ---
  453. n_params = dmat.shape[0] * embedding_dim + dmat.shape[0] + 1 - embedding_dim*(embedding_dim - 1)/2 # N*D + N + 1 - D*(D-1)/2
  454. fit['bic'] = calculate_bic(fit['dmat'], fit['emb_mat'], n_params,
  455. is_hyperbolic=True, lambda_val=fit['lambda'],
  456. sig_vals=fit['sig'], dmat_unc=dmat_unc)
  458. # --- Save Results to Disk ---
  459. if output_path:
  460. try:
  461. # Ensure the directory exists
  462. os.makedirs(os.path.dirname(output_path), exist_ok=True)
  463. with open(output_path, 'wb') as f_out:
  464. pickle.dump(fit, f_out)
  465. print(f"\nEmbedding results successfully saved to: {output_path}")
  466. except Exception as e:
  467. print(f"\nWarning: Could not save results to '{output_path}'. Error: {e}")
  469. return fit
  471. def run_embedding_trials(dmat, embedding_dim, n_trials=5, dmat_unc=None, verbose=False,
  472. output_path=None):
  473. """
  474. Runs the hyperbolic embedding multiple times and aggregates lambda across trials.
  476. The Stan optimizer uses a random initialization on each call, so lambda can
  477. vary across runs. This function collects all trial results, reports mean ± std,
  478. and returns the trial with the best (lowest) BIC as the representative fit.
  480. Args:
  481. dmat (np.ndarray): The input distance matrix.
  482. embedding_dim (int): Target embedding dimension.
  483. n_trials (int): Number of independent runs. Default 5.
  484. dmat_unc (np.ndarray, optional): Uncertainty matrix, same shape as dmat.
  485. verbose (bool): If True, prints Stan output for each trial.
  486. output_path (str, optional): If provided, saves the best-BIC fit to this path.
  488. Returns:
  489. dict with keys:
  490. 'best_fit' : full fit dict for the trial with the lowest BIC
  491. 'all_fits' : list of all n_trials fit dicts
  492. 'lambda_mean' : mean lambda across trials
  493. 'lambda_std' : std of lambda across trials
  494. 'lambda_all' : array of per-trial lambda values
  495. 'bic_all' : array of per-trial BIC values
  496. """
  497. lambda_vals = []
  498. bic_vals = []
  499. all_fits = []
  501. for i in range(n_trials):
  502. print(f"\n{'='*60}")
  503. print(f"Trial {i + 1} / {n_trials}")
  504. fit = run_embedding(dmat, embedding_dim, dmat_unc=dmat_unc, verbose=verbose)
  505. lambda_vals.append(fit['lambda'])
  506. bic_vals.append(fit['bic'])
  507. all_fits.append(fit)
  509. lambda_vals = np.array(lambda_vals)
  510. bic_vals = np.array(bic_vals)
  511. best_idx = int(np.argmin(bic_vals))
  512. best_fit = all_fits[best_idx]
  514. print(f"\n{'='*60}")
  515. print(f"Trial Summary ({n_trials} runs, dim={embedding_dim})")
  516. print(f"{'Trial':<8} {'lambda':>10} {'BIC':>14}")
  517. print("-" * 34)
  518. for i, (lam, bic) in enumerate(zip(lambda_vals, bic_vals)):
  519. marker = " <-- best BIC" if i == best_idx else ""
  520. print(f"{i + 1:<8} {lam:>10.4f} {bic:>14.2f}{marker}")
  521. print("-" * 34)
  522. print(f"{'Mean':<8} {lambda_vals.mean():>10.4f}")
  523. print(f"{'Std':<8} {lambda_vals.std():>10.4f} "
  524. f"({100 * lambda_vals.std() / lambda_vals.mean():.1f}% of mean)")
  526. if output_path:
  527. try:
  528. os.makedirs(os.path.dirname(output_path), exist_ok=True)
  529. with open(output_path, 'wb') as f_out:
  530. pickle.dump(best_fit, f_out)
  531. print(f"\nBest-BIC fit saved to: {output_path}")
  532. except Exception as e:
  533. print(f"\nWarning: Could not save results to '{output_path}'. Error: {e}")
  535. return {
  536. 'best_fit' : best_fit,
  537. 'all_fits' : all_fits,
  538. 'lambda_mean' : lambda_vals.mean(),
  539. 'lambda_std' : lambda_vals.std(),
  540. 'lambda_all' : lambda_vals,
  541. 'bic_all' : bic_vals,
  542. }
  545. def surrogate_distance_matrix(cov_mat, corr_unc=None, seed=None, distance_method='chord'):
  546. """
  547. Constructs a null-model distance matrix that preserves the eigenvalue spectrum
  548. of a covariance matrix but randomises its geometric structure.
  550. Procedure
  551. ---------
  552. 1. Eigendecompose the real covariance matrix: Cov = V @ diag(w) @ V.T
  553. 2. Draw a random orthonormal basis Q via QR decomposition of a Gaussian matrix.
  554. 3. Reconstruct a surrogate covariance matrix: Cov_surr = Q @ diag(w) @ Q.T
  555. 4. Normalise to a correlation matrix: Corr_surr[i,j] = Cov_surr[i,j] / sqrt(Cov_surr[i,i] * Cov_surr[j,j])
  556. 5. Convert to a distance matrix and, if corr_unc is provided, propagate
  557. the correlation uncertainty to a distance uncertainty via corr_unc_to_dist_unc.
  559. Args:
  560. cov_mat (np.ndarray): MxM real symmetric positive-semidefinite covariance matrix.
  561. corr_unc (np.ndarray, optional): MxM matrix of correlation uncertainties.
  562. The same uncertainty matrix is used for the surrogate (the surrogate
  563. randomises geometry, not measurement precision).
  564. seed (int, optional): Random seed for reproducibility.
  565. distance_method (str): 'chord' (default) or 'linear'.
  567. Returns:
  568. dict with keys:
  569. 'corr_surrogate' : MxM surrogate correlation matrix
  570. 'dmat_surrogate' : MxM surrogate distance matrix
  571. 'dmat_unc' : MxM propagated distance uncertainty (None if corr_unc not given)
  572. 'eigenvalues' : sorted eigenvalues (descending) from the real covariance matrix
  573. 'Q' : the random orthonormal basis used
  574. """
  575. rng = np.random.default_rng(seed)
  576. M = cov_mat.shape[0]
  578. # Step 1 — sanitize and symmetrise the input before eigendecomposition.
  579. if not np.all(np.isfinite(cov_mat)):
  580. raise ValueError("Input covariance matrix contains NaN or Inf values.")
  581. cov_in = cov_mat.astype(np.float64, copy=True)
  582. cov_sym = (cov_in + cov_in.T) / 2.0
  584. # Step 2 — eigendecompose; clip negative eigenvalues to 0 (numerical noise).
  585. w, _ = np.linalg.eigh(cov_sym) # ascending order
  586. n_neg_clipped = int(np.sum(w < 0))
  587. w = np.clip(w, 0.0, None)
  588. w = w[::-1].copy() # descending, contiguous for BLAS
  590. if n_neg_clipped > 0:
  591. print(f"Note: {n_neg_clipped} negative eigenvalue(s) clipped to 0 "
  592. f"(input matrix is not perfectly PSD).")
  593. print(f"Eigenvalue range: [{w[-1]:.4f}, {w[0]:.4f}], sum={w.sum():.4f}")
  595. # Step 3 — random orthonormal basis via QR of a Gaussian matrix.
  596. G = rng.standard_normal((M, M))
  597. Q, _ = np.linalg.qr(G)
  599. # Step 4 — surrogate covariance matrix.
  600. cov_surr = np.einsum('ik,k,jk->ij', Q, w, Q)
  602. # Step 5 — convert surrogate covariance to correlation matrix.
  603. std = np.sqrt(np.diag(cov_surr))
  604. # Guard against zero variance (degenerate directions).
  605. std = np.where(std == 0, 1.0, std)
  606. corr_surr = cov_surr / np.outer(std, std)
  608. # Numerical cleanup: clip to [-1, 1] and reset diagonal.
  609. corr_surr = np.clip(corr_surr, -1.0, 1.0)
  610. np.fill_diagonal(corr_surr, 1.0)
  611. if not np.all(np.isfinite(corr_surr)):
  612. n_bad = int(np.sum(~np.isfinite(corr_surr)))
  613. raise ValueError(f"Surrogate correlation matrix contains {n_bad} non-finite "
  614. f"entries after reconstruction. Check your input covariance matrix.")
  616. # Step 6 — correlation -> distance.
  617. dmat_surr = corr_to_distance(corr_surr, method=distance_method)
  619. # Step 7 — propagate uncertainty if provided.
  620. dmat_unc = None
  621. if corr_unc is not None:
  622. dmat_unc = corr_unc_to_dist_unc(corr_surr, corr_unc, method=distance_method)
  624. print(f"Surrogate built: eigenvalue range [{w[-1]:.4f}, {w[0]:.4f}], "
  625. f"distance range [{dmat_surr[dmat_surr > 0].min():.4f}, {dmat_surr.max():.4f}]")
  627. return {
  628. 'corr_surrogate': corr_surr,
  629. 'dmat_surrogate': dmat_surr,
  630. 'dmat_unc': dmat_unc,
  631. 'eigenvalues': w,
  632. 'Q': Q,
  633. }
  636. def corr_to_distance(corr_mat, method='chord'):
  637. """
  638. Converts a correlation matrix to a distance matrix.
  640. Two methods are available:
  642. 'chord' (default) — D_ij = sqrt(2 * (1 - C_ij))
  643. The chord distance between unit vectors on a hypersphere.
  644. Satisfies the triangle inequality; a proper metric.
  645. Recommended for HMDS because the embedding model assumes a valid metric space.
  647. 'linear' — D_ij = 1 - C_ij
  648. A simpler dissimilarity. Does NOT satisfy the triangle inequality in general
  649. (e.g. three vectors at 60°/60°/120° give D(u,w)=1.5 > D(u,v)+D(v,w)=1.0).
  650. Using this with HMDS is technically misspecified; it compresses mid-range
  651. correlations and can inflate the fitted lambda relative to 'chord'.
  653. Args:
  654. corr_mat (np.ndarray): MxM correlation matrix with values in [-1, 1].
  655. method (str): 'chord' (default) or 'linear'.
  657. Returns:
  658. np.ndarray: MxM distance matrix with zeros on the diagonal.
  659. """
  660. if method == 'chord':
  661. dmat = np.sqrt(np.clip(2.0 * (1.0 - corr_mat), 0.0, None))
  662. elif method == 'linear':
  663. dmat = 1.0 - corr_mat
  664. else:
  665. raise ValueError(f"Unknown method '{method}'. Choose 'chord' or 'linear'.")
  666. np.fill_diagonal(dmat, 0.0)
  667. return dmat
  670. def corr_unc_to_dist_unc(corr_mat, corr_var, method='chord', reg=1e-8):
  671. """
  672. Propagates per-element variance from a correlation matrix to a distance matrix
  673. using first-order error propagation (the delta method).
  675. The input and output are both **variances**, matching what the Stan model expects:
  676. seff = sqrt(sig[i]^2 + sig[j]^2 + deltaij_unc[i,j])
  677. where deltaij_unc is treated as variance (added to squared sigma terms).
  679. Because the linear distance D_linear = 1 - C, Var(C) = Var(D_linear), so you
  680. can pass the loaded variance matrix (e.g. var_dij_laserOn) directly as corr_var.
  682. Chord method — D_chord = sqrt(2*(1 - C))
  683. dD/dC = -1 / D_chord
  684. Var(D_chord) = (dD/dC)^2 * Var(C) = Var(C) / D_chord^2
  686. Diverges as D_chord → 0 (C → 1). The denominator is regularised:
  687. Var(D_chord) = Var(C) / max(D_chord, reg)^2
  689. Linear method — D_linear = 1 - C
  690. dD/dC = -1 → Var(D_linear) = Var(C) (pass-through, reg unused).
  692. Args:
  693. corr_mat (np.ndarray): MxM correlation matrix (values in [-1, 1]).
  694. corr_var (np.ndarray): MxM matrix of correlation variances (>= 0).
  695. Equal to Var(D_linear) since D_linear = 1 - C.
  696. method (str): 'chord' (default) or 'linear'.
  697. reg (float): Regularisation floor for the chord-distance denominator.
  698. Entries with D_chord < reg are treated as D_chord = reg.
  699. Default 1e-8 (negligible on the [0, 2] distance scale).
  701. Returns:
  702. np.ndarray: MxM distance-variance matrix. Diagonal is 0.
  703. Pass directly as dmat_unc to run_embedding().
  704. """
  705. if corr_var.shape != corr_mat.shape:
  706. raise ValueError("corr_var must have the same shape as corr_mat.")
  708. if method == 'chord':
  709. dmat = np.sqrt(np.clip(2.0 * (1.0 - corr_mat), 0.0, None))
  710. n_reg = np.sum((dmat < reg) & ~np.eye(dmat.shape[0], dtype=bool))
  711. if n_reg > 0:
  712. print(f"Note: {n_reg} off-diagonal pairs have D_chord < {reg:.0e}; "
  713. f"denominator floored at {reg:.0e} for variance propagation.")
  714. dmat_var = corr_var / np.maximum(dmat, reg) ** 2
  715. elif method == 'linear':
  716. dmat_var = corr_var.copy()
  717. else:
  718. raise ValueError(f"Unknown method '{method}'. Choose 'chord' or 'linear'.")
  720. np.fill_diagonal(dmat_var, 0.0)
  721. return dmat_var
  724. def outlier_sensitivity_analysis(dmat, embedding_dim, removal_fractions=(0.05, 0.10, 0.20),
  725. dmat_unc=None, verbose=False):
  726. """
  727. Assesses whether the hyperbolic fit is driven by a tail of outlier neurons.
  729. Neurons are ranked by their mean pairwise distance to all others (their
  730. "centrality" in distance space). The top removal_fractions of most-extreme
  731. neurons are removed one fraction at a time, and the hyperbolic embedding is
  732. refit on each pruned submatrix. The returned lambda values reveal whether
  733. curvature is a distributed population property or an artefact of a few outliers.
  735. Interpretation guide
  736. --------------------
  737. - lambda drops sharply (>30%) after removing 5-10%: curvature is outlier-driven;
  738. the hierarchical interpretation is weak.
  739. - lambda stays within ~20-30% of the full-data value through 20% removal:
  740. curvature reflects a distributed property of the population.
  742. Args:
  743. dmat (np.ndarray): The MxM distance matrix (raw, un-normalised).
  744. embedding_dim (int): Target dimension for each hyperbolic embedding.
  745. removal_fractions (tuple of float): Fractions of neurons to remove (cumulative).
  746. dmat_unc (np.ndarray, optional): Uncertainty matrix, same shape as dmat.
  747. verbose (bool): Whether to print Stan optimisation output.
  749. Returns:
  750. dict: Keys are 'full' and each fraction string (e.g. '5%').
  751. Values are dicts with 'lambda', 'n_neurons', 'removed_indices', 'fit'.
  752. """
  753. M = dmat.shape[0]
  755. # Rank neurons by mean pairwise distance (excluding self-distance on diagonal)
  756. mean_dist = (dmat.sum(axis=1) - np.diag(dmat)) / (M - 1)
  757. rank_order = np.argsort(mean_dist)[::-1] # descending: most extreme first
  759. results = {}
  761. # --- Full-data baseline ---
  762. print(f"\n{'='*60}")
  763. print(f"[Baseline] Fitting on all {M} neurons...")
  764. fit_full = run_embedding(dmat, embedding_dim, dmat_unc=dmat_unc, verbose=verbose)
  765. lambda_full = fit_full['lambda']
  766. results['full'] = {
  767. 'lambda': lambda_full,
  768. 'n_neurons': M,
  769. 'removed_indices': np.array([], dtype=int),
  770. 'fit': fit_full,
  771. }
  772. print(f"[Baseline] lambda = {lambda_full:.4f}")
  774. # --- Pruned subsets ---
  775. for frac in removal_fractions:
  776. n_remove = int(np.round(frac * M))
  777. removed_idx = rank_order[:n_remove]
  778. keep_idx = np.sort(np.setdiff1d(np.arange(M), removed_idx))
  779. n_keep = len(keep_idx)
  780. label = f"{int(frac * 100)}%"
  782. print(f"\n{'='*60}")
  783. print(f"[Remove top {label}] Removing {n_remove} most-extreme neurons, "
  784. f"refitting on {n_keep} neurons...")
  786. sub_dmat = dmat[np.ix_(keep_idx, keep_idx)]
  787. sub_unc = dmat_unc[np.ix_(keep_idx, keep_idx)] if dmat_unc is not None else None
  789. fit_pruned = run_embedding(sub_dmat, embedding_dim, dmat_unc=sub_unc, verbose=verbose)
  790. lambda_pruned = fit_pruned['lambda']
  791. pct_change = 100.0 * (lambda_pruned - lambda_full) / lambda_full
  793. results[label] = {
  794. 'lambda': lambda_pruned,
  795. 'n_neurons': n_keep,
  796. 'removed_indices': removed_idx,
  797. 'fit': fit_pruned,
  798. }
  799. print(f"[Remove top {label}] lambda = {lambda_pruned:.4f} "
  800. f"({pct_change:+.1f}% vs baseline)")
  802. # --- Summary table ---
  803. print(f"\n{'='*60}")
  804. print("Outlier Sensitivity Summary")
  805. print(f"{'Condition':<18} {'N neurons':>10} {'lambda':>10} {'Δλ vs full':>12}")
  806. print("-" * 52)
  807. print(f"{'Full data':<18} {M:>10} {lambda_full:>10.4f} {'—':>12}")
  808. for frac in removal_fractions:
  809. label = f"{int(frac * 100)}%"
  810. r = results[label]
  811. pct_change = 100.0 * (r['lambda'] - lambda_full) / lambda_full
  812. print(f"{'Remove top ' + label:<18} {r['n_neurons']:>10} "
  813. f"{r['lambda']:>10.4f} {pct_change:>+11.1f}%")
  815. stable = all(
  816. abs(results[f"{int(f*100)}%"]['lambda'] - lambda_full) / lambda_full <= 0.30
  817. for f in removal_fractions
  818. )
  819. print()
  820. if stable:
  821. print("Verdict: lambda is STABLE (<= 30% change). Curvature appears to be a "
  822. "distributed population property.")
  823. else:
  824. print("Verdict: lambda is UNSTABLE (> 30% change for at least one removal). "
  825. "The hyperbolic fit may be driven by a small tail of outlier neurons.")
  827. return results
  830. def run_euclidean_embedding(dmat, embedding_dim, dmat_unc=None):
  831. """
  832. Performs classical multidimensional scaling (MDS) for a Euclidean embedding.
  834. Args:
  835. dmat (np.ndarray): The input distance matrix.
  836. embedding_dim (int): The target dimension for the Euclidean embedding.
  837. dmat_unc (np.ndarray, optional): The matrix of uncertainties for each distance.
  839. Returns:
  840. dict or None: A dictionary containing the embedding results, or None if
  841. scikit-learn is not installed.
  842. """
  843. if MDS is None:
  844. print("\nWarning: scikit-learn is not installed. Skipping Euclidean MDS. Please run 'pip install scikit-learn'.")
  845. return None
  847. print(f"\nStarting Euclidean MDS embedding into {embedding_dim} dimensions...")
  849. # Normalize the distance matrix for fair comparison with hyperbolic results
  850. dmat = 2.0 * dmat / np.max(dmat)
  852. # scikit-learn's MDS expects dissimilarities, which is what dmat is.
  853. # We use metric=True for classical MDS.
  854. mds = MDS(n_components=embedding_dim, dissimilarity='precomputed', metric=True,
  855. random_state=0, normalized_stress=False)
  857. # Fit the model and get the embedded coordinates
  858. coords = mds.fit_transform(dmat)
  860. # Calculate the pairwise distances in the new Euclidean embedding
  861. embedded_dmat = pairwise_distances(coords, metric='euclidean')
  863. # --- Calculate BIC ---
  864. n_params = dmat.shape[0] * embedding_dim + 1 # N*D + 1 (for variance)
  865. bic = calculate_bic(dmat, embedded_dmat, n_params, dmat_unc=dmat_unc)
  868. return {'dmat': dmat, 'emb_mat': embedded_dmat, 'coords': coords, 'bic': bic}
  871. if __name__ == '__main__':
  872. # --- Code Quality & Clarity Suggestion: Using argparse ---
  873. # This makes the script more reusable by allowing you to pass file paths
  874. # and parameters from the command line instead of editing the code.
  875. parser = argparse.ArgumentParser(description="Run Hyperbolic MDS on a distance matrix from a .mat file.")
  876. parser.add_argument("mat_file", help="Path to the input .mat file.")
  877. parser.add_argument("matrix_name", help="Name of the distance matrix variable within the .mat file.")
  878. parser.add_argument("--unc-name", help="Name of the distance uncertainty matrix variable within the .mat file (optional).")
  879. parser.add_argument("-d", "--dim", type=int, default=3, help="Target dimension for the embedding (default: 3).")
  880. parser.add_argument("--no-plot", action="store_true", help="Suppress the Shepard diagram plot.")
  881. parser.add_argument("-o", "--output", help="Output file prefix to save plots as PDF files instead of displaying them.")
  882. parser.add_argument("--shepard-sample", type=float, default=1.0, help="Fraction of points to sample for the Shepard diagram (e.g., 0.1 for 10%). Default is 1.0 (all points).")
  884. args = parser.parse_args()
  886. try:
  887. # Command-line workflow
  888. print(f"Loading data from: {args.mat_file}")
  889. distance_matrix = load_dmat_from_mat(args.mat_file, args.matrix_name)
  890. distance_unc_matrix = None
  891. if args.unc_name:
  892. print(f"Loading uncertainty data from variable: {args.unc_name}")
  893. distance_unc_matrix = load_dmat_from_mat(args.mat_file, args.unc_name)
  895. # --- Hyperbolic Embedding ---
  896. fit_results = run_embedding(distance_matrix, args.dim,
  897. dmat_unc=distance_unc_matrix,
  898. verbose=True)
  899. print(f"Hyperbolic Embedding BIC: {fit_results['bic']:.2f}")
  901. if not args.no_plot:
  902. # Shepard diagram for Hyperbolic embedding
  903. plot_shepard_diagram(fit_results['dmat'], fit_results['emb_mat'],
  904. fit_results['lambda'], title="Shepard Diagram: Hyperbolic Embedding",
  905. output_file=f"{args.output}_hyperbolic_shepard.pdf" if args.output else None,
  906. sample_fraction=args.shepard_sample)
  907. # Visualize the embedding, using PCA if dimension > 3
  908. visualize_embedding(fit_results, output_prefix=args.output)
  910. # --- Euclidean Embedding ---
  911. euclidean_results = run_euclidean_embedding(distance_matrix, args.dim)
  912. print(f"Euclidean Embedding BIC: {euclidean_results['bic']:.2f}")
  914. if euclidean_results and not args.no_plot:
  915. # Shepard diagram for Euclidean embedding
  916. plot_shepard_diagram(euclidean_results['dmat'], euclidean_results['emb_mat'],
  917. title="Shepard Diagram: Euclidean Embedding",
  918. output_file=f"{args.output}_euclidean_shepard.pdf" if args.output else None,
  919. sample_fraction=args.shepard_sample)
  921. except MatReadError as e:
  922. print(f"Error: {e}", file=sys.stderr)
  923. sys.exit(1)

analysis_from_mat.py at commit e69dc9c, under MIT · at the source

Overview

Authors: Kexin Qi1,2, Yuming Chai1,2, Guodong Tan1,2, Daguang Li1,2, Quan Wen1,2,3
ORCID iDs: Yuming Chai, Quan Wen
  1. Division of Life Sciences and Medicine, University of Science and Technology of China, Hefei, China
  2. Hefei National Laboratory for Physical Sciences at the Microscale, Center for Integrative Imaging, University of Science and Technology of China, Hefei, China
  3. School of Data Science, University of Science and Technology of China, Hefei, China
Institutions: University of Science and Technology of China (China)
Journal: eLife, volume 15, article RP110370
Dates: published online 20 July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.7554/elife.110370 · PMID 42474031 · PMCID PMC13384494 · OpenAlex W7154474291
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: zebrafish (organism), systems (subfield)
Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Preprocessing, Machine learning, fMRI & imaging
Keywords: Zebrafish
MeSH: Dorsal Raphe Nucleus*, Motor Activity*, Serotonergic Neurons*, Serotonin*, Sleep*, Zebrafish*, Acoustic Stimulation, Animals, Larva, Optogenetics (* major topic)
Topic: Zebrafish Biomedical Research Applications (Cell Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: STI2030-Major Projects (2022ZD0211900)
Citations: not cited yet (Europe PMC); 73 references in the paper

Abstract

The dorsal raphe nucleus (DRN) serotonergic (5-HT) system has been implicated in regulating sleep and motor control; however, its specific role remains controversial. In this study, we found that optogenetic activation of DRN 5-HT neurons in larval zebrafish induced a quiescent state and a reduced response to acoustic stimuli. Unlike sleep, the induced quiescent state was not accompanied by a loss of postural control, and nighttime activation of DRN 5-HT neurons led to a subsequent sleep rebound. Whole brain light field imaging combined with demixed principal component analysis (dPCA) revealed distinct neural subspaces related to DRN activation, sound responses, and motor activity. DRN 5-HT activation selectively modulated the motor-related subspace while leaving the sound-evoked subspace unaffected. Unlike DRN activation, sleep induced by mepyramine significantly altered sound-evoked neuronal activity patterns. Further analysis demonstrated that serotonin had a graded effect on the motor subspace, wherein downstream neurons responsible for particular bout types were more significantly influenced. Embedding motor population activity in a curved geometric space revealed that the degree of curvature scales with behavioral suppression across animals, providing a quantitative signature of the quiescent state. Together, these results elucidate that serotonergic modulation promotes behavioral quiescence through selective regulation of motor populations.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 15 matches between paragraphs and lines of code.

wenquan/BayesianHMDS

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: e69dc9c35f04a35c0fafddd4a37bb14a29690007, 15 July 2026
Languages: Python (5), Jupyter (3), Shell (2)
Size: 28 files, 10 scripts
Software Heritage: not archived
Found in: “Data availability”
Holds: README, license file, environment (requirements.txt), 3 notebooks
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (8 files), Matplotlib (4 files), scikit-learn (2 files), SciPy (2 files), h5py (1 file), Stan (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
12 files
At the source: github.com/wenquan/BayesianHMDS

kexin2016/zebrafish-quiescent-state

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 05465637c1ffe7cab8a07ab1c5292799fd4431e0, 24 September 2026
Languages: MATLAB (47)
Size: 60 files, 47 scripts
Software Heritage: archived
Found in: the references
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Statistics and Machine Learning Toolbox (10 files), Image Processing Toolbox (5 files), Parallel Computing Toolbox (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
49 files
At the source: archive.softwareheritage.org/swh:1:dir:3bd9cda314fe354f3357a4418411457bcb080cc6;origin=https://github.com/kexin2016/zebrafish-quiescent-state;visit=swh:1:snp:b6ea94b82eb3be943e58206a2bb74ebe0e35a9a5;anchor=swh:1:rev:0211bfd1dc727c6222efd33971567e7976f7b75a

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 57 scripts, each with its path and the digest of its content;
  • 15 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data availability

The Bayesian multi-dimensional scaling (MDS) code used in this study is publicly available at https://github.com/wenquan/BayesianHMDS, copy archived at Wen, 2026. Code used for figure generation and data visualization is available at https://github.com/kexin2016/zebrafish-quiescent-state, copy archived at Qi, 2026. The behavioral and processed datasets generated and analyzed during this study have been deposited in Figshare and are publicly available at https://doi.org/10.6084/m9.figshare.32121937.

The following dataset was generated:

Qi K, Chai Y, Tan G, Li D, Wen Q. 2026. Serotonergic modulation of motor subspace dynamics drives a sleep-independent quiescent state. figshare.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 5 authors, 1 keyword, 10 MeSH terms, 1 funder, 69 references.

Cite

This paper

Qi, K., Chai, Y., Tan, G., Li, D., & Wen, Q. (2026). Serotonergic modulation of motor subspace dynamics drives a sleep-independent quiescent state. eLife, 15, RP110370. https://doi.org/10.7554/elife.110370

BibTeX

@article{qi2026serotonergic,
author = {Qi, Kexin and Chai, Yuming and Tan, Guodong and Li, Daguang and Wen, Quan},
title = {{Serotonergic modulation of motor subspace dynamics drives a sleep-independent quiescent state}},
journal = {eLife},
year = {2026},
month = jul,
volume = {15},
pages = {RP110370},
publisher = {eLife Sciences Publications, Ltd},
issn = {2050-084X},
doi = {10.7554/elife.110370},
url = {https://doi.org/10.7554/elife.110370},
pmid = {42474031},
pmcid = {PMC13384494}
}

RIS

TY - JOUR
AU - Qi, Kexin
AU - Chai, Yuming
AU - Tan, Guodong
AU - Li, Daguang
AU - Wen, Quan
TI - Serotonergic modulation of motor subspace dynamics drives a sleep-independent quiescent state
T2 - eLife
J2 - Elife
PY - 2026
DA - 2026/07/20
VL - 15
SP - RP110370
SN - 2050-084X
PB - eLife Sciences Publications, Ltd
DO - 10.7554/elife.110370
UR - https://doi.org/10.7554/elife.110370
LA - en
ER -

CSL-JSON

{
"id": "10.7554/elife.110370",
"type": "article-journal",
"title": "Serotonergic modulation of motor subspace dynamics drives a sleep-independent quiescent state",
"container-title": "eLife",
"author": [
{
"family": "Qi",
"given": "Kexin"
},
{
"family": "Chai",
"given": "Yuming"
},
{
"family": "Tan",
"given": "Guodong"
},
{
"family": "Li",
"given": "Daguang"
},
{
"family": "Wen",
"given": "Quan"
}
],
"container-title-short": "Elife",
"volume": "15",
"page": "RP110370",
"DOI": "10.7554/elife.110370",
"PMID": "42474031",
"PMCID": "PMC13384494",
"ISSN": "2050-084X",
"publisher": "eLife Sciences Publications, Ltd",
"URL": "https://doi.org/10.7554/elife.110370",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
20
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41467-026-75490-y [code]
Topographically organized dorsal raphe activity modulates forebrain sensory-motor representations and contributes to defensive behaviors.
Journal: Nature communications
In common: Parallel Computing Toolbox, Image Processing Toolbox, Statistics and Machine Learning Toolbox, zebrafish, 13 references
[2] doi:10.1038/s41467-026-72222-0 [code]
Eye movement kinematics reveal novel circadian organization of sleep substates.
Journal: Nature communications
In common: zebrafish, systems, 12 references
[3] doi:10.1016/j.neuron.2026.07.016 [code]
Inferring brain-wide interactions using data-constrained recurrent neural network models.
Journal: Neuron
In common: Parallel Computing Toolbox, Image Processing Toolbox, Statistics and Machine Learning Toolbox, 2 other tools, zebrafish, systems, 4 references
[4] doi:10.1038/s41467-026-76242-8 [code]
Whole-brain, all-optical interrogation of neuronal dynamics underlying gut and vascular interoception in zebrafish.
Journal: Nature communications
In common: Parallel Computing Toolbox, h5py, Image Processing Toolbox, 4 other tools, zebrafish, systems, 1 reference
[5] doi:10.1038/s41467-026-71725-0 [code]
Interactions across hemispheres in prefrontal cortex reflect global cognitive processing.
Journal: Nature communications
In common: Parallel Computing Toolbox, Image Processing Toolbox, Statistics and Machine Learning Toolbox, 4 other tools, 2 references
[6] doi:10.1016/j.neuron.2026.03.034 [code]
Dentate gyrus interneurons modulate winner-take-all network dynamics in freely behaving mice.
Journal: Neuron
In common: Parallel Computing Toolbox, h5py, Image Processing Toolbox, 5 other tools, systems
[7] doi:10.1038/s41467-026-73106-z [code]
Respiratory pauses highlight sleep architecture in mice.
Journal: Nature communications
In common: Parallel Computing Toolbox, h5py, Image Processing Toolbox, 4 other tools, 1 reference
[8] doi:10.1038/s41467-026-71568-9 [code]
Convergent and selective representations of pain, appetitive processes, aversive processes, and cognitive control in the insula.
Journal: Nature communications
In common: Parallel Computing Toolbox, Image Processing Toolbox, Statistics and Machine Learning Toolbox, 4 other tools, 1 reference
[9] doi:10.1038/s41592-026-03154-2 [code]
Simultaneous single-cell calcium imaging of neuronal population activity and brain-wide BOLD fMRI.
Journal: Nature methods
In common: Parallel Computing Toolbox, h5py, Image Processing Toolbox, 5 other tools
[10] doi:10.1038/s41586-026-10631-3 [code]
A prognostic human brain network for diffuse midline glioma.
Journal: Nature
In common: Parallel Computing Toolbox, h5py, Image Processing Toolbox, 5 other tools

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.