Accumulative Difference Image Protocol for Particle Tracking in Fluorescence Microscopy Tested in Mouse Lymphonodes

The basic research in cell biology and in medical sciences makes large use of imaging tools mainly based on confocal fluorescence and, more recently, on non-linear excitation microscopy. Substantially the aim is the recognition of selected targets in the image and their tracking in time. We have developed a particle tracking algorithm optimized for low signal/noise images with a minimum set of requirements on the target size and with no a priori knowledge of the type of motion. The image segmentation, based on a combination of size sensitive filters, does not rely on edge detection and is tailored for targets acquired at low resolution as in most of the in-vivo studies. The particle tracking is performed by building, from a stack of Accumulative Difference Images, a single 2D image in which the motion of the whole set of the particles is coded in time by a color level. This algorithm, tested here on solid-lipid nanoparticles diffusing within cells and on lymphocytes diffusing in lymphonodes, appears to be particularly useful for the cellular and the in-vivo microscopy image processing in which few a priori assumption on the type, the extent and the variability of particle motions, can be done.

• numero_t: number of images on which the number of targets segmented (default is 5) is averaged.
• (x0,y0) -(x1,y1): coordinates of the box on which the computation of the average number of segmented targets is done (default is (30,120) -(80,170)); • num_iter: number of values of the threshold on which (default is 20); • soglia: is the starting value of the threshold; • incrementosoglia: increment of the variable soglia.
While running, the routine writes on the prompt screen the actual value of the threshold used. The user can visually choose which one is the right choice. In the results reported in the Ms. we usually preferred to visually count the target in a specific area of interest and check the identified target.
Output: the routine stores in the matrix bw4(256,256,numero,num_iter) the segmented images.
The i-th segmented image obtained with the threshold value, soglia + k*incrementosoglia, is stored in the matrix bw4 (:,:,i,k). This routine opens two figures that reports the first image (colormap=jet) to be segmented and the corresponding image segmentation obtained with the value of the threshold = soglia + incrementosoglia*num_iter/2 (colormap=hot).
The estimate of the threshold can also be gained from the plot of the vector testsoglia(num_iter,2). This vector is plotted in the figure 5 at the end of the run.
• testsoglia(num_iter,2): contains the average number of segmented targets on the first numero images of the stack (testsoglia(:,2)) as a function of the threshold (testsoglia(:,1)). • Bw4(:,:,i,k): segmented images. The i-th segmented image obtained with the threshold value, soglia + k*incrementosoglia, is stored in the matrix bw4(:,:,i,k).      A_draw_binaries.m. To test the output of A_test_thresholds.m the user runs this routine which provides a series of figures related to the last image segmented in the series, that is the image number=numero of the stack. The variable numero, that assign the total number of images to be read and analyzed, is set by the user in the first lines of the A_imagefilters.m.
• Fig.1: Image after the Minimum filter application.
• Fig.2: Image after the threshold filter

Input.
g: total number of images of which the ADI matrix must be built. The default value is g = numero. It can be changed in the first lines of the routine.

Output:
vel (256,256,g): this is the ADI matrix for the segmented targets and is built from the binary image stack bw3.

C_COORD_VECTOR.M
This routines actually creates the matrix in which the trajectories of all the segmented targets are stored.
These are written as (x,y) coordinates and the corresponding passage time, for each target. The targets that are followed in time are those which are segmented every incrementostep images in the original time stack. This option is essential for the pathological cases of loops or intersecting or interwound trajectories. These cases are then built by composing sub-trajectories detected from different starting images chosen every incrementostep.
During this routine the matrix vel is duplicated into vel2, since the vel matrix is destroyed during the tracking procedure.

Input.
lunghezza: this parameter sets a threshold to the minimum length in pixel of the trajectories that are stored. The default value is 20.
incrementostep: The targets which are followed in time are those which are segmented every incrementostep images in the original time stack. If incrementostep >= numero, the tracking algorithm is performed only once through the whole stack. If incrementostep < numero, the tracking algorithm is performed (numero/incrementostep -1) times, starting from one every incrementostep frames through the whole stack.

Output.
The main output of these routines is the matrix "fine". This is the matrix where the software stores all the data that define the trajectories of all the identified targets. The number of segmented targets is stored in the variable contatore1, that is also the screen output of the routine. The fine matrix is stored for the i-th target (1 ≤ i ≤ contatore1) as: • fine(:,1,i): passage times on the (x,y) pixel.
• fine(:,2,i): X coordinate of the pixel at which the target was found at this passage time.
• fine(:,3,i): Y coordinate of the pixel at which the target was found at this passage time.
The number of items in the first dimension of the matrix fine depends on the target (and therefore on the index i), since it is the length of its trajectory. NOTE: It is important to notice that not all the targets found by this routine corresponds to really different targets in the image, since all the targets that are found every incrementostep images are followed in the ADI image. The variable contatore1 stores the number of all the targets that have been segmented on all the numero/incrementostep images, irrespective of the control on the length of the trajectory. This can be a large number. However the size of the third component of the matrix fine is the value contatore, that is typically much smaller than contatore1, due to the constraint imposed on the minimum length of the trajectory (that can be relaxed by the user acting on the variable lunghezza).

Related routines (phase C).
C_Test_Visualize_Vel.m. This routine provides the user with a fast parallel picture of the trajectories of all the targets that have been segmented, whose number is stored in the variable contatore. This image is opened in the window entitled "Parallel Visualization of the trajectories". The color codes, in increasing levels of green, the passage time of the target in the pixels, normalized to the length of the trajectory for each target. The magenta point on each trajectory indicates the starting pixel of a trajectory that origins from the very first image of the stack. Two figures are opened with a different choice of background for better evaluating the whole development of the trajectory. In this figures no number of targets is reported (see phase D for target number encoding).
C_Test_length.m. The length of each target's trajectory (in pixels) is provided to the user by the routine lunghezzaFINE.m. that can be run after having the compilation of the matrix fine. The length of the trajectory for the i-th target is stored, after the run of lunghezzaFINE.m, in the variable contatorepunti(i).

Input:
max_duration: this is the maximum duration of the trajectory that is tested. This value can be equal to the total number of frames (numero).

Output:
The figure reporting the plot of the trajectory duration as a function of the number of the target.
C_ee_FINE.m. This routine computes for each target the end-to-end distance of the trajectory and stores it in the vector contatorepunti2. This is the most valuable parameter for detecting a putative loop in the trajectory. Those trajectories to which corresponds a small value of the end-to-end distance, are likely to be entangled.

Output:
The routine provides a plot of the end-to-end distance as a function of the index of the trajectory/target.

D_BETA_ADD_LABEL.M
D_add_label.m. This routine creates a display of the trajectories (the matrix fine) together with the number of the segmented target superimposed on the image itself. The color codes (from dark to light) the passage time on the specific pixel of the target whose index number is written on the image close to the position of the first appearance of the target in the time stack (i.e. the frame from which the target has been followed along the stack). It opens a new figure window entitled "Time projection of the dADI".
The total number of the displayed target numbers can be limited, with respect to the total contatore value, by setting a minimum value of the trajectory length to be displayed.

Input:
lenght_traj: this is the minimum value of the length of the trajectory for which the corresponding target number is displayed on the figure.
T_F_size: the size of the target numbers marked on the image, default is 10. The user should resize the image in order to read correctly the target numbers in highly dense images.
D_T_F_size: the increment of the size of the target numbers marked on the image. The default value is 0. Additional to this feature, the user should resize the image in order to read correctly the target numbers in highly dense images. The algorithm for the text sizing is: colore=max(fine(1,1,:)); % colore = maximum occupancy time in the trajectory T_F_size = T_F_size_0 + D_T_F_size*fine(1,1,i)/colore; %% fine(1,1,i) represents the i-th occupancy time.

Output:
the figure that reports the vel (or ADI) matrix superimposed on the first binarized image of the stack and with the number of the corresponding target marked on the image itself. The color of the target is chosen between red and blue (RGB) according to the image from which the tracking has started. The algorithm for the color coding of the number of the target on the image is: Therefore the most red numbers on the images corresponds to targets that have been segmented on the first image of the stack and the most blue ones corresponds to those segmented on the last image in the stack from which a tracking has been performed. The ADI image is coded in green levels (passage time from dark to light green).
The routine performs also a search for possible superposition of text on the image. In the case that two indices corresponds to the same pixel position (this may occur for incrementostep<numero), the routine displaces the largest occupancy time in the closest nearest neighbors and changes the color of the displaced text to white.
Please notice that the text of the index is written by MATLAB in such a way that the most significant digit lies on the assigned pixel. However, for highly dense images, the user should zoom the image in order to see correctly placed the indices. See also

Regarding_D_add_label.pdf and
AnimatedZoom_forTextPositioning.gif in the OAM, to this regard.

D_DRAW_TARGET_RMS.M
This routine provides the user with a visualization of the trajectory of the selected target. It is suggested that the user run "D_add_label.m" before this step in order to identify the most interesting targets.
The input parameter to be set in the routine is the order number of the selected target (variable: punto) and the output is a series of images and plots.

Input:
• punto: is the variable in which to store the number of the target of which one wants to display the trajectory.
• High_T: is the value of the threshold to be used to display the images in Fig.3 (see output below. This value enhance the details of the images and is used only for display purposes. Output: • Fig.1: this figure, whose title is "TRACE" displays the trajectory of the selected target (specified by the variable punto) in jet colormap. The color codes for the increasing passage time along the trajectory, normalized to the length of the trajectory for the specified target.
• Fig.2. this figure, whose title is "rms displacement" ( displays the root mean square displacement of the selected target.

E_CHECK_CROSSINGS.M
This routine has been devised to detect the entanglement points in trajectories. The user must select a target and the routine analyzes its trajectory looking for space-time points (x,y,t) of other targets' trajectories that lie "close" to the selected trajectory. This test is performed by selecting a (3D) volume whose size, ∆x∆y∆t, is set by the user in the routine. The output is a table of targets whose trajectories lies within the preset xyt volume with the trajectory of the target chosen to be analyzed.

Input.
Punto_iniziale: this variable contains the number of the target whose trajectory must be scanned for possible interactions with other targets.
delta_t: this is the half of the time axis of the hypersphere. The default value is 3. delta_x: this is the half of the x-axis of the hypersphere. The default value is 2. delta_y: this is the half of the y-axis of the hypersphere. The default value is 2.
The routine searches for the instances that corresponds to trajectories that, at times t 1 and t 2 that lies within 2×delta_t +1, passes for points (x 1 ,y 1 ) and (x 2 ,y 2 ) in the image that lie within a cube of size (2×delta_x+1)×( 2×delta_y+1). The number of couples of such points (not trajectories but points along the trajectories) is stored in the variable controllo_contatore.

Output.
The output of the routine is the matrix controllo_tabella whose size is controllo_contatore ×8. The format is the following.
An example of the controllo_tabella is given hereafter.
Therefore in the above example, the targets number 594, 604, 612, 615, 921, 923 and 940 may have interactions with the target 290.

F_CROSS_ADD_LABEL.M
This routine displays the trajectories of those targets that present some interactions with the target whose index is Punto_iniziale, as selected in phase E.

Input.
PuntoIniziale: this value is the index of the trajectory/target of which we would like to consider possible superpositions with other targets' trajectories. This input is set in the E_check_crossings.m.

Output.
A figure is opened reporting the ADI image on which the indices of the trajectories that may have some superposition with the chosen target's trajectory, as set by the content of the matrix controllo_tabella, are displayed. If the value of PuntoIniziale does not correspond to a target with interactions with other targets only the trajectory of the target PuntoIniziale will be shown on the image.

G_CREATE_MOVIE.M
This routine creates a movie of the trajectory of the selected target and allows to perform simple movie operations, such as zoom, back and forward, pause and frame-by-frame visualization of all the trajectories that are linked to the selected target by the matrix controllo_tabella. This routine must then be run after E_check_crossings.m (that creates the controllo_tabella table) and after setting the proper colormap.
Each of the targets present in the table controllo_tabella, and therefore possibly related to the target with index Punto_iniziale, is reported in the images of the stack as a moving circle with colors that falls within the colormap cmap2. All the circles are reported on the images: they start moving when the corresponding target moves on the image and disappear when the target has been lost on the image. The image is displayed instead as a 8 bits color map according to the colormap cmap1.

Input:
High_T: this is the threshold value used to amplify the contrast on the image.
a: is the actual number of images to be mounted in the movie. The default is a=numero, i.e. all the numero images are mounted.
Max_tab: is the maximum number of targets linked to the target analyzed (whose index is Punto_iniziale) to be displayed on the movie. This number should be larger than the rows in the table controllo_tabella. The default value is Max_tab=30, which more than sufficient to treat all the linked targets. It can be lowered by the user if only some of the linked targets are to be followed on the movie.
cmap and cmap2. These are the color maps for the image and for the circles that describe the position of the targets. The default values are: jet(128) for the image and hsv(128) for the circles. This parameter can be changed by acting on the input parameter lut_flag in the following way: lut_flag = 1; % lut_flag = 0 --> jet, lut_flag = 1 --> hot

Output:
The movie is mounted by copying the frames initially displayed on the screen on an avi file, whose name is created as 'Gmovie.avi'.
The Movie is also opened on the screen by means of the Matlab encoded avi player, implay.

H_TRACK_3D.M
This routine allows the user to superimpose the trajectories that refer to the targets that have been found to have interactions with the analyzed target (whose index is Punto_iniziale). It displays a 3D plot in which the trajectory found for each of the targets listed in the table controllo_tabella, are displayed on different z planes. The z coordinate corresponds to the order with which the target appears in the controllo_tabella. The plot is opened first on the z projection: all the trajectories are superimposed. The user can then use the "rotate 3D" button on the image window menu to visualize different projections.
The starting point of each trajectory is marked by a red sphere. The actual target number that corresponds to a trajectory reported in this plot can be gained by the inspection of the controllo_tabella. If we take as an example the controllo_tabella reported above (see description of phase E), the trajectory reported at z=3 by the H_track_3D.m routine corresponds to the target 612.   In order to build the most complete paragon table for one target followed along the stack and to be fed to the phase K of the algorithm, the user should perform the phases E I for all the targets that on the image appear to have possible interactions, and combine results similar to those shown in table I-2 obtained for different values of Punto_iniziale.

K_PREPARE_PARAGON.M, K_FINAL_MOVIE.M AND K_BETA_FINAL_MOVIE_TRAJ.M
These routines help the user in building up a final trajectory movie and to compute the x,y trajectory to be used for further physical analysis. The K_final_movie.m and K_beta_final_movie_traj.m routines, that differ only for the display mode of the trajectory on the movie, take as an input the table paragon, created by the phase I of the algorithm and mount a movie with the trajectory superimposed on the image stack.
K_final_movie.m reports the trajectory as a collection of circles that follow the target in its motion.
K_final_movie_traj.m reports the trajectory as a trace that is superimposed on all the frames of the final movie.
The user must validate the table paragon before using these routines. The validation is performed by looking at the continuity of the trajectory on the table as described in the phase I. In order to perform this task the user can inspect visually the table paragon, as described in the phase I, or use the K_prepare_paragon.m routine.

K_prepare_paragon.m
The routine looks for the continuity in the x component of adjacent rows of the table paragon according to the parameter threshold_paragon.
The routine builds a matrix trj(max_trj_build, max_len_trj,6) in which at most max_trj_build sub-matrices extracted from the table paragon, according to the continuity test on the x coordinate of the trajectory, are stored sub-sequentially. The maximum length (number of elements) of the subtrajectories is max_len_trj.
Since the visualization of a 3D matrix is not easy in the MATLAB editor, the routine provides, at each run, one of the continuous submatrices in a separate 2D matrix, trj(max_len_trj, 6).
out_trj: is the index number of the sub-trajectory to be stored in the output 2D matrix trj_cur(max_len_trj,6) for the user inspection. Default value is out_trj = 1. This parameter should be changed by the user before running the routine in order to obtain a readable 2D matrix to be used to run K_final_movie.m or K_final_movie_traj.m. One needs only to replace the matrix paragone with one of the trj_cur output matrices obtained by this routine (K_prepare_paragon.m).
max_len_trj. Is the maximum length of the sub-trajectories to be built. Default value is max_len_trj = 200;

Output:
trj_cur. This is a 2D matrix that stores a subset of the table paragon that satisfies the continuity constraint. This table represents the projection of the trj(max_trj_build,max_len_trj,6) on one of the indices (chosen by setting the index out_trj before running the routine) of the first dimension: trj_cur( : , : ) = trj(out_trj, : , : ) The trj_cur matrix has the same dimension(paragone(max_len_trj,6)) as paragone and can therefore be used as it is to replace paragone before running K_final_movie.m or K_final_movie_traj.m. The whole sets of continuous sub-trajectories is stored in trj (max_trj_build, max_len_trj,6). In order to get all the subtrajctories in a format Let us take first a reference image. This is taken from the stack of images reporting the motion of the SLN nanoparticles within the cells, reported in Fig.2 and discussed in the text.
The image reported below corresponds to the bottom row of images in Fig.2 D_beta _add_label.m implements the possibility to choose the text font size. We could also give different sizes of the numbers reported in the image since the size is chosen according to two parameters T_F_size_0 = 10; % initial font size D_T_F_size = 0; % increment on the font size and the following algorithm: T_F_size = T_F_size_0 + D_T_F_size*fine(1,1,i)/colore; where fine(1,1,i) is the occupancy time of the pixel and the variable colore reports the maximum value of the occupancy time for the trajectory.

Color map of the trajectory.
Let us consider first the issues related to the choice of the colors of the trajectory pixels. This is a green colormap as reported in the Figures. However we also encode the fact that a trajectory starts from the very first frames of the stack, by coloring magenta that pixel.
Take the magenta pixel numbered "18". This corresponds to the 18th trajectory. It corresponds to the pixel (215,65) with RGB values (output of the D_add_label.m) (R,G,B) = (0.975,0.0167,0.975). By debugging the run we get that this pixel corresponds to the starting point of a trajectory very long (it is mapped on 116 out of 120 frames total) and starts from the very first frame.
The algorithm for building up the colors in the trajectory is implemented in D_add_label.m as follows: for k=1: 256 for j=1: 256 if abs(bw3(k,j,1)-bw3(k,j,i+1))>0.5, Htempo3(k,j,1)=Htempo3(k,j,1)+1; if Htempo3(k,j,2)==0, Htempo3(k,j,2)=Htempo3(k,j,2)+i;end; Htempo3(k,j,3)=Htempo3(k,j,3)+1; end end end end; The Htempo3(x,y,:) matrix contains the RGB components of the (x,y) pixel. This matrix contains, at the end of the run, a picture of the motion. The piece of code reported above, instructs the program to assign to the G component (i.e. Htempo(x,y,2)) of the (x,y) pixel the value of that frame ("i" in the code) on which some motion has been detected for the first time on the specified pixel. The R and B components (Htempo(x,y,1) and Htempo(x,y,3)) are continuously incremented each time that the pixel has been interested by a motion. So, at the end, these two values report the number of times we have detected a motion on the specified pixel. Therefore those long trajectories that starts from the first frames corresponds to RGB triples in which the G value is much less than the R and B ones, and result in a magenta color.
Take the dark green pixel numbered "337". This corresponds to the 337th trajectory. It corresponds to the pixel (225,78) with RGB values (output of the D_add_label.m) (R,G,B) = (0.0417,0.392,0.0417). By debugging the run we get that on this pixel a change occurred for the first time on the 47 th frame and that only on 5 frames total a change occurred on this pixel. Therefore the Htempo3(225,78,2)= 47 and Htempo3(225,78,1)=Htempo3(225,78,3) =5. In such a case the R and B values are small (0.0417) and the G value is relatively large (0.392) but definitely smaller than one.
Take the pixel (219,78). This pixel does not correspond to the starting position of any trajectory: its position does not correspond to (x,y) = (fine(1,3,i), fine(1,2,i)) for any value of i=1..punto_tot (punto_tot is the total number of trajectories). The content of the image on this pixel corresponds to RGB values (output of the D_add_label.m) (R,G,B) = (0.00833,0.95,0.00833).
By debugging the run we get that on this pixel a change occurred for the first time on the 114 th frame (out of 120 analyzed) and that only on 1 frame a change occurred on this pixel. Therefore the Htempo3(225,78,2)= 114 and Htempo3(225,78,1)= Htempo3(225,78,3) = 1. In such a case the R and B values are small (0.00833) and the G value is very large (0.95), so that the pixel is bright.

Positioning of the text on the image.
Regarding the position of the target number on the output image of D_add_label.m we note that the number is always set in such a way that the most significant digit lies at the center of the pixel and the other digits are then found on the right side of the pixel. As an example of this algorithm we report the following image: A specific parameter (T_F_size=12) in the D_add_label routine specifies for the font size. It may help to lower the value of this parameter to avoid the partial superposition of targets' number on the image.
For example if we run D_add_label with length_traj=2; T_F_size_0 = 12; We obtain the following (zoomed) image on which some of the targets' number are superimposed. This is due to the large font size and the high density of targets due to the low constraints imposed by the algorithm (particularly low value of the threshold on the length of the trajectory to be displayed, length_traj=2).
If on the contrary we run the D_add_label with a lower font size, We obtain for the same area of the image the following result in which no superposition of targets' indices can be found. For example if we run D_add_label with length_traj=2; T_F_size_0 = 8; In addition to the choice of the font size we have now added to the D_add_label.m routine a feature that check for possible real superposition of the numbers on the image. This may occur for entangled trajectories, since the same pixel may belong to more than one trajectory and may be tracked from different initial frames (if incrementostep is chosen less than the total number of frames). If a superposition occurs the routine displays by one unit in the nearest neighbors the text reporting the target number and uses for this text the white color, in order to notify the user of this instance.
It is very important to notice that Matlab does not scale the size of the text on an image when the user zoom in or out the image itself. The location of the numbers with respect to starting pixel of the trajectory depends then on the zoom of the image. However, by zooming in the image the location of the numbers with respect to the starting pixels becomes more and more accurate, as can be gained by the following sequence of images taken at increasing zoom level (see also the animated gif image: AnimatedZoom_forTextPositioning.gif as OAM).