Question regarding max_velocity_detector and velocity correction

[antortjim]

Dear Gilestro lab

Thank you for developing sleepr and making it publicly available for free. It’s great!

I had a question regarding your max_velocity_detector annotation function, which you use in sleep_annotation to determine the level of activity of a fly in time_window_length seconds.
The logic you implement consists of retrieving the distance computed live on the ethoscopes, and you then divide by the time passed since the previous frame. This time difference is expressed in the dt column, which you compute in the beginning of the function

sleepr/R/motion_detectors.R

Line 44 in 59949e3

d[,dt := c(NA,diff(t))]

With the distance and the time difference, you then can compute a velocity

sleepr/R/motion_detectors.R

Line 47 in 59949e3

d[,velocity := dist/dt]

In the next lines you simulate beam crosses and mask behavior if an interaction occurs, which is not relevant for my question.
The next relevant line is the following:

sleepr/R/motion_detectors.R

Line 74 in 59949e3

d[, velocity_corrected := velocity * dt /a]

Which according to the supplementary material in your PLOS One article https://journals.plos.org/plosbiology/article/file?id=10.1371/journal.pbio.2003026.s005&type=supplementary is designed to correct for the problem

that resulting velocity computation depended on the frame rate of the processed video

i.e. you correct the velocity taking into account the FPS.
However, by having a close look to the code, I would argue that what you effectively do is instead getting the distance again and then normalizing with a. The result of this operation is called velocity_corrected, but it is basically a distance normalized by a. You then select the maximum value in time_window_length and if the maximum is bigger than 1 i.e. the original distance dist is greater than a, the fly is annotated as moving for that block of time, and not moving otherwise.
Wouldn’t it be more correct to omit line 74 and just check if the velocity is greater than a? Or maybe divide again by dt (and not multiply). Under my reasoning, the fact that the velocity is multiplied with dt does not only not correct the velocity but also makes the time between frames irrelevant because it takes the magnitudes to the distance space, and not the velocity space.

This question is related to issue gilestrolab/ethoscope#97 in the ethoscope repository. We are trying to find what is the problem and misannotation could be one of the culprits.

Once again, thank you very much!

Best,
Antonio

[antortjim]

On a sidenote, can this a (velocity_correction_coef), which is effectively used as a distance threshold in my reasoning, be interpreted as being in units of width of the ROI? i.e. a coefficient of 0.003 means 0.3% percent of the width of the tube? I come to this conclusion since a is compared with the distance, which is an euclidean distance between two points that have been normalized with the width of the ROI

https://github.com/gilestrolab/ethoscope/blob/0567585158ff3a51ea998faa97a8d60285711163/src/ethoscope/trackers/adaptive_bg_tracker.py#L468
https://github.com/gilestrolab/ethoscope/blob/0567585158ff3a51ea998faa97a8d60285711163/src/ethoscope/trackers/adaptive_bg_tracker.py#L470

Correct me if I am wrong, but that would then mean that different ROIs require different correction coefficients?
Thanks!

[qgeissmann]

Hi @antortjim,
Sorry if this is confusing! After we got the PLoS paper out, I think we had a new version of raspberry pis. Despite trying to account for FPS in the original calculation, we noticed some discrepancies, therefore we patched our calculation by coming up with an empirical correction. This is explained better in section 2.8.3 of my thesis. I would recommend you to make a histogram of your corrected max velocities overall and check it is bimodalish and that the value 1 is after the low-velocity peak (thesis, fig 2.9C).

Yes, the threshold/correction will need to be scaled by ROI size. ROI size is the only stable reference we had since ethoscopes can be used at higher/lower magnifications (so it is hard to convert pixels to actual distances).

I hope it makes some sort of sense!

[antortjim]

Hi @qgeissmann
Thank you very much for your reply.
It makes sense. But since you are indeed multiplying the velocity by the time difference and dividing by the velocity coefficient, your corrected velocities can also be called 'scaled distances'. So the activity scoring is effectively done based on distances (because you observed that the velocity was sensitive to the time difference). I find it surprising since velocity is intuitively the best activity marker.

Regarding the addition of an extra pixel, since the ROI width is around 457 pixels (as shown in the ROI_MAP table), that means that adding 1 pixel to all measurements implies the ROI-width-scaled distance gets increased 1/457 = 0.0021. Since the threshold is 0.003 (default velocity correction coeff.), this means that the activity threshold on the actual activity of the fly is 0.003 - 0.0021 < 0.001, i.e. < 0.1% of the width of the ROI.

We are trying now to compute the distance also offline in R, using the X and Y columns of the ROI_X tables. I wrote a new class in variables.py that admits float data. Using these X, Y data with float precision, the distance can be computed like this:

euclidean_distance <- function(x, y) {
  square_diffs_x <- (x[-1]-x[1:(length(x)-1)])**2 # horizontal side of each triangle squared
  square_diffs_y <- (y[-1]-y[1:(length(y)-1)])**2 # vertical side of each triangle squared
  result <- c(NA, sqrt(square_diffs_x +  square_diffs_y)) # sqrt of the sum of squares
  return(result)
}


dt_raw <- load_ethoscope(metadata,
                              # change reference hour so
                              # it reflects the incubator settings
                              reference_hour = reference_hour
  )

dt_raw[, xy_dist_log10x1000 := 1000*log10(euclidean_distance(x, y))]
dt <- sleep_annotation(dt_raw)

Using this approach, the results should be identical if I decrease the velocity correction coefficient to match 0.003 - one pixel, right?

Please let me know if my conclusions are wrong!

Thank you very much,

Antonio