Improving Multicolor Colocalization in Single-Vesicle Flow Cytometry with Vesicle Transit Time

Immunophenotyping of vesicles, such as extracellular vesicles (EVs), is essential to understanding their origin and biological role. We previously described a custom-built flow analyzer that utilizes a gravity-driven flow, high numerical aperture objective, and micrometer-sized flow channels to reach the sensitivity needed for fast multidimensional analysis of the surface proteins of EVs, even down to the smallest EVs (e.g., ∼30–40 nm). It is difficult to flow focus small EVs, and thus, the transiting EVs exhibit a distribution in particle velocities due to the laminar flow. This distribution of vesicle velocities leads to potentially incorrect results when immunophenotyping nanometer-sized vesicles using cross-correlation analysis (Xcorr), as the order of appearance of the vesicles might not be the same at different spatially offset laser excitation regions. Here, we describe an alternative cross-correlation analysis strategy (Scorr), which uses information on particle transit time across the laser excitation beam width to improve multicolor colocalization in single-vesicle immunoprofiling. We tested the performance of the algorithm for colocalization analysis of multicolor nanobeads and EVs experimentally and via simulations and found that Scorr improved both the efficiency and accuracy of colocalization versus Xcorr. As shown from Monte Carlo simulations, Scorr provided an ∼1.2–4.7-fold increase in the number of colocalized peaks and ensured negligible colocalization of peaks. In silico results were in good agreement with experimental data, which showed an increase in colocalized peaks of ∼1.3–2.5-fold and ∼1.2–2-fold for multicolor beads and EVs, respectively.


1) Monte Carlo simulation
The flow trajectory of channel #1 was generated as follows: i) From experimental data on multi-color beads, we calculcated the distribution of distances between consecutive peaks at different particles densities (see Fig. S1a). We used these distributions as weights during random selection of peak occurrences in the blue trajectory. This allowed as to spread a pre-defined number of peaks (580, 3400 or 30200, for the low, medium or high density scenario, respectively) along the flow trajectory in a way more representative for the real scanarios (see Fig. S1b-c).
ii) Locations of brighter peaks (either 5% or 50% of total peak, with 3-fold higher intensity) were selected randomly, to resemble the real case of particles of different sizes flowing within the microfluidic channel.
iii) From experimental data on multi-color beads, we derived the spatial distribution of particles within the 2x2 μm channel, by rearranging Poiseuille's equation (see Eq. S1).
Eq. S1 R refers to the channel radius in µm, PTmin refers to the smallest passing time of particles (i.e., fastest particle recorded) whereas PTi refers to the passing time of the i-th particle. Our results (see Fig. 1e) were in agreement with the literature, 1 which shows a focusing of particles at ~ 0.6 times the channel size due to the balancing between two forces: a shear gradient (from the centre of the channel to the wall) and a lift force (wall effect) from the wall to the centre. Thus, we used the experimental distribution of particles' location as weights to randomly assign a a given travelling speed to individual peaks. This ensured similar PT distributions between experimental and simulation data (see Fig S1f- The flow trajectory of channel #2 was generated as follows: i) each peak from channel #1 was shifted forward in time according to Eq. S5, S8 (see below), using the preseleceted parameters for laser widths and laser-to-laser distance (see main text) and the pre-assigned peakspecific travelling speed (i.e., transit time).
ii) To each peak-specific travelling speed (PT), a random ±10% error was added, to resemble the real case of particles flowing within our µm-sized channel (see below). 1

2) Algorithm workflow for Scorr colocalization
Let's consider the simple case of particles flowing in a 2-color system, where laser #2 is placed after laser #1, along the flow direction. Therefore, two distinct flow trajectories (T1 and T2) with same number of peaks (N) but different locations in time (p1 or p2, with peaks in T2 being delayed in time compared to their respective peaks in T1) will be measured. To achieve single-vesicle colocalization, the algorithm performs a 4-step process (see Fig. S2a): (i) a "coarse" shift of peaks in T1 forward in time which generates new peak locations (p'1), to identify initial colocalized particles according to a user-defined time window (w1). Specifically, p'1 refers to the theoretical location which a peak in trajectory #1 would have in trajectory #2 according to its travelling speed.
(ii) a "fine" shift of peaks in T1 forward in time which generates a new set of "refined" peak locations (p'1), to re-colocalized particles according to a user-defined time window (w2).
(iii) removal of any residual double colocalization of peaks in T1 with peaks in T2.
(iv) removal of any residual double colocalization of peaks in T2 with peaks in T1.
In each stage, the algorithm generates an NxN matrix where peaks in T1 and T2 run along rows or columns, respectively. After coarse shift, elements in the matrix (Matrix0, see Fig. S2b) are either 0 or 1, depending on whether the two peaks are closer (1) or farther (0) in time than w1. Matrix1 (see Fig. S2c) is generated after fine shift and is built in the same way as Matrix0 but considering w2, instead. Then, 1s are replaced by the Figure S2. a, Schematic of the algorithm workflow. b-e, matrices generated during single vesicle-based colocalization.
corresponding p2-p'1 values (see Fig. S2d), and the algorithm proceeds with removing double colocalization on both channels (see Fig. S2e). In case of multiple colocalization, the algorithm will consider as "true" colocalized those peaks which show the smallest p2-p'1 value (i.e., which are closer in time), or (for same p2p'1 values) those which show the smallest change in travelling speed between T1 and T2 (i.e., particles that are flowing in the same laminar sheet).

3) Derivation of equations for coarse and fine shifting
ΔX1, ΔX2 and ΔX3 represent the blue (488nm) laser width, the red (640nm) laser width and the distance between lasers, whereas PT1, PT2 and PT3 the times required to cover those respective lengths (see Fig. S3a). The travelling speed for each j-th particle is calculated as vi = ΔXi / PTi where i refers to the specific time-distance pair selected (i.e., i = 1 refers to the beam width and PT for the blue laser etc.).
Coarse shift: relies on a first assumption (HP1) that each i-th vesicle-travelling along the microfluidic channel at its own speed (vi)-does not change its velocity (for instance, due to lateral diffusion). Therefore, according to HP1, vi remains constant throughout the channel length and we can write: where 640 is the time at which the n-th particle crosses the red laser and 488 ′ is the theoretical time at which the same particle should cross the blue lasers, if travelling with constant speed v (i.e., if HP1 holds true). Solving Eq. S3 for PT3 and combining it with Eq. S4 we obtain an equation (see Eq. S5) which allows to calculate 488 ′ for each individual vesicle (see Fig. S3b). Eq. S5 Fine shift: Instead, if we consider that the travelling speed of each vesicle might slightly change while crossing the two lasers, we can assume (HP2) that the velocity calculated from the different occurrence in time between peaks in the two trajectories (PT3) and the laser-to-laser distance (ΔX3) is an average of the two velocities estimated from the crossing time at the two lasers (PT1 and PT2). Therefore, we can rewrite Eq. S2 and S3 as following: Eq. S6 Eq. S8

4) Estimation of parameters in Eq. S1 and S3-S8 and vesicle transit time across the laser beam width
The parameters used in Eq. S4-S7 (i.e., the laser beam widths and the pairwise distances) where estimate from a picture of the three lasers used in our setup, which was taken by using a CCD camera, and by performing line profile in ImageJ. The transit times of vesicles (i.e., PTi) across the laser beam width were estimated as previously reported. 2 Briefly, each peak was fit to a gaussian-like curve, and PTi was calculated as six times the peak width (σ). Instead, the error on PTi was calculated from the confidence intervals, namely six times the 95% confidence intervals of σ, obtained from the fitting.

5) DeltaT distributions for low-and medium-density simulated beads
Similar to the case of high-density multicolor beads (see main text), colocalization via Scorr function greatly improved the accuracy of colocalized peaks (i.e., narrower deltaT distributions). For both low-and mediumdensity beads, Scorr yielded ~61% and 89% of peaks colocalized within ± 0.25 ms and ± 0.5 ms, respectively.
Xcorr led similar results only for low density particles, whereas values for colocalized peaks dropped to only ~21% and 37 % within the same range, for medium density particles. Figure S4. a-b, Distribution of deltaT |t'488 -t488| of colocalized peaks with t'488 calculated according to Eq. 3-4 (see main text), obtained from low (a-b) or medium (c-d) simulated data (1:1 dim-bright mixture), setting W1 equal to either 2 ms (wide bars) or 5 ms (thin bars) and W2 = 4 ms. DeltaT are obtained after shifting using Xcorr (a-c) or Scorr, (b-d).
Similar to the scenario of 50% fraction of bright-over-dim peaks, Scorr ensured the false colocalizations remained low (~1%) even at high-density particles. Instead, Xcorr failed to prevent mis-colocalization, yielding up to ~5-fold higher false colocalization for the case of high-density particles (dark red curve). Figure S5. Percentage of mis-colocalized simulated peaks (i.e., dim peaks with brighter ones) as a function of W2 (W1 = 2 ms), at different particles densities (see figure legend) and with 5% fraction of bright-over-dim peaks. Each line refers to shifting via either Xcorr (light-to-dark red lines) or Scorr (light-to-dark blue lines).

Figure S6.
Percentage of true colocalized peaks as a function of W2, obtained from multicolour beads analyzed at either medium-(a-c) or high-density (d-f). Plots refer to different pairwise colocalization: blue-red channels (a, d), blue-green channels (b, e) and green-red channels (c, f). Figure S7. Distribution of peak intensities for multicolor beads at low (green curve) and high (blue curve) particles density. Solid lines represent the respective lognormal fitting.