In rogerswt/flowMie: Modeling a Detector Response to Scattering of (Multi-layer) Spherical Particles

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
# pre-load libraries
library(flowCore)
library(Rscattnlay)
library(flowMie)
library(wadeTools)
library(fields)
library(KernSmooth)

Introduction

There has been a great deal of interest in recent years in Extracellular Vesicles (EVs), including exosomes and microvesicles. There are a number of methods of detecting and characterizing EVs. Flow cytometry is capable of multi-parameter EV detection and characterization, but not all instruments are capable of this difficult feat, and of those that are, none of them are capable of seeing the smallest EVs. That is, flow cytometry has yet to plumb the small-size depths of EV size distributions. Nevertheless, flow cytometry is an extremely important platform for EV characterization. Thus, it's important to quantitatively understand the capabilities and limitations of the instruments that are used for these applications.

flowMie provides tools to describe small (single- or multi-layer) particles, including calibration beads as well as EVs, along with an optical description of the flow cytometer laser and detector systems, that enable a complete and accurate modeling environment to calculate and characterize the response of your flow cytometer to these very small and dim particles. We also implement the Mie Transform which uses calibrated flow cytometry data to directly compute the diameter of each event in an FCS file.

Describe the Flow Cytometer

The first step is to describe the physical characteristics of the flow cytometer that affect its response to the scattering of light. These characteristics include:

The nominal scattering angle of the detector
The scattering detector half-angle of acceptance with respect to the objective's focal plane
Possible variation of efficiency across the face of the objective
The laser polarization

Please note that as of now flowMie assumes the detector aperture to be a round disk. Oddly, the laser wavelength isn't modeled here, in the detector, but rather in the particle description (see below, due to the organization of objects in Rscattnlay).

# create a detector, 
#    - located at 90 degrees with respect to the incident laser light,
#    - with a detector half-angle of 60 degrees,
#    - laser is linearly polarized perpendicular to the plane of incidence
detector = create_detector(theta_0 = 90, alpha = 60, psi_0 = 90, pol = 1.0, eta = eta_uniform)

Describe Particles

The function create_particle() flexibly describes a particle characterized by:

How many layers it has
The radius of each layer
the refractive index of each layer

as well as:

the refractive index of the medium in which it's suspended
the wavelength of the incident light

# create a 2-layer particle
#    - the first layer has a radius of 95 nm and RI of 1.38 (the core)
#    - the second layer has a radius of 100 nm and RI of 1.46 (the membrane)
# The particle is in a fluid whose RI is 1.34, corresponding to normal saline
# The incident wavelenth is 488 (the blue laser, off which the SSC detector is located)
vesicle = create_particle(medium = 1.34, lambda = 488, n_layers = 2, r = c(95, 100), n = c(1.38, 1.46))

In addition to this general-purpose function there are a few convenience functions with built-in refractive indices, and which accept the size description as a diameter rather than a radius (in keeping with the way calibration beads are typically specified by their manufacturers). These functions are simply wrappers for create_particle():

# a 100 nm diameter polystyrene bead in saline, illuminated by the blue laser
bead_polystyrene = create_PS(d = 100)

# a 200 nm diameter silica bead in saline, illuminated by the blue laser
bead_silica = create_SI(d = 200)

# an EV in saline, illuminated by the blue laser
ev = create_EV(d = 200)

The best information we have found on refractive indices of various materials were obtained from the following:

EV cytosol and membrane indices, obtained from fluorescence lifetime measurements:
- Henk-Jan van Manen et. al., Biophys J. 2008 Apr 15;94(8):L67-9. doi: 10.1529/biophysj.107.127837. Epub 2008 Jan 25.
Silica:
- https://refractiveindex.info/?shelf=glass&book=fused_silica&page=Malitson
Polystyrene
- https://refractiveindex.info/?shelf=organic&book=polystyren&page=Sultanova
Gold
- https://refractiveindex.info/?shelf=main&book=Au&page=Johnson

Indices included as default values are all for 488 nm radiation. If you use a different laser (e.g. 405 nm) for light scatter measurements please override the default RI values n in the appropriate functions.

Calculate Scattering Response

Now that we've described the instrument and the particle we can perform scattering calculations.

# how big is the scattering signal with our default detector due to a 200 nm diameter silica bead?
val = calculate_detector_response(particle = create_SI(d = 200), detector = create_detector())
val

Let's look at the size-dependence of several types of particles, with a realistic instrument, and with a hypothetical instrument with a very narrow acceptance half-angle:

# make a hypothetical instrument with narrow acceptance angle for SSC
instr_narrow = create_detector(alpha = 1)

# make an instrument that corresponds to a BD FACSCanto
instr_canto = create_detector()

dia = seq(20, 1000, by = 5)
ssc_au   = ssc_ps   = ssc_si   = ssc_ev   = vector('numeric')
ssc_au_n = ssc_ps_n = ssc_si_n = ssc_ev_n = vector('numeric')
for (i in 1:length(dia)) {
  # gold nanoparticles.  Note the complex refractive index, characteristic of metals
  gold = create_particle(n_layers = 1, r = dia[i] / 2, n = 1.1271+1.8382i)
  ssc_au[i]   = calculate_detector_response(gold, detector = instr_canto)
  ssc_au_n[i] = calculate_detector_response(gold, detector = instr_narrow)

  # polystyrene
  ps = create_PS(d = dia[i])
  ssc_ps[i]   = calculate_detector_response(ps, detector = instr_canto)
  ssc_ps_n[i] = calculate_detector_response(ps, detector = instr_narrow)

  # silica
  si = create_SI(d = dia[i])
  ssc_si[i]   = calculate_detector_response(si, detector = instr_canto)
  ssc_si_n[i] = calculate_detector_response(si, detector = instr_narrow)

  # EV
  ev = create_EV(d = dia[i])
  ssc_ev[i]   = calculate_detector_response(ev, detector = instr_canto)
  ssc_ev_n[i] = calculate_detector_response(ev, detector = instr_narrow) 
}

# plot the results
par(mar = c(5, 5, 4, 1))
decades = c(1e-8, 1e1)
cex = 1.5
lwd = 2
plot(dia, ssc_au, type = 'l', log = 'y', ylim = decades, xlim = c(0, 1000),
     xlab = "Particle Diameter (nm)", ylab = "Scattering Intensity (Arb. Units)",
     yaxt = 'n', cex.axis = cex, cex.lab = cex, lwd = lwd)
axis(side = 2, at = 10^(-8:1), labels = tick.labels(-8, 1), cex.axis = cex)
lines(dia, ssc_au_n, col = 'black', lty = 'dotdash', lwd = lwd)

lines(dia, ssc_ps, col = 'indianred2', lwd = lwd)
lines(dia, ssc_ps_n, col = 'indianred2', lty = 'dotdash', lwd = lwd)

lines(dia, ssc_si, col = 'dodgerblue2', lwd = lwd)
lines(dia, ssc_si_n, col = 'dodgerblue2', lty = 'dotdash', lwd = lwd)

lines(dia, ssc_ev, col = 'springgreen4', lwd = lwd)
lines(dia, ssc_ev_n, col = 'springgreen4', lty = 'dotdash', lwd = lwd)

comp_dia = 100
idx = which(dia == comp_dia)
text(x = comp_dia, y = ssc_au[idx], labels = "Au", col = 'black', cex = cex, pos = 4)
text(x = comp_dia, y = ssc_ps[idx], labels = "PS", col = 'indianred2', cex = cex, pos = 4)
text(x = comp_dia, y = ssc_si[idx], labels = "SI", col = 'dodgerblue2', cex = cex, pos = 4)
text(x = comp_dia, y = ssc_ev[idx], labels = "EV", col = 'springgreen4', cex = cex, pos = 4)

ssc_au[idx] / ssc_ps[idx]
ssc_ps[idx] / ssc_si[idx]
ssc_si[idx] / ssc_ev[idx]
ssc_au[idx] / ssc_ev[idx]

There are a few things to note about this picture.

First and foremost, the y-axis is logarithmic and covers 10 decades of dynamic range (you just don't see this very often!).
Once the particle size drops below about 200-300 nm (or about half the wavelength of the incident light, which in this example is 488 nm), the curves become ridiculously steep. This means that a small reduction in the particle size results in a large decrease in the signal it will generate. This sets a pretty hard stop on the size limit of detectability in most cytometers.
The material of which the particle is made has a profound effect on the amount of light it scatters. Gold particles (black curves) >> polystyrene (red curves) >> silica (blue curves) >> EVs (green curves). Comparing at 100 nm diameter,
- gold scatters almost 50 times as much as polystyrene
- polystyrene scatters about 5 times as much as silica
- silica scatters about 4 times as much as an EV
- gold scatters almost 1000 times as much as an EV
This dependence on the particle material is why the choice of reference materials is so important. Artificial particles whose size and refractive index are most similar to EVs make the most appropriate reference materials.
Comparing the solid versus the dot-dashed curves, we see the smoothing effect due to the large acceptance angle of real-life flow cytometers. They are the way they are because high numerical aperture optics are generally part of the cytometer design in order to collect as much light as possible to optimize sensitivity. This means that the side scatter detector doesn't just measure 90 deg. scattering, but a range of +- 60 or so degrees around that 90 deg nominal angle. The wiggles, or "Mie Resonances", are very pronounced with a narrow detector, and not so much with a realistically described flow cytometer.
For instruments with large numerical aperture, the Mie scattering signal is pretty much monotonic up to about 500-600 nm. This is important for computing the Mie Transform (see the next section). If the scattering signal is NOT monotonic, then a given scattering signal could correspond to two (or more) vastly different particle sizes (draw a horizontal line through the Mie Resonances in the above figure to see what I mean). For example, EV Mie Resonances intersect the scattering intensity of 10^-3^ at the following EV sizes: 135, 315, 370, 585, 615 and 850 nm. Fortunately, for particle sizes of interest (<500 nm typically) and large numerical aperture, this isn't the case.

Polarization Effects

All other treatments of Mie Theory for flow cytometry of which we are aware have not taken into account the polarization state of the incident light. This is a small but significant effect. Given that most flow cytometers use lasers as the light source, and since lasers nearly always produce linearly polarized light, it seems to us that this should be taken into account.

det_unpol = create_detector(theta_0 = 90, alpha = 60, psi_0 = 90, pol = 0)
det_pol   = create_detector(theta_0 = 90, alpha = 60, psi_0 = 90, pol = 1.0)
det_pol_in   = create_detector(theta_0 = 90, alpha = 60, psi_0 = 0, pol = 1.0)

dia = seq(20, 1000, by = 5)
ssc_si    = vector('numeric')
ssc_si_un = vector('numeric')
ssc_si_in = vector('numeric')
for (i in 1:length(dia)) {
  si = create_SI(d = dia[i])
  ssc_si[i]   = calculate_detector_response(si, detector = det_pol)
  ssc_si_un[i] = calculate_detector_response(si, detector = det_unpol)
  ssc_si_in[i] = calculate_detector_response(si, detector = det_pol_in)
}

# plot the results
par(mar = c(5, 5, 4, 1))
decades = c(1e-3, 1e-1)
cex = 1.5
lwd = 2
plot(dia, ssc_si, type = 'l', lty = 'dotdash', log = 'y', ylim = decades, xlim = c(100, 300),
     xlab = "Particle Diameter (nm)", ylab = "Scattering Intensity (Arb. Units)",
     yaxt = 'n', cex.axis = cex, cex.lab = cex, lwd = lwd)
axis(side = 2, at = 10^(-8:1), labels = tick.labels(-8, 1), cex.axis = cex)
lines(dia, ssc_si_un, col = 'black', lwd = lwd)
lines(dia, ssc_si_in, col = 'black', lty = 'dotted', lwd = lwd)
x.text = 200
y.text = 5e-3
text(x = x.text, y = y.text, labels = "Solid: Unpolarized", cex = 2, pos = 4)
text(x = x.text, y = y.text / 1.5, labels = "DotDash: Perpendicular Polarized", cex = 2, pos = 4)
text(x = x.text, y = y.text / 1.5^2, labels = "Dotted: Parallel Polarized", cex = 2, pos = 4)

mean(ssc_si / ssc_si_in)
max(ssc_si / ssc_si_in)
min(ssc_si / ssc_si_in)

The differences, while relatively small, are not negligible, and can lead to a systematic errors in the estimation of EV sizes. Differences between polarization in-plane versus perpendicular to the plane range from about 12% to 32%, with an average about 20%. Polarized light perpendicular to the plane of incidence yields higher scattering, whereas polarization in the plane of incidence produces lower scattering. Unpolarized light is somewhere in between.

Polarization effects are a function of the detector numerical aperture. For lower numerical aperture instruments, polarization effects can lead to differences of several-fold, as opposed to 10 - 30% for high numerical aperture machines. Thus, care should be taken to accurately characterize the detector, in terms of both detector numerical aperture (or equivalently the half-angle of acceptance) as well as the polarization state of the illuminating laser light.

We welcome instrument manufacturers to submit specifications for inclusion in this package. Please send an email to Wade Rogers to start a discussion.

The Mie Transform

We often look at EVs in a flow cytometry dataset by plotting a fluorescence parameter against, say, SSC-A. However, we now have the means, via the Mie Transform, to directly invert the signal in the data to EV diameter. The function mie_transform() performs this inversion. We will use this function to add a physical size parameter to a flowFrame so that particle size distributions can be directly studied.

library(flowCore)    # core functions for flow cytometry files
library(wadeTools)   # for some convenient functions, like doTransform(), bx(), ibx(), etc.
library(fields)      # xline, yline
library(KernSmooth)  # bkde

# load bead data.  The beads are assumed to be 200 nm Polystyrene
fn_bead = system.file("extdata", "beads.fcs", package = "flowMie")
ff_bead = read.FCS(fn_bead)

# load triton data.  We'll use this to gate CD41a+ EVs from background
fn_triton = system.file("extdata", "sampled_triton.fcs", package = "flowMie")
ff_triton = read.FCS(fn_triton)

# load the EV data.
fn_ev = system.file("extdata", "sampled_ev.fcs", package = "flowMie")
ff_ev = read.FCS(fn_ev)

# biexp transform the data
ff_bead = doTransform(ff_bead, cols = 1:4, method = 'biexp')
ff_triton = doTransform(ff_triton, cols = 1:4, method = 'biexp')
ff_ev = doTransform(ff_ev, cols = 1:4, method = 'biexp')

# find the bead MFI on SSC-H
kde_bead = bkde(exprs(ff_bead)[, "SSC-H"], bandwidth = 0.05)
bead_mfi = ibx(find.local.maxima(kde_bead)$x)

# calibrate the gain factor.  We'll use this later for mie_transform().
detector = calibrate_ssc_conversion(detector = create_detector(), 
                                     bead = create_PS(d = 200), 
                                     bead_mfi = bead_mfi)

# Compute the Mie Transform and add a size parameter to the EV data
# Note that we specify transformation = "biexp" to indicate that the bead data
# were biexponentially transformed.
ffs = mie_transform(ff = ff_ev, ssc_param = "SSC-H", transformation = 'biexp', detector = detector)

# gate the EV data based on the triton negative control (for fluorescence) and the calibration beads (for size)
kde_triton = bkde(exprs(ff_triton)[, "780/60 Red-A"], bandwidth = 0.05)
kde_triton$y = bx(1e5) * kde_triton$y / max(kde_triton$y)
ckde = cumm.kde(kde_triton)
thresh = kde_triton$x[min(which(ckde$y > 0.995))]
ffs_gated = Subset(ffs, rectangleGate("780/60 Red-A" = c(thresh, Inf),
                                      "SSC-H" = c(-Inf, bx(bead_mfi)))
)
kde_size = bkde(exprs(ffs_gated)[, "Size (nm)"], bandwidth = 5.0)
kde_size$y = bx(1e5) * kde_size$y / max(kde_size$y)

# Visualize results
par(mfrow = c(2, 2))

# first plot: the beads
pplot(ff_bead, c("780/60 Red-A", "SSC-H"), main = "Calibration Beads")
lines(kde_bead$y, kde_bead$x)
yline(bx(bead_mfi), col = 'red', lty = 'dotdash', lwd = 2)

# second plot: the triton control and fluorescence threshold
pplot(ff_triton, c("780/60 Red-A", "SSC-H"), main = "Triton Negative Control")
lines(kde_triton)
xline(thresh, col = 'red', lty = 'dotdash', lwd = 2)

# third plot: the raw EV data
pplot(ffs, c("780/60 Red-A", "SSC-H"), main = "EV Data")
xline(thresh, col = 'red', lty = 'dotdash', lwd = 2)
yline(bx(bead_mfi), col = 'red', lty = 'dotdash', lwd = 2)
kde_size$y = bx(1e5) * kde_size$y / max(kde_size$y)

# fourth plot: gated EV size distribution
pplot(ffs_gated, c("780/60 Red-A", "Size (nm)"), ty = 'linear', main = "Gated EV Data by Size")
points(exprs(ffs_gated)[, c("780/60 Red-A", "Size (nm)")], pch = 20, cex = .2)
lines(kde_size$y, kde_size$x)
x.text = bx(1e3)
y.text = 400
text(x = x.text, y = y.text, 
     labels = sprintf("Median EV size = %.1f nm", median(exprs(ffs_gated)[, "Size (nm)"])),
     pos = 4, cex = 2)
text(x = x.text, y = y.text - 30, 
     labels = sprintf("Minimum EV size = %.1f nm", min(exprs(ffs_gated)[, "Size (nm)"])), 
     pos = 4, cex = 2)
text(x = x.text, y = y.text - 60, 
     labels = sprintf("Stand Dev EV size = %.1f nm", sd(exprs(ffs_gated)[, "Size (nm)"])),
     pos = 4, cex = 2)

The upper left panel shows the analysis of the calibration beads.
The upper right panel deduces a fluorescence threshold based on a tritonized sample, in which EVs are wiped out but sources of background remain.
The lower left panel shows the EV sample, analyzed under identical conditions as both the bead and triton runs. The red lines indicate cut-offs used to gate the EV region.
Finally, the lower right panel shows the gated EVs, with the Y-axis indicating particle size in nm as computed with the Mie Transform. We note that the sharp rise of the distribution at about 150 nm is an instrument effect, and not to the actual underlying EV size distribution. In this case the instrument was triggering on SSC, and the trigger threshold was high enough to reduce excess event rates, thereby limiting the lower size limit of this run. Triggering on a fluorescence signal would presumably lower the size limit of detection.

Concluding Remarks

I hope you find this package useful in understanding your flow vesiculometry data.

This package is in alpha development. I encourage you to contact me if you have any comments, bug reports, or suggestions for improvements and extensions.

Wade Rogers

Still Pond Cytomics

Please cite this package if you use it in a publication.

rogerswt/flowMie documentation built on Jan. 31, 2021, 8 a.m.

rdrr.io home R language documentation Run R code online

CRAN packages Bioconductor packages R-Forge packages GitHub packages

Note that we can't provide technical support on individual packages. You should contact the package authors for that.

Tweet to @rdrrHQ

GitHub issue tracker

ian@mutexlabs.com