OSCR

Embodied behavioural complexity in a ciliated microorganism.

Code ↔ Paper

4 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 4 matches
  1. [1] § Methods › Wavelet analysis and multidimensional scaling analysis ↔ wavelet_analysis.jl, lines 189–253 · score 0.70 · overlapping points, wavelet vector, MDS coordinates, Euclidean, distances
  2. [2] § Results › Active ciliary bundle oscillations drive distinct behaviours ↔ Example_Codes.py, lines 97–166 · score 0.69 · Ciliary waveforms, tangent angle, arc length, cell body, axis, frame
  3. [3] § Methods › Wavelet analysis and multidimensional scaling analysis ↔ wavelet_analysis.jl, lines 128–187 · score 0.55 · dispersion heatmap, wavelet transform, multidimensional
  4. [4] § Methods › Cilium tracking and reconstruction ↔ Example_Codes.py, lines 48–94 · score 0.53 · Cartesian space, arc length, tangent, tracking, reconstruction, Cilium

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

Julia · 438 lines · 19 KB · CC-BY-4.0 · 2 matches

  1. cd(@__DIR__); import Pkg; Pkg.activate("")
  2. # Load relevant packages
  3. using LinearAlgebra, Statistics
  4. using ClassicalOrthogonalPolynomials, FFTW, RecurrenceRelationships
  5. using Optim
  6. using DelimitedFiles
  7. using ProgressMeter
  8. using CairoMakie, GLMakie; CairoMakie.activate!()
  9. using Roots
  10. using Metaheuristics
  11. using FFTW
  12. using Distances, Clustering, MultivariateStats
  13. import TOML
  15. include("chebyshev.jl")
  16. include("utils.jl")
  17. include("wavelet.jl")
  18. include("plotting.jl")
  20. ## Analysis parameters
  21. cwtfreqs = 2.0 .^LinRange(log2(5), log2(300), 101) # Log spaced frequencies to evaluate the wavelet transform at
  22. f = 1.0 # Morelet scale factor
  23. max_chebind = 20 # Maximum chebyshev index
  24. crplen = 20 # Length to crop wavelet transform at the boundaries
  26. # Experiments to pick out as examples
  27. target_names = ["F26(6)", "F20(4)", "F20(29)", "F20(16)", "F20(20)", "F20(25)", "F15(16)", "F7(4)", "F46(71)", "F27(4)", "F15(30)"]
  29. ##
  30. # Read in the data
  31. all_data = read_coefficients("LocalDataPath.toml")
  32. # Get all the nonstationary datasets
  33. nonstatInds = unique(push!(findall(.!in(["Stationary"]), all_data.datasettypes), findall(all_data.datasetnames .∈ (target_names, ))...))
  34. nonstat_data = all_data[nonstatInds]
  35. target_data = nonstat_data[target_names]
  36. # Get all the stationary datasets
  37. statInds = findall(in(["Stationary"]), all_data.datasettypes)
  38. stat_data = all_data[statInds]
  40. # Evaluate the wavelet transform
  41. nonstat_wavelets = CiliaWaveletTransforms(nonstat_data, cwtfreqs, f; max_cwtind = max_chebind, fps = 1000, zp = true)
  42. target_wavelets = CiliaWaveletTransforms(target_data, cwtfreqs, f; max_cwtind = max_chebind, fps = 1000, zp = true)
  44. ## Calculate the disperion plot
  45. cwtfreqs = nonstat_wavelets[1].cwtfreqs
  46. dispmat = Matrix{Float64}(undef, length(cwtfreqs), max_chebind)
  47. for i = 2:max_chebind + 1
  48. dispmat[:, i - 1] = mean(vcat([cdata.abs_wavelet_transform[crplen + 1:end - crplen, :, i] for cdata in nonstat_wavelets]..., target_wavelets[8].abs_wavelet_transform[crplen + 1:end - crplen, :, i]), dims = 1)[:]
  49. end
  50. # normalized version of the plot
  51. normdispmat = dispmat ./ sqrt.(sum(abs2, dispmat, dims = 1))
  52. # Calculate the Chebyshev wavenumber
  53. mL = Statistics.median(vcat(nonstat_data.Ls...))
  54. chebwavlen = 0.370*mL*[2; [mean(abs.(cos.(((0:n - 2) .+ 0.5) * pi / n) * (cos.(pi/n) - 1) .- sin.(((0:n - 2) .+ 0.5) * pi / n) * sin(pi/n))) for n = 2:max_chebind]]
  55. chebwavnum = 2*pi ./ chebwavlen
  57. # Get the trajectories on the wavelet transform
  58. dispkind = "softmax"
  59. if dispkind === "softmax"
  60. target_disptrajs = dispersion_trajectory_softmax.(target_wavelets, (chebwavnum, ), filtlen = 50)
  61. nonstat_disptrajs = dispersion_trajectory_softmax.(nonstat_wavelets, (chebwavnum, ), filtlen = 50)
  62. elseif dispkind === "mean"
  63. target_disptrajs = dispersion_trajectory.(target_wavelets, (chebwavnum, ))
  64. nonstat_disptrajs = dispersion_trajectory.(nonstat_wavelets, (chebwavnum, ))
  65. else
  66. error("Oops")
  67. end
  69. target_disptrajs_mean = dispersion_trajectory.(target_wavelets, (chebwavnum, ))
  70. nonstat_disptrajs_mean = dispersion_trajectory.(nonstat_wavelets, (chebwavnum, ))
  72. target_disptrajs_softmax = dispersion_trajectory_softmax.(target_wavelets, (chebwavnum, ), filtlen = 50)
  73. nonstat_disptrajs_softmax = dispersion_trajectory_softmax.(nonstat_wavelets, (chebwavnum, ), filtlen = 50)
  74. #target_disptrajs = dispersion_trajectory.(target_wavelets, (chebwavnum, ))
  75. #nonstat_disptrajs = dispersion_trajectory.(nonstat_wavelets, (chebwavnum, ))
  77. #= Uncomment if you want to save the analysis results
  78. # Save the results
  79. savelength = maximum(length.(getindex.(target_disptrajs, 1)))
  80. for i = 1:length(target_dispfs)
  81. push!(target_disptrajs[i][1], NaN*ones(savelength - length(target_disptrajs[i][1]))...)
  82. push!(target_disptrajs[i][2], NaN*ones(savelength - length(target_disptrajs[i][1]))...)
  83. end
  84. writedlm("dispersion_trajectories_v2.csv", [reshape([reshape(target_names, 1, :); reshape(target_names, 1, :)][:], 1, :); repeat(["wavenumber" "freq"], 1, 11); reshape([hcat(getindex.(target_disptrajs, 1)...); hcat(getindex.(target_disptrajs, 2)...)], savelength, :)], ',')
  85. writedlm("dispersion_04_08_25.csv", [["freq (Hz) / chebwavnum (um^-1)"; nonstat_wavelets.cwtfreqs] [chebwavnum'; dispmat]], ',')
  86. =#
  88. ##
  89. # Plot the result for the dispersion subsamples
  90. fig = Figure(size = (1000, 480), fontsize = 24)
  91. # Plot the dispersion heatmap
  92. LOGAX = false; LOGCB = false
  93. if LOGAX
  94. ax = Axis(fig[1, 1], aspect = AxisAspect(1), yscale = log10); ax.xlabel = "k (μm⁻¹)"; ax.ylabel = "f (Hz)" # yscale = log10
  95. else
  96. ax = Axis(fig[1, 1], aspect = AxisAspect(1)); ax.xlabel = "k (μm⁻¹)"; ax.ylabel = "f (Hz)" # yscale = log10
  97. end
  98. if LOGCB
  99. crange = (-2.35, -0.5)
  100. pltmat = log10.(normdispmat')
  101. hmp = heatmap!(ax, chebwavnum, nonstat_wavelets.cwtfreqs, pltmat, colormap = :nuuk, colorrange = crange)
  102. cb = Colorbar(fig[1, 2], hmp, label = "Log10 Normalized power")
  103. else
  104. crange = (0.0, 0.3)
  105. pltmat = normdispmat'
  106. hmp = heatmap!(ax, chebwavnum, nonstat_wavelets.cwtfreqs, pltmat, colormap = :nuuk, colorrange = crange)
  107. cb = Colorbar(fig[1, 2], hmp, label = "Normalized power")
  108. end
  109. # Reconstructions
  110. xs = LinRange(-1, 1, 1001)
  111. V = zeros(length(xs), 126); eval_vandermonde!(V, xs)
  112. for i = 1:9
  113. dispks = target_disptrajs[i][1]
  114. dispfs = target_disptrajs[i][2]
  115. crange = (1, size(target_data[i].allcfsx, 1))
  116. lines!(ax, dispks, dispfs, colormap = cmaps[i][2:end], color = 1:length(dispks), linewidth = 1.5, colorrange = crange)
  117. xrecon = V[:, 2:end] * target_data[i].allcfsx[:, 2:end]'
  118. yrecon = V[:, 2:end] * target_data[i].allcfsy[:, 2:end]'
  119. axtmp = Axis(fig[1, 3][divrem(i - 1, 3) .+ 1...], aspect = DataAspect()); hidedecorations!(axtmp)
  120. ylims!(axtmp, -100, 100); xlims!(axtmp, -100, 100)
  121. nplt = div(size(target_data[i].allcfsx, 1), 25) # Plot every 25 frames
  122. for j = size(target_data[i].allcfsx, 1):-nplt:1
  123. lines!(axtmp, Point2.(xrecon[:, j], yrecon[:, j]), color = j, colormap = cmaps[i][2:end], colorrange = crange, linewidth = 0.5)
  124. end
  125. end
  126. display(fig)
  128. ## Plot the results corresponding to the SI movies
  129. # Plot the result for the dispersion subsamples
  130. fig = Figure(size = (1000, 480), fontsize = 24)
  131. # Plot the dispersion heatmap
  132. ax = Axis(fig[1, 1], aspect = AxisAspect(1), yscale = log10); ax.xlabel = "k (μm⁻¹)"; ax.ylabel = "f (Hz)" # yscale = log10
  133. hmp = heatmap!(ax, chebwavnum, nonstat_wavelets.cwtfreqs, normdispmat', colormap = :nuuk, colorrange = (0, 0.3))
  134. cb = Colorbar(fig[1, 2], hmp, label = "Normalized power")
  135. # Reconstructions
  136. xs = LinRange(-1, 1, 1001)
  137. V = zeros(length(xs), 126); eval_vandermonde!(V, xs)
  138. cmapind = [1, 6]
  139. for i2 = 1:2
  140. i = i2 + 9
  141. dispks = target_disptrajs[i][1]
  142. dispfs = target_disptrajs[i][2]
  143. crange = (1, size(target_data[i].allcfsx, 1))
  144. lines!(ax, dispks, dispfs, colormap = cmaps[cmapind[i2]][2:end], color = 1:length(dispks), linewidth = 1.5, colorrange = crange)
  145. xrecon = V[:, 2:end] * target_data[i].allcfsx[:, 2:end]'
  146. yrecon = V[:, 2:end] * target_data[i].allcfsy[:, 2:end]'
  147. axtmp = Axis(fig[1, 3][i2, 1], aspect = DataAspect()); hidedecorations!(axtmp)
  148. ylims!(axtmp, -100, 100); xlims!(axtmp, -100, 100)
  149. nplt = div(size(target_data[i].allcfsx, 1), 25) # Plot every 25 frames
  150. for j = size(target_data[i].allcfsx, 1):-nplt:1
  151. lines!(axtmp, Point2.(xrecon[:, j], yrecon[:, j]), color = j, colormap = cmaps[cmapind[i2]][2:end], colorrange = crange, linewidth = 0.5)
  152. end
  153. end
  154. display(fig)
  156. ## Multidimensional scaling
  157. reorientation_mats = [hcat([cwav.abs_wavelet_transform[crplen + 1:end - crplen, :, i] for i = 2:max_chebind + 1]...) for cwav in nonstat_wavelets]
  158. for rmat = reorientation_mats
  159. rmat ./= sum(rmat, dims = 2)
  160. end
  162. # Split the data and perform MDS on each split
  163. using Random; Random.seed!(123456)
  164. # Total number of time points
  165. n = sum(size.(reorientation_mats, 1))
  166. # Map each time point to a unique linear time index (index of coli in expi = offset[i] + coli)
  167. offsets = [0; cumsum(size.(reorientation_mats, 1))]
  168. # Overlap
  169. c = 1666
  170. # Number of independent MSDs
  171. nfolds = 3
  172. m, r = divrem(n - c, nfolds)
  173. l = m + c
  174. # Maybe there is a small fold is not a perfect divisor
  175. np = ceil(Int, 1 + (n - l)/(l - c))
  176. # Generate a random permutation of the data
  177. rp = randperm(n)
  178. # Get the columns in each fold
  179. partitionInds = [rp[1:l], [rp[l + (i - 1)*(l - c) + 1:l + (i)*(l - c)] for i = 1:np - 2]..., rp[l + (np - 2)*(l - c) + 1:end]]
  180. # Map each linear time index to an experiment ind
  181. expindvec = vcat(map(i ->i * ones(Int64, size(reorientation_mats[i], 1)), 1:length(reorientation_mats))...)
  182. # Map each linear time index to an experiment colum
  183. colindvec = vcat(map(i -> 1:size(reorientation_mats[i], 1), 1:length(reorientation_mats))...)
  184. # Map each index in the fold to a partition
  185. partindvec = [ones(Int, l); [i*ones(Int, l - c) for i = 2:np]...]
  186. # Map each index in the folds to a column in a given partition
  187. partcolindvec = [1:l; [1:(l - c) for i = 2:np]...]
  189. # A note on indexing: the ith column of partition j in the splits is mapped to experiment expindvec[partitionInds[j][i]] and row within that experiment colindvec[partitionInds[j][i]]
  191. # Build the wavelet vectors in each fold
  192. X1 = zeros(size(reorientation_mats[1], 2), l)
  193. for (ind, i) in enumerate(partitionInds[1])
  194. X1[:, ind] = reorientation_mats[expindvec[i]][colindvec[i], :]
  195. end
  196. # Choose the overlap points randomly from the first fold
  197. cinds = randperm(l)[1:c]
  198. Xc = X1[:, cinds]
  199. Xps = [X1, ]
  200. # Loop to extract each of the partition matrics
  201. for j = 2:np - 1
  202. Xj = zeros(size(reorientation_mats[1], 2), l)
  203. for (ind, i) in enumerate(partitionInds[j])
  204. Xj[:, ind] = reorientation_mats[expindvec[i]][colindvec[i], :]
  205. end
  206. # Append the overlap region
  207. Xj[:, l - c + 1:end] .= Xc
  208. push!(Xps, Xj)
  209. end
  210. # Build the possibly small excess set of points
  211. Xj = zeros(size(reorientation_mats[1], 2), length(partitionInds[np]) + c)
  212. for (ind, i) in enumerate(partitionInds[np])
  213. Xj[:, ind] = reorientation_mats[expindvec[i]][colindvec[i], :]
  214. end
  215. # Append the overlap region
  216. Xj[:, length(partitionInds[np]) + 1:end] .= Xc
  217. push!(Xps, Xj)
  219. ## Calculate the pairwise distance matrices
  220. Dps = similar(Xps)
  221. @time for i = 1:np
  222. Dps[i] = pairwise(SqEuclidean(1e-8), Xps[i])
  223. println("Evaluated distance matrix $(i) / $(np)")
  224. end
  226. # Calculate the MDS on each split
  227. λs = Vector{Float64}[]
  228. Vs = Matrix{Float64}[]
  229. Xs = Matrix{Float64}[]
  230. # With the current set up this takes about 8 minutes to run! Increase the number of folds to speed up
  231. @time for i = 1:np
  232. # Current MDS Matrix
  233. D2 = Dps[i]
  234. # Calculate the Gram matrix
  235. G = 0.5 * (-D2 .+ mean(D2, dims = 1) .+ mean(D2, dims = 2) .- mean(D2))
  236. # Eigen decomposition of the Gram matrix
  237. E = eigen(Symmetric(G))
  238. # Make sure the eigenvalue are ordered appropriately
  239. sp = sortperm(E.values, rev = true)
  240. λ = E.values[sp]
  241. # Flip if necessary
  242. if abs(λ[end]) > abs(λ[1])
  243. λ = -λ[end:-1:1]
  244. sp = reverse(sp)
  245. end
  246. # We keep save first 5 MDS coordinates
  247. V = E.vectors[:, sp[1:5]]
  248. X = V' .* sqrt.(λ[1:5])
  249. push!(λs, λ)
  250. push!(Vs, V)
  251. push!(Xs, X)
  252. println("Evaluated MDS $(i) / $(np)")
  253. end
  255. # Plot the variance explained
  256. fig = Figure()
  257. ax = Axis(fig[1, 1])
  258. for i = 1:np
  259. scatter!(ax, cumsum(abs2.(λs[i]) ./ sum(abs2, λs[i]))[1:14], color = plotGray)
  260. end
  261. ax.xticks = 1:14
  262. xlims!(ax, 0.5, 14.5)
  263. lines!(ax, [0.5, 14.5], [0.9, 0.9], color = plotGray, linestyle = :dash)
  264. lines!(ax, [0.5, 14.5], [0.95, 0.95], color = plotGray, linestyle = :dash)
  265. lines!(ax, [0.5, 14.5], [0.99, 0.99], color = plotGray, linestyle = :dash)
  266. mλ = mean(λs[1:end - 1])
  267. scatter!(ax, cumsum(abs2.(mλ) ./ sum(abs2, mλ)), color = plotBlue, markersize = 25)
  268. ax.xlabel = "Number of dimensions"
  269. ax.ylabel = "Variance explained"
  270. #save("/Users/alaasdairhastewell/Dropbox (MIT)/Research/Pterosperma/Pterosperma Plots/all_attactor_dimension.pdf", fig)
  271. display(fig)
  273. # Align the seperate MDS results
  274. Xtarget = Xs[1][:, cinds]' # X values from
  275. mTarget = mean(Xtarget, dims = 1)
  276. sclXtarget = Xtarget .- mTarget
  277. Xs_aligned = [Xs[1], ]
  278. for i = 2:np
  279. Xcur = Xs[i][:, size(Xs[i], 2) - c + 1:end]'
  280. mCur = mean(Xcur, dims = 1)
  281. #sclXcur = Xcur .- mCur
  282. U, S, V = svd(sclXtarget'*Xcur)
  283. rotMat = V*U'
  284. transVec = mTarget[:] - rotMat' * mCur[:]
  285. Xtrans = Xs[i][:, 1:size(Xs[i], 2) - c]'
  286. XcurAligned = Matrix((Xtrans * rotMat .+ transVec')')
  287. push!(Xs_aligned, XcurAligned)
  288. @show sqrt(mean(abs2, (Xcur * rotMat .+ transVec') - Xtarget)) / sqrt(mean(abs2, Xtarget))
  289. end
  291. # Get the MDS coordinates for each experiment
  292. globalMDScoords = Matrix{Float64}[]
  293. for expIndCur = 1:length(reorientation_mats)
  294. tmp = zeros(size(Xs_aligned[1], 1), size(reorientation_mats[expIndCur], 1))
  295. for (ind, partinds) in enumerate(partitionInds)
  296. for (ind2, i) in enumerate(partinds)
  297. if expindvec[i] == expIndCur
  298. tmp[:, colindvec[i]] .= Xs_aligned[ind][:, ind2]
  299. end
  300. end
  301. end
  302. push!(globalMDScoords, tmp)
  303. end
  305. ## Color by MDS split
  306. using GLMakie; GLMakie.activate!()
  307. pltaz = 0.4809765624999961
  308. pltel = 0.32449374999999997
  309. #CairoMakie.activate!()
  310. fig = Figure(size = (800, 500), fontsize = 15)
  311. ax = Axis3(fig[1, 1], aspect = :data, xgridvisible = false, zgridvisible = false, ygridvisible = false)
  312. for i = 1:np
  313. scatter!(ax, -Xs_aligned[i][1, :], -Xs_aligned[i][2, :], -Xs_aligned[i][3, :], color = cmaps[[1, 5, 9][i]][end], markersize = 2)
  314. end
  315. Legend(fig[1, 2], [MarkerElement(color = cmaps[1][end], markersize = 15, marker = :circle), MarkerElement(color = cmaps[5][end], markersize = 15, marker = :circle), MarkerElement(color = cmaps[9][end], markersize = 15, marker = :circle)], ["Split 1", "Split 2", "Split 3"], framevisible = false)
  316. ax.xlabel = "MDS 1"
  317. ax.ylabel = "MDS 2"
  318. ax.zlabel = "MDS 3"
  319. ax.elevation = pltel
  320. ax.azimuth = pltaz
  321. #save("MDS_split.pdf", fig)
  322. fig
  324. ## Color by wavenumber
  325. fig = Figure(size = (800, 500), fontsize = 15)
  326. ax = Axis3(fig[1, 1], aspect = :data, xgridvisible = false, zgridvisible = false, ygridvisible = false)
  327. for i = 1:length(globalMDScoords)
  328. scatter!(ax, -globalMDScoords[i][1, :], -globalMDScoords[i][2, :], -globalMDScoords[i][3, :], color = nonstat_disptrajs[i][1][crplen + 1:end - crplen], markersize = 2, colorrange = (0.05, 0.2), colormap = :bamako)
  329. end
  330. Colorbar(fig[1, 2], colorrange = (0.025, 0.075), label = "Wavenumber μm⁻¹", colormap = :bamako, height = Relative(0.5))
  331. ax.xlabel = "MDS 1"
  332. ax.ylabel = "MDS 2"
  333. ax.zlabel = "MDS 3"
  334. ax.elevation = pltel
  335. ax.azimuth = pltaz
  336. #save("MDS_wavenumber.pdf", fig)
  337. fig
  339. ## Color by the frequency
  340. fig = Figure(size = (800, 500), fontsize = 15)
  341. ax = Axis3(fig[1, 1], aspect = :data, xgridvisible = false, zgridvisible = false, ygridvisible = false)
  342. for i = 1:length(globalMDScoords)
  343. scatter!(ax, -globalMDScoords[i][1, :], -globalMDScoords[i][2, :], -globalMDScoords[i][3, :], color = nonstat_disptrajs[i][2][crplen + 1:end - crplen], markersize = 2, colorrange = (6, 125), colormap = :batlow)
  344. end
  345. Colorbar(fig[1, 2], colorrange = (6, 125), label = "Frequency (Hz)", colormap = :batlow, height = Relative(0.5))
  346. ax.xlabel = "MDS 1"
  347. ax.ylabel = "MDS 2"
  348. ax.zlabel = "MDS 3"
  349. ax.elevation = pltel
  350. ax.azimuth = pltaz
  351. #save("MDS_frequency.pdf", fig)
  352. fig
  353. fig
  355. ## Color by MDS4
  356. fig = Figure(size = (800, 500), fontsize = 15)
  357. ax = Axis3(fig[1, 1], aspect = :data, xgridvisible = false, zgridvisible = false, ygridvisible = false)
  358. for i = 1:length(globalMDScoords)
  359. scatter!(ax, -globalMDScoords[i][1, :], -globalMDScoords[i][2, :], -globalMDScoords[i][3, :], color = globalMDScoords[i][4, :], markersize = 2, colorrange = (-0.075, 0.075), colormap = :bam)
  360. end
  361. Colorbar(fig[1, 2], colorrange = (-0.075, 0.075), label = "MDS 4", colormap = :bam, height = Relative(0.5))
  362. ax.xlabel = "MDS 1"
  363. ax.ylabel = "MDS 2"
  364. ax.zlabel = "MDS 3"
  365. ax.elevation = pltel
  366. ax.azimuth = pltaz
  367. #save("MDS_mds4.pdf", fig)
  368. fig
  370. ##
  371. # Equal count Histogram
  372. num_per_bin = 1000
  373. indata = vcat(getindex.(nonstat_disptrajs, 1)...)
  374. sp = sortperm(indata)
  375. edgeinds = sp[1:num_per_bin:end]
  376. binedges = indata[edgeinds]
  377. binheights = (num_per_bin / length(sp)) ./ diff(binedges)
  378. if (sp[1:num_per_bin:end])[end] != sp[end]
  379. tmpcount = length(sp) - findfirst(==(edgeinds[end]), sp)
  380. push!(binedges, indata[sp[end]])
  381. push!(binheights, (tmpcount / length(sp)) / (binedges[end] - binedges[end - 1]))
  382. end
  383. fig = Figure()
  384. ax = Axis(fig[1, 1])
  385. barplot!(ax, (binedges[1:end - 1] .+ binedges[2:end])/2, binheights, strokewidth = 0, gap = 0, width = (binedges[2:end] .- binedges[1:end - 1]), color = plotGray)
  386. ylims!(ax, 0.0, 20.0)
  387. xlims!(ax, 0.0, 0.55)
  388. ax.xlabel = "Wavenumber μm⁻¹"
  389. ax.ylabel = "PDF"
  390. #save("wavenumer_histogram.pdf", fig)
  391. bincenters_wavenumber = (binedges[1:end - 1] + binedges[2:end])/2
  392. binheights_wavenumber = copy(binheights)
  393. fig
  395. ##
  396. # Equal count Histogram
  397. num_per_bin = 1000
  398. indata = vcat(getindex.(nonstat_disptrajs, 2)...)
  399. sp = sortperm(indata)
  400. edgeinds = sp[1:num_per_bin:end]
  401. binedges = indata[edgeinds]
  402. binheights = (num_per_bin / length(sp)) ./ diff(binedges)
  403. if (sp[1:num_per_bin:end])[end] != sp[end]
  404. tmpcount = length(sp) - findfirst(==(edgeinds[end]), sp)
  405. push!(binedges, indata[sp[end]])
  406. push!(binheights, (tmpcount / length(sp)) / (binedges[end] - binedges[end - 1]))
  407. end
  408. fig = Figure()
  409. ax = Axis(fig[1, 1])
  410. barplot!(ax, (binedges[1:end - 1] .+ binedges[2:end])/2, binheights, strokewidth = 0, gap = 0, width = (binedges[2:end] .- binedges[1:end - 1]), color = plotGray)
  411. ylims!(ax, 0.0, 0.025)
  412. xlims!(ax, 0.0, 245)
  413. ax.xlabel = "Frequency Hz"
  414. ax.ylabel = "PDF"
  415. bincenters_frequency = (binedges[1:end - 1] + binedges[2:end])/2
  416. binheights_frequency = copy(binheights)
  417. #save("frequency_histogram.pdf", fig)
  418. fig
  420. ##
  421. hist_savemat = [reshape(["Freq bin center", "Freq bin count", "Wavenum bin center", "Wavenum bin count"], 1, :); [bincenters_frequency binheights_frequency bincenters_wavenumber binheights_wavenumber]]
  422. writedlm("dispersion_trajectory_histograms_$(dispkind)_$(today()).csv", hist_savemat, ',')
  424. mds_savemat
  425. labels = ["Experiment"; "Frame"; ("MDS" .* string.(1:5))...; "partition"; "softmax_ks"; "softmax_fs"; "mean_ks"; "mean_fs"]
  426. disptrajmat = [vcat(map(x -> x[1][crplen + 1:end - crplen], nonstat_disptrajs_softmax)...) vcat(map(x -> x[2][crplen + 1:end - crplen], nonstat_disptrajs_softmax)...) vcat(map(x -> x[1][crplen + 1:end - crplen], nonstat_disptrajs_mean)...) vcat(map(x -> x[2][crplen + 1:end - crplen], nonstat_disptrajs_mean)...)]
  427. exppartvec = zeros(size(disptrajmat, 1))
  428. exppartvec[partitionInds[1]] .= 1
  429. exppartvec[partitionInds[2]] .= 2
  430. exppartvec[partitionInds[3]] .= 3
  432. mds_savemat = [reshape(labels, 1, :); [nonstat_data.datasetnames[expindvec] colindvec vcat(globalMDScoords'...) exppartvec disptrajmat]]
  435. writedlm("mds_variance_explained_$(today()).csv", cumsum(abs2.(mλ) ./ sum(abs2, mλ)), ',')
  437. writedlm("mds_coordinates_with_metadata_$(today()).csv", mds_savemat, ',')

wavelet_analysis.jl, under CC-BY-4.0 · at the source

Overview

Authors: Alexander K Boggon1, Alasdair D Hastewell2,3, Jörn Dunkel3, Kirsty Y Wan1
ORCID iDs: Alexander K Boggon, Alasdair D Hastewell, Kirsty Y Wan
  1. Living Systems Institute & Department of Mathematics and Statistics, University of Exeter, Exeter, UK
  2. NSF-Simons National Institute for Theory and Mathematics in Biology, Chicago, IL, USA
  3. Department of Mathematics, Massachusetts Institute of Technology, Cambridge, MA USA
Institutions: University of Exeter (United Kingdom); Massachusetts Institute of Technology (United States)
Journal: Nature communications, volume 17, issue 1, article 8445
Dates: received 24 December 2025; accepted 19 June 2026; published online 8 July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41467-026-75076-8 · PMID 42420268 · PMCID PMC13478290 · OpenAlex W4413470643
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.08.20.671109
Status: code verified
Methods: Spectral & time-frequency
Keywords: Biological physics, Cell migration
MeSH: Cilia*, Movement (* major topic)
Topic: Protist diversity and phylogeny (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: EC | Horizon 2020 Framework Programme (853560); Schmidt Sciences, LLC (Polymath award) MIT MathWorks Professorship Fund; NSF-Simons National Institute for Theory and Mathematics in Biology (NITMB) Fellowship supported via grants from the NSF
Citations: cited by 1 paper (Europe PMC); 74 references in the paper

Abstract

Most animals coordinate behaviour using neural computations. Yet, single-celled organisms also exhibit stimulus-responsive, even cognitive, actions. To understand how a single cell can coordinate and drive complex behaviours without any neural encoding, we study an algal protist – a motile cell with four extremely long cilia. The organism displays a surprisingly rich locomotor repertoire, emerging from the intricate dynamics of the cilia, which form a tight bundle when swimming. We use high-speed quantitative live imaging to extract the spectrum of possible ciliary beating patterns and derive a dispersion relation coupling the temporal frequency and spatial wavelength of cilia oscillations. We further reconstruct the manifold embedded in the behavioural space, showing that despite the range and complexity of ciliary beating modes, the underlying behavioural manifold is intrinsically low-dimensional with non-trivial topological structure. Dynamic transitions in motility patterns are encoded as trajectories in this space, thus reducing the macroscopic behavioural states to the underlying microscopic dynamics of the cilia.

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 4 matches between paragraphs and lines of code.

figshare 29941784

License: CC-BY-4.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Languages: Julia (6), Python (1)
Size: 14 files, 7 scripts
Software Heritage: not checked
Found in: “Code availability”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Makie (2 files), h5py (1 file), Matplotlib (1 file), NumPy (1 file), pandas (1 file), SciPy (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
8 files
At the source:

Code availability

Custom tracking and analysis codes associated with this manuscript are found at 10.6084/m9.figshare.29941784.

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

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;
  • 7 scripts, each with its path and the digest of its content;
  • 4 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

No dataset and no data link were found in the paper.

Data availability

All datasets necessary to reproduce the work detailed in this manuscript are found at 10.6084/m9.figshare.29941784. Source Data associated with all main and Supplementary Figs. are provided with this paper.

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, 2 keywords, 2 MeSH terms, 3 funders, 55 references.

Cite

This paper

Boggon, A. K., Hastewell, A. D., Dunkel, J., & Wan, K. Y. (2026). Embodied behavioural complexity in a ciliated microorganism. Nature communications, 17(1), 8445. https://doi.org/10.1038/s41467-026-75076-8

BibTeX

@article{boggon2026embodied,
author = {Boggon, Alexander K and Hastewell, Alasdair D and Dunkel, Jörn and Wan, Kirsty Y},
title = {{Embodied behavioural complexity in a ciliated microorganism}},
journal = {Nature communications},
year = {2026},
month = jul,
volume = {17},
number = {1},
pages = {8445},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-75076-8},
url = {https://doi.org/10.1038/s41467-026-75076-8},
pmid = {42420268},
pmcid = {PMC13478290}
}

RIS

TY - JOUR
AU - Boggon, Alexander K
AU - Hastewell, Alasdair D
AU - Dunkel, Jörn
AU - Wan, Kirsty Y
TI - Embodied behavioural complexity in a ciliated microorganism
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/07/08
VL - 17
IS - 1
SP - 8445
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-75076-8
UR - https://doi.org/10.1038/s41467-026-75076-8
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-75076-8",
"type": "article-journal",
"title": "Embodied behavioural complexity in a ciliated microorganism",
"container-title": "Nature communications",
"author": [
{
"family": "Boggon",
"given": "Alexander K"
},
{
"family": "Hastewell",
"given": "Alasdair D"
},
{
"family": "Dunkel",
"given": "Jörn"
},
{
"family": "Wan",
"given": "Kirsty Y"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "8445",
"DOI": "10.1038/s41467-026-75076-8",
"PMID": "42420268",
"PMCID": "PMC13478290",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-75076-8",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
8
]
]
}
}

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.1126/sciadv.adv3770 [code]
Evolution of a central dopamine circuit underlies adaptation of a light-evoked sensorimotor response in the blind cavefish.
Journal: Science advances
In common: h5py, pandas, SciPy, 2 other tools, 1 reference
[2] doi:10.1162/imag.a.1341 [code]
Massively parallelized brain tractography using compute clusters, supercomputers, and graphics processing units.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Makie, pandas, SciPy, 2 other tools
[3] doi:10.1038/s41593-026-02342-9 [code]
Interpretable abstractions of artificial neural networks predict behavior and neural activity during human information gathering.
Journal: Nature neuroscience
In common: Makie, pandas, SciPy, 2 other tools
[4] doi:10.1002/mrm.70378 [code]
Investigating the Sensitivity of the Diffusion MRI Signal to Magnetization Transfer and Permeability via Monte-Carlo Simulations.
Journal: Magnetic resonance in medicine
In common: Makie, pandas, SciPy, 2 other tools
[5] doi:10.7554/elife.89629 [code]
Active dendrites enable robust spiking computations despite timing jitter.
Journal: eLife
In common: Makie, pandas, Matplotlib, 1 other tool
[6] doi:10.1038/s41467-026-75352-7 [code]
Mechanosensory encoding of surface mechanics optimizes locomotion.
Journal: Nature communications
In common: pandas, SciPy, Matplotlib, 1 other tool, 1 reference
[7] doi:10.1038/s41467-026-76522-3 [code]
Transcriptomic analysis of spinal V1 interneurons informs their multifunctional role in motor output.
Journal: Nature communications
In common: pandas, SciPy, Matplotlib, 1 other tool, 1 reference
[8] doi:10.1002/advs.77986 [code]
DUET-seq: An Open-Source Droplet Platform for High-Fidelity Joint Chromatin and Transcriptome Profiling Reveals Temporal Regulatory Decoupling in Single Cells.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: h5py, pandas, SciPy, 2 other tools
[9] doi:10.7554/elife.110588 [code]
Opening the black box toward a modular approach to spike sorting.
Journal: eLife
In common: h5py, pandas, SciPy, 2 other tools
[10] doi: [code]
Real-time closed-loop feedback system for mouse mesoscale cortical signal and movement control
Journal: eLife
In common: h5py, pandas, SciPy, 2 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.