OSCR

PumpKin: A machine-learning pipeline for automatically tracking localized kinematics in freely moving C. elegans.

Code ↔ Paper

7 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 7 matches
  1. [1] § Materials and methods › Design of a custom Butterworth filter for high-frequency noise reduction ↔ library/pumpkin.py, lines 383–448 · score 0.82 · low pass Butterworth, signed magnitude, motion vectors, grinder points, signal, filter
  2. [2] § Results › Overview of the PumpKin pipeline ↔ library/pumpkin.py, lines 383–448 · score 0.74 · genetic strain, low pass filter, cutoff frequency, expected pumping rate, signal, food
  3. [3] § Results › Overview of the PumpKin pipeline ↔ library/pumpkin.py, lines 450–501 · score 0.74 · grinder motion signal, cutoff frequency, valid pumps, binning, smoothing, optional
  4. [4] § Materials and methods › Training the faster R-CNN networks ↔ library/training.py, lines 22–40 · score 0.71 · COCO v1, CNN model, FPN, torchvision, weights, Faster
  5. [5] § Materials and methods › Motion compensation for reduction of noise from body motion ↔ library/pumpkin.py, lines 187–234 · score 0.64 · affine transformation, matched feature, RANSAC, flow, motion, frame
  6. [6] § Materials and methods › Use of peristimulus time histogram (PSTH) to estimate a continuous pumping rate from event times ↔ library/pumpkin.py, lines 450–501 · score 0.63 · trough occurs, detected pump, histogram, binned, pumping rate, filtered
  7. [7] § Materials and methods › Training the faster R-CNN networks ↔ cropROI.ipynb, lines 18–79 · score 0.53 · ROI centered, pharyngeal bulb, cropped, FRCNN, track, videos

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 · 501 lines · 22 KB · MIT · 5 matches

  1. ################################################################################
  2. # pumpkin.py
  3. # Written by Erin Shappell for Lu Lab
  4. #
  5. # This module contains all functions used to generate a continuous pumping rate
  6. # estimate from an FRCNN-tracked grinder in freely moving C. elegans.
  7. #
  8. ################################################################################
  9. # Imports
  10. import re
  11. import os
  12. import cv2
  13. import numpy as np
  14. from scipy.signal import iirfilter, find_peaks, filtfilt
  15. from scipy.ndimage import gaussian_filter1d
  17. ################################################################################
  18. def sample_dark_pts(frame, point, box_size=(15, 15), grid_size=(3, 3)):
  19. """
  20. Sample the darkest pixels from a grid of subregions within a box centered at a given point.
  22. Inputs:
  23. frame (ndarray): Input video frame (color image as a NumPy array).
  24. point (tuple): (x, y) coordinates representing the center of the sampling box.
  25. box_size (tuple): (width, height) of the box around the center within which to sample.
  26. Default is (15, 15).
  27. grid_size (tuple): (rows, cols) specifying how the box is divided into subregions.
  28. Default is (3, 3).
  30. Output:
  31. list: A list of (x, y) tuples corresponding to the darkest pixel in each subregion,
  32. adjusted to the coordinates of the original frame.
  33. Returns an empty list if the region is invalid.
  34. """
  35. cx, cy = map(int, point) # Ensure integer coordinates
  36. w, h = box_size # Width and height of the bounding box
  38. # Compute bounding box coordinates
  39. x1, y1 = max(cx - w // 2, 0), max(cy - h // 2, 0)
  40. x2, y2 = min(cx + w // 2, frame.shape[1]), min(cy + h // 2, frame.shape[0])
  42. # Crop the region of interest (ROI)
  43. roi = frame[y1:y2, x1:x2]
  45. # Convert to grayscale if the image is color--if invalid, return empty list
  46. if roi is not None: roi = cv2.cvtColor(roi, cv2.COLOR_BGR2GRAY)
  47. else: return []
  49. # Grid division
  50. rows, cols = grid_size
  51. step_y, step_x = (y2 - y1) // rows, (x2 - x1) // cols # Compute step size
  52. sampled_points = []
  54. # Loop over grid cells
  55. for i in range(rows):
  56. for j in range(cols):
  57. # Define cell boundaries
  58. y_start, y_end = y1 + i * step_y, y1 + (i + 1) * step_y
  59. x_start, x_end = x1 + j * step_x, x1 + (j + 1) * step_x
  61. # Extract the sub-region
  62. cell = roi[(i * step_y):(i + 1) * step_y, (j * step_x):(j + 1) * step_x]
  64. # Ensure cell has valid size before processing
  65. if cell.size == 0: continue
  67. # Find the darkest point in this cell
  68. min_val, _, min_loc, _ = cv2.minMaxLoc(cell)
  70. # Adjust coordinates to global frame
  71. dark_x = x_start + min_loc[0]
  72. dark_y = y_start + min_loc[1]
  74. sampled_points.append((dark_x, dark_y))
  76. return sampled_points
  78. ################################################################################
  79. def calc_motion_vecs(dark_points_t1, dark_points_t2, threshold=200):
  80. """
  81. Calculate motion vectors by matching dark points between two consecutive frames.
  83. Inputs:
  84. dark_points_t1 (list): List of (x, y) coordinates from frame t-1.
  85. dark_points_t2 (list): List of (x, y) coordinates from frame t.
  86. threshold (float): Maximum distance allowed to consider two points as a match. Default is 200.
  88. Output:
  89. ndarray: Array of motion vectors [(dx, dy), ...] representing the displacement
  90. of matched points from t-1 to t. Only vectors below the threshold are included.
  91. """
  92. # Initialize array of motion vectors
  93. motion_vectors = []
  95. for pt1 in dark_points_t1:
  96. # Calculate distances from pt1 to all points in t2
  97. distances = np.linalg.norm(np.array(dark_points_t2) - np.array(pt1), axis=1)
  99. # Find the closest point in t2
  100. min_distance = np.min(distances)
  101. if min_distance <= threshold:
  102. # Get the index of the closest point
  103. closest_idx = np.argmin(distances)
  104. pt2 = dark_points_t2[closest_idx]
  106. # Calculate the motion (difference in coordinates)
  107. dx = pt2[0] - pt1[0]
  108. dy = pt2[1] - pt1[1]
  109. motion_vectors.append((dx, dy))
  111. return np.array(motion_vectors)
  113. ################################################################################
  114. def find_troughs_and_peaks(signal, distance=10, height=[0.1,5]):
  115. """
  116. Detect peaks (local maxima) that are followed by troughs (local minima) within a specified distance.
  118. Inputs:
  119. signal (ndarray): 1D array of signal values.
  120. distance (int): Maximum number of frames allowed between a peak and its corresponding trough.
  121. Default is 10.
  122. height (list): Two-element list specifying the minimum and maximum height of peaks/troughs.
  123. Default is [0.1, 5].
  125. Outputs:
  126. tuple:
  127. - troughs (ndarray): Indices of troughs that follow valid peaks within the specified distance.
  128. - peaks (ndarray): Indices of peaks that are followed by a trough within the specified distance.
  129. """
  130. # Find troughs (invert signal to detect minima)
  131. troughs, _ = find_peaks(-signal, distance=distance/2, height=(height[0], height[1]))
  133. # Find peaks
  134. peaks, _ = find_peaks(signal, distance=distance/2, height=(height[0], height[1]))
  136. # Keep only peaks that come after a trough
  137. valid_peaks = []
  138. valid_troughs = []
  140. for peak in peaks:
  141. # Find the first trough after the peak
  142. following_troughs = troughs[troughs > peak]
  144. # Remove any already matched troughs to ensure exclusivity
  145. following_troughs = np.array([t for t in following_troughs if t not in valid_troughs])
  147. if len(following_troughs) > 0 and np.abs(following_troughs[0] - peak) < distance:
  148. valid_troughs.append(following_troughs[0])
  149. valid_peaks.append(peak)
  151. return np.array(valid_troughs), np.array(valid_peaks)
  153. ################################################################################
  154. def make_bin_edges(bin_width, dt, total_time):
  155. """
  156. Generate timestamp edges to segment a time series into bins of fixed width.
  158. Inputs:
  159. bin_width (float): Desired width of each bin in seconds.
  160. dt (float): Time step between consecutive samples in seconds.
  161. total_time (float): Total duration of the signal in seconds.
  163. Outputs:
  164. ndarray: Array of timestamps representing the edges of each bin, starting at 0 and
  165. covering the total_time in increments of bin_width.
  166. """
  167. # Compute the timesteps (same code as previous)
  168. ts = np.arange(0, total_time, dt)
  170. # Num timesteps
  171. n_timesteps = np.prod(ts.shape)
  173. # Warn if binsize doesn't divide the timestep evenly
  174. if (bin_width % dt) >= dt :
  175. print("Warning: bin_width doesn't evenly divide the timestep")
  176. print(bin_width % dt)
  178. # How many timesteps per bin?
  179. steps_per_bin = np.around(bin_width / dt)
  181. # Compute the bin edges using np.arange
  182. # remember that arange is an open interval at the end
  183. # so add one to the number of timesteps
  184. bin_edges = np.arange(0, n_timesteps+1, steps_per_bin) * dt
  185. return bin_edges
  187. ################################################################################
  188. def motion_comp(prev_frame, curr_frame, num_points=500, points_to_use=500):
  189. """
  190. Estimate the motion transformation matrix between two sequential image frames.
  191. Contains code adapted from [https://github.com/itberrios/CV_projects/tree/main]
  193. Inputs:
  194. prev_frame (ndarray): The first image frame (RGB).
  195. curr_frame (ndarray): The second sequential image frame (RGB).
  196. num_points (int): Number of feature points to detect in the first frame.
  197. points_to_use (int): Number of matched points to use for motion estimation.
  199. Outputs:
  200. A (ndarray or None): Estimated affine transformation matrix (2x3) or None if estimation fails.
  201. prev_points (ndarray): Feature points detected in the previous frame.
  202. curr_points (ndarray): Corresponding matched points in the current frame.
  203. """
  204. # Convert to grayscale
  205. prev_gray = cv2.cvtColor(prev_frame, cv2.COLOR_RGB2GRAY)
  206. curr_gray = cv2.cvtColor(curr_frame, cv2.COLOR_RGB2GRAY)
  208. # Get features for first frame
  209. features = cv2.goodFeaturesToTrack(prev_gray, num_points, qualityLevel=0.01, minDistance=10)
  211. # If no feature points exist, exit
  212. if features is None: return None, [], []
  214. # Get matching features in next frame with Sparse Optical Flow Estimation
  215. matched_features, status, _ = cv2.calcOpticalFlowPyrLK(prev_gray, curr_gray, features, None)
  217. # Reformat previous and current feature points
  218. prev_points = features[status==1]
  219. curr_points = matched_features[status==1]
  221. # Subsample number of points so we don't overfit
  222. if points_to_use > prev_points.shape[0]: points_to_use = prev_points.shape[0]
  224. index = np.random.choice(prev_points.shape[0], size=points_to_use, replace=False)
  225. prev_points_used = prev_points[index]
  226. curr_points_used = curr_points[index]
  228. # If no feature points exist, exit
  229. if len(prev_points_used) == 0 or len(curr_points_used) == 0: return None, [], []
  231. # Find transformation matrix from frame 1 to frame 2
  232. A, _ = cv2.estimateAffine2D(prev_points_used, curr_points_used, method=cv2.RANSAC)
  234. return A, prev_points, curr_points
  236. ################################################################################
  237. def get_grinder_motion(prev_frame, curr_frame, prev_com, curr_com):
  238. """
  239. Detects matched grinder keypoints and computes motion vectors between two frames.
  240. Contains code adapted from [https://github.com/itberrios/CV_projects/tree/main]
  242. Inputs:
  243. prev_frame (ndarray): Previous video frame.
  244. curr_frame (ndarray): Current video frame.
  245. prev_com (tuple): Centroid (x, y) of the grinder in the previous frame.
  246. curr_com (tuple): Centroid (x, y) of the grinder in the current frame.
  248. Outputs:
  249. prev_grinder_pts (ndarray): Transformed grinder points from previous frame.
  250. curr_grinder_pts (ndarray): Grinder points sampled in current frame.
  251. motion (ndarray): Motion vectors between matched points.
  252. magnitude (ndarray): Magnitude of motion vectors.
  253. angle (ndarray): Angles (radians) of motion vectors.
  254. """
  256. ### Get affine transformation for frame alignment
  257. # Get frame info
  258. h, w, _ = prev_frame.shape
  260. # Get affine transformation matrix for motion compensation between frames
  261. A, prev_pts, curr_pts = motion_comp(prev_frame, curr_frame, num_points=10000, points_to_use=5000)
  263. ### Transform previous frame's points using affine transformation
  264. # First, check that A was obtained correctly
  265. if A is None: return [],[],[],[],[]
  267. # Get transformed grinder points from the previous frame using A
  268. A = np.vstack((A, np.zeros((3,)))) # get 3x3 matrix to transform points
  269. prev_grinder_pts = sample_dark_pts(prev_frame, prev_com)
  270. curr_grinder_pts = sample_dark_pts(curr_frame, curr_com)
  272. # Check if points aren't found
  273. if not prev_grinder_pts or not curr_grinder_pts: return [],[],[],[],[]
  275. # Compensate the previous frame's points for motion
  276. comp_pts = np.hstack((prev_grinder_pts, np.ones((len(prev_grinder_pts), 1)))) @ A.T
  277. comp_pts = comp_pts[:, :2]
  279. ### Obtain motion information about the grinder points
  280. motion = calc_motion_vecs(comp_pts, curr_grinder_pts)
  281. if len(motion) < 2: return [],[],[],[],[]
  283. magnitude = np.linalg.norm(motion, ord=2, axis=1)
  284. angle = np.arctan2(motion[:, 0], motion[:, 1])
  286. return comp_pts, curr_grinder_pts, motion, magnitude, angle
  288. ################################################################################
  289. def process_video(video_path=None, grinder_coms=None):
  290. """
  291. Processes a video to compute motion vectors between grinder positions frame-by-frame.
  293. Inputs:
  294. video_path (str): Path to the input video file.
  295. grinder_coms (ndarray): Array of grinder center-of-mass positions per frame, shape (n_frames, 2).
  297. Outputs:
  298. mags (list of ndarray): List of arrays containing motion vector magnitudes per frame.
  299. angs (list of ndarray): List of arrays containing motion vector angles (radians) per frame.
  300. """
  301. ### Check that a video path and grinder CoMs were both provided
  302. if video_path is None: raise ValueError("No video path was provided.")
  303. if grinder_coms is None: raise ValueError("No grinder CoMs were provided.")
  305. ### Load video
  306. video_name = os.path.splitext(os.path.basename(video_path))[0]
  307. vid = cv2.VideoCapture(video_path)
  308. width = int(vid.get(cv2.CAP_PROP_FRAME_WIDTH))
  309. height = int(vid.get(cv2.CAP_PROP_FRAME_HEIGHT))
  310. fps = vid.get(cv2.CAP_PROP_FPS)
  311. tot_frames = int(vid.get(cv2.CAP_PROP_FRAME_COUNT))
  313. ### Initialize lists for motion information and labeled frames
  314. mags, angs = [], []
  315. saved_pts, frames = [], []
  317. # Initialize frames and coms
  318. success, prev_frame = vid.read()
  319. prev_com = grinder_coms[0]
  320. i = 1
  321. while success and i < grinder_coms.shape[0]:
  322. success, curr_frame = vid.read()
  323. curr_com = grinder_coms[i] if success else None
  324. curr_pts = []
  326. if success:
  327. # If a grinder was detected, calculate the motion
  328. if any(prev_com) and any(curr_com):
  329. # Get the distances between the grinder points in the previous and current frames
  330. prev_pts, curr_pts, motion, mag, ang = get_grinder_motion(prev_frame, curr_frame, prev_com, curr_com)
  332. # Save the motion information from the current frame
  333. mags.append(mag)
  334. angs.append(ang)
  336. # Draw detected motion vectors
  337. if curr_pts and prev_pts.shape == motion.shape and save:
  338. plot_frame = plot_vecs(prev_frame.copy(), np.hstack([prev_pts,prev_pts+motion]))
  340. # If a grinder was NOT detected, fill with blanks and continue
  341. else:
  342. mags.append([])
  343. angs.append([])
  345. # Save previous frame and grinder CoM for next iteration
  346. prev_frame = curr_frame.copy()
  347. prev_com = curr_com
  349. # Update iterator
  350. i += 1
  352. # Save the cluster locations from each frame
  353. saved_pts.append(curr_pts)
  355. return mags, angs
  357. ################################################################################
  358. def get_strain_cond(video_name):
  359. """
  360. Parses a video filename to extract the strain and condition identifiers.
  362. Inputs:
  363. video_name (str): Full name of the video file (e.g., 'n2_ff_00001.wmv' or 'n2_ff_00001_cropped.wmv').
  365. Outputs:
  366. strain (str): Extracted strain identifier from the filename (e.g., 'n2').
  367. condition (str): Extracted condition identifier from the filename (e.g., 'ff').
  369. Raises:
  370. ValueError: If the filename does not match the expected format 'STRAIN_CONDITION_#####.ext'.
  371. """
  372. base = os.path.splitext(video_name)[0] # Remove extension, e.g., '.wmv'
  374. # Match STRAIN_CONDITION_##### optionally followed by '_cropped'
  375. match = re.match(r'^([^_]+)_([^_]+)_\d+(?:_cropped)?$', base)
  377. if not match:
  378. raise ValueError("Input string must be in format 'STRAIN_CONDITION_#####.ext'")
  380. strain, condition = match.group(1), match.group(2)
  381. return strain, condition
  383. ################################################################################
  384. def get_filtered_motion(strain_name=None, cond_name=None, mags=None, angs=None, fps=20):
  385. """
  386. Calculates the signed and filtered motion magnitude of grinder points across frames.
  388. Inputs:
  389. strain_name (str): Name of the genetic strain.
  390. cond_name (str): Name of the experimental condition.
  391. mags (list of ndarray): List of motion magnitude arrays per frame.
  392. angs (list of ndarray): List of motion angle arrays per frame (radians).
  393. fps (float): Sampling frequency of the video frames (default is 20).
  395. Outputs:
  396. signed_mags_filt (ndarray): Low-pass filtered signed average motion magnitude signal.
  397. fc (float): Cutoff frequency used for filtering based on the strain and condition.
  398. """
  399. ### Check that a strain name, condition name, and motion information were all provided
  400. if strain_name is None: raise ValueError("No strain name was provided.")
  401. if cond_name is None: raise ValueError("No condition name was provided.")
  402. if mags is None: raise ValueError("No magnitudes were provided.")
  403. if angs is None: raise ValueError("No angles were provided.")
  405. ### Determine if any of the grinder points have moved significantly
  406. # Store magnitudes
  407. max_mag = 5 # ignore motion vectors exceeding this threshold
  408. n_frames = len(mags)
  409. avg_mags = np.zeros(n_frames) # unsigned magnitude of motion
  410. avg_angs = np.zeros(n_frames) # average angle of motion
  411. signs = np.zeros(n_frames) # negative/positive motion
  412. avg_signed_mags = np.zeros(n_frames) # signed magnitude of motion
  414. # Main loop for calculating signed motion magnitude
  415. for i in range(n_frames):
  416. # Save the average angle values (used to determine sign)
  417. if len(angs[i]) > 0: avg_angs[i] = np.mean(angs[i])
  418. else: avg_angs[i] = 0
  420. # Determine the sign based on the angle value
  421. if avg_angs[i] > 0: signs[i] = 1
  422. elif avg_angs[i] < 0: signs[i] = -1
  424. # Save the largest magnitude + average magnitude from each frame
  425. if len(mags[i]) > 0 and len([x for x in mags[i] if x < max_mag]) > 0:
  426. avg_mags[i] = np.mean([x for x in mags[i] if x < max_mag])
  427. else: avg_mags[i] = 0
  429. # Save the signed magnitude
  430. avg_signed_mags[i] = avg_mags[i]*signs[i]
  432. ### Build low-pass filter and apply to signal
  433. ### NOTE: if using a new strain, you will need to add it as follows:
  434. ### elif strain_name == 'new_strain' and (cond_name == 'ff' or cond_name == 'sf'):
  435. ### fc = 2*expected_max_rate_for_new_strain
  436. # Use a cutoff frequency (fc) equal to 2*max_expected_pumping_rate
  437. if strain_name == 'n2' and (cond_name == 'ff' or cond_name == 'sf'):
  438. fc = 10 # N2 on food has a higher expected pumping rate
  439. elif strain_name == 'eat2' and (cond_name == 'ff' or cond_name == 'sf'):
  440. fc = 4 # eat-2 on food have lower expected pumping rate
  441. elif cond_name == 'fs' or cond_name == 'ss':
  442. fc = 2 # worms off food have signficantly lower expected pumping rate
  444. Wn = fc / (fps/2) # fps is the sampling frequency
  445. b, a = iirfilter(N=2, Wn=Wn, btype='low', ftype='butter') # low-pass Butterworth filter
  446. signed_mags_filt = filtfilt(b, a, avg_signed_mags) # filtered signed magnitude
  448. return signed_mags_filt, fc
  450. ################################################################################
  451. def get_pumping_rate(video_path=None, grinder_motion=None, fps=20, fc=None,
  452. min_height=0.5, max_height=5, sigma=1, save=False, save_path=None):
  453. """
  454. Estimates the continuous and discrete pumping rate of a grinder based on motion signals.
  456. Inputs:
  457. video_path (str): Path to the input video file.
  458. grinder_motion (ndarray): Array of grinder motion magnitudes per frame.
  459. fps (float): Frame rate of the video.
  460. fc (float): Cutoff frequency used for filtering based on the strain.
  461. min_height (float): Minimum height threshold for peak detection.
  462. max_height (float): Maximum height threshold for peak detection.
  463. sigma (float): Standard deviation for Gaussian smoothing of discrete pumping rate.
  464. save (bool): Whether to save the pumping rate and peak times to disk.
  465. save_path (str): Directory path to save output files (required if save=True).
  467. Outputs:
  468. pr_cont (ndarray): Continuous pumping rate signal (Gaussian smoothed).
  469. pr_disc (ndarray): Discrete pump counts per time bin.
  470. peaks (ndarray): Frame indices where valid pump peaks occur.
  471. troughs (ndarray): Frame indices where valid pump troughs occur.
  472. """
  473. ### Check that a video path, grinder motion, and grinder CoMs were both provided
  474. if video_path is None: raise ValueError("No video path was provided.")
  475. if grinder_motion is None: raise ValueError("No grinder motion signal was provided.")
  477. ### Detect peaks only if they are followed by a corresponding trough
  478. dist = int(fps/(fc/2)) # minimum number of frames it takes to complete a pump
  479. height_range = [min_height, max_height]
  480. troughs, peaks = find_troughs_and_peaks(grinder_motion, distance=dist, height=height_range)
  482. ### Bin the detected pumps (i.e., peaks followed by troughs)
  483. bins = make_bin_edges(fps,1/fps,grinder_motion.shape[0])
  484. pr_disc, bin_edges = np.histogram(peaks, bins=bins)
  486. ### Filter the binned counts to obtain continuous estimate of pumping rate
  487. pr_cont = gaussian_filter1d(pr_disc.astype('float'), sigma, mode='nearest')
  489. ### OPTIONAL: save motion and pumping rate estimates
  490. if save:
  491. if save_path is None: raise ValueError("No save path was provided.")
  492. video_name = os.path.splitext(os.path.basename(video_path))[0]
  493. save_rate_path = save_path + video_name.replace('_cropped', '') + '_PumpKin_rate.csv'
  494. np.savetxt(save_rate_path, pr_cont)
  495. print('Continuous pumping rate saved to ' + save_rate_path)
  497. save_times_path = save_path + video_name.replace('_cropped', '') + '_PumpKin_times.csv'
  498. np.savetxt(save_times_path, peaks/fps) # saving pump times in seconds, not frames
  499. print('Pumping times saved to ' + save_times_path)
  501. return pr_cont, pr_disc, peaks, troughs

pumpkin.py at commit 35b4184, under MIT · at the source

Overview

Authors: Erin Shappell1,2, Debra Buggs3, Jennah Walcott4, Hang Lu1,3
ORCID iDs: Hang Lu
  1. Interdisciplinary Program in Bioengineering, Georgia Institute of Technology, Atlanta, GeorgiaUnited States of America
  2. School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, GeorgiaUnited States of America
  3. School of Chemical & Biomolecular Engineering, Georgia Institute of Technology, Atlanta, GeorgiaUnited States of America
  4. Coulter Department of Biomedical Engineering, Georgia Institute of Technology, Atlanta, GeorgiaUnited States of America
Institutions: Georgia Institute of Technology (United States); The Wallace H. Coulter Department of Biomedical Engineering (United States)
Journal: PLoS computational biology, volume 22, issue 7, article e1014489
Dates: received 19 January 2026; accepted 22 June 2026; published online 17 July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1371/journal.pcbi.1014489 · PMID 42467749 · PMCID PMC13399524 · OpenAlex W7169503192
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: other (modality), human (organism), C. elegans (organism), methods / tools (subfield)
Methods: Spectral & time-frequency, Machine learning, fMRI & imaging, Single-unit activity, calcium imaging
MeSH: Caenorhabditis elegans*, Machine Learning*, Animals, Behavior, Animal, Biomechanical Phenomena, Computational Biology, Convolutional Neural Networks, Image Processing, Computer-Assisted, Locomotion, Motion Capture, Neural Networks, Computer, Software (* major topic)
Topic: Genetics, Aging, and Longevity in Model Organisms (Aging, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: National Science Foundation (DGE-2039655, 1648035, 1764406); NIA NIH HHS (R01 AG082039); Foundation for the National Institutes of Health (NIH R01AG082039)
Citations: not cited yet (Europe PMC); 42 references in the paper

Abstract

One of the many goals of neuroscience is to understand how the brain encodes and transforms sensory information into behavior. These animal behaviors can be studied at the level of multi-limb poses or through the focused analysis of individual body parts. Techniques for tracking animal pose, such as DeepLabCut and SLEAP, enable detailed studies of large-scale multi-limb behaviors but show reduced accuracy when used for single-keypoint tracking, where insufficient spatial context leads to increased drift and instability in tracking (Arent I, Schmidt FP, Botsch M et al. Marker-less motion capture of insect locomotion with deep neural networks pre-trained on synthetic videos. Frontiers in Behavioral Neuroscience. Vol. 15. 2021. Tang G, Han Y, Sun X, et al. Anti-drift pose tracker (ADPT), a transformer-based network for robust animal pose estimation cross-species. eLife. Vol. 13. 2025). More general techniques, such as Faster Region-based Convolutional Neural Network (Faster R-CNN) and You Only Look Once (YOLO), have also been used to track location-based behaviors such as center-of-mass position and velocity. However, behaviors localized to a single body structure, such as the pharyngeal pumping (i.e., feeding) in the microscopic roundworm Caenorhabditis elegans (C. elegans), are particularly sensitive to noise from moving non-target body parts. This limitation cannot be resolved by simply adding more training data, as doing so often leads to overfitting rather than improved robustness, and instead requires additional processing beyond existing object tracking packages. To address these challenges, we present a fast, automated method that reliably measures pumping in freely moving C. elegans by combining a state-of-the-art object detector (Faster R-CNN) with a tunable noise filter in a technique we call PumpKin. To validate its performance, we demonstrate both its speed (average of 0.4 seconds/frame) and its robust estimation capabilities through application to eight different experimental conditions that encompass both satiety and genetically-driven changes to feeding. PumpKin accurately estimates average pumping rates under eight different experimental conditions, which are positively correlated with the estimates of two expert annotators. Furthermore, PumpKin provides reliable estimates of the instantaneous pumping rate dynamics, achieving an average overlap that exceeds the human–human agreement measured via leave-one-out analysis. Applying PumpKin to conditions differing in satiety revealed a shared basal pumping rate of 0.5 Hz across all worm groups recorded off food, regardless of genetic background or satiety state. Together, these findings highlight PumpKin’s ability to accurately isolate and estimate the motion of a single body part during locomotion. Although we present results specific to C. elegans, we anticipate that PumpKin will generalize to behaviors localized to a single body structure in other systems.

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

Repository

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

lu-lab/PumpKin

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 35b41846face65d5dca3c04ca4ba68f520e9f219, 11 August 2026
Languages: Python (9), Jupyter (5)
Size: 39 files, 14 scripts
Software Heritage: not archived
Found in: “Data Availability”
Holds: README, license file, 5 notebooks
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (5 files), OpenCV (5 files), Matplotlib (4 files), PyTorch (4 files), SciPy (3 files), pandas (2 files), scikit-image (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
16 files
At the source: github.com/lu-lab/PumpKin

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:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 14 scripts, each with its path and the digest of its content;
  • 7 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

Data supporting the findings of this study are available at https://osf.io/79hfv. The PumpKin package is available for download on our GitHub: https://github.com/lu-lab/PumpKin.

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, issue, pages, dates, 4 authors, 12 MeSH terms, 3 funders, 37 references.

Cite

This paper

Shappell, E., Buggs, D., Walcott, J., & Lu, H. (2026). PumpKin: A machine-learning pipeline for automatically tracking localized kinematics in freely moving C. elegans. PLoS computational biology, 22(7), e1014489. https://doi.org/10.1371/journal.pcbi.1014489

BibTeX

@article{shappell2026pumpkin,
author = {Shappell, Erin and Buggs, Debra and Walcott, Jennah and Lu, Hang},
title = {{PumpKin: A machine-learning pipeline for automatically tracking localized kinematics in freely moving C. elegans}},
journal = {PLoS computational biology},
year = {2026},
month = jul,
volume = {22},
number = {7},
pages = {e1014489},
publisher = {PLOS},
issn = {1553-734X},
doi = {10.1371/journal.pcbi.1014489},
url = {https://doi.org/10.1371/journal.pcbi.1014489},
pmid = {42467749},
pmcid = {PMC13399524}
}

RIS

TY - JOUR
AU - Shappell, Erin
AU - Buggs, Debra
AU - Walcott, Jennah
AU - Lu, Hang
TI - PumpKin: A machine-learning pipeline for automatically tracking localized kinematics in freely moving C. elegans
T2 - PLoS computational biology
J2 - PLoS Comput Biol
PY - 2026
DA - 2026/07/17
VL - 22
IS - 7
SP - e1014489
SN - 1553-734X
PB - PLOS
DO - 10.1371/journal.pcbi.1014489
UR - https://doi.org/10.1371/journal.pcbi.1014489
LA - en
ER -

CSL-JSON

{
"id": "10.1371/journal.pcbi.1014489",
"type": "article-journal",
"title": "PumpKin: A machine-learning pipeline for automatically tracking localized kinematics in freely moving C. elegans",
"container-title": "PLoS computational biology",
"author": [
{
"family": "Shappell",
"given": "Erin"
},
{
"family": "Buggs",
"given": "Debra"
},
{
"family": "Walcott",
"given": "Jennah"
},
{
"family": "Lu",
"given": "Hang"
}
],
"container-title-short": "PLoS Comput Biol",
"volume": "22",
"issue": "7",
"page": "e1014489",
"DOI": "10.1371/journal.pcbi.1014489",
"PMID": "42467749",
"PMCID": "PMC13399524",
"ISSN": "1553-734X",
"publisher": "PLOS",
"URL": "https://doi.org/10.1371/journal.pcbi.1014489",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
17
]
]
}
}

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-72057-9 [code]
Sex-specific behavioral feedback modulates sensorimotor processing and drives flexible social behavior.
Journal: Nature communications
In common: OpenCV, scikit-image, PyTorch, 4 other tools, 2 references
[2] doi:10.1038/s41586-026-10348-3 [code]
An enteric neuron ionotropic receptor regulates salt stress resistance.
Journal: Nature
In common: OpenCV, scikit-image, PyTorch, 4 other tools, C. elegans, 1 reference
[3] doi:10.1038/s41467-026-72710-3 [code]
A modular multi-color fluorescence microscope for simultaneous tracking of cellular activity and behavior.
Journal: Nature communications
In common: OpenCV, scikit-image, pandas, 3 other tools, C. elegans, 2 references
[4] doi:10.1038/s41593-026-02262-8 [code]
Cheese3D enables sensitive detection and analysis of whole-face movement in mice.
Journal: Nature neuroscience
In common: OpenCV, scikit-image, pandas, 3 other tools, 3 references
[5] doi:10.1016/j.isci.2026.116825 [code]
Social hierarchy shapes behavioral and transcriptional responses to chronic stress and ketamine in male mice.
Journal: iScience
In common: OpenCV, scikit-image, PyTorch, 4 other tools, 2 references
[6] doi:10.1038/s41593-026-02232-0 [code]
Entorhinal cortex represents task-relevant remote locations independently of CA1.
Journal: Nature neuroscience
In common: OpenCV, scikit-image, PyTorch, 4 other tools, 2 references
[7] doi:10.1038/s41467-026-72709-w [code]
An epifluorescence microscope design for naturalistic behavior and cellular activity in freely moving Caenorhabditis elegans.
Journal: Nature communications
In common: OpenCV, scikit-image, PyTorch, 4 other tools, C. elegans, 1 reference
[8] doi:10.1371/journal.pcbi.1013441 [code]
Large vision model framework for automated C. elegans analysis: From static morphometry to dynamic neural activity.
Journal: PLoS computational biology
In common: OpenCV, scikit-image, PyTorch, 4 other tools, C. elegans, methods / tools
[9] doi:10.1038/s41467-026-73476-4 [code]
Developmental molecular signatures define de novo cortico-brainstem circuit for skilled forelimb movement.
Journal: Nature communications
In common: OpenCV, scikit-image, PyTorch, 3 other tools, 2 references
[10] doi:10.1126/sciadv.aed3650 [code]
Truthful visualizations for mass spectrometry imaging enable high-spatial-resolution interactive &lt;i&gt;m/z&lt;/i&gt; mapping and exploration.
Journal: Science advances
In common: OpenCV, scikit-image, PyTorch, 4 other tools, other, methods / 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.