# SoftwareVignettes

# Biomedical Imaging Group Software Overview

Over the last three decades, the BIG has written and validated a large collection of software tools
(hundreds of programs) for image acquisition, processing, analysis and visualization, image restoration,
computer vision, and the modeling and simulation of biological systems. The bulk of this software is written in
C, as well as in Fortran, C++, Perl, Java and Matlab. These programs range from the simple (e.g., add/subtract/multiply/divide
images; perform global thresholding, etc.) to the complex (e.g., deformable models to identify the plasma
membrane of a cell). The following is a description of many of the programs, roughly sorted by category.
Some of the larger or more important programs are listed at the end with a short description and a link to
a separate web page which provides further information.
*Note: there is additional BIG software described here.*
This software was written, primarily, by BIG members Kevin Fogarty
Lawrence Lifshitz (this includes an overview of Dr. Lifshitz's research areas and a full bibliography),
and Karl Bellve.

## Generating Random Data/Simulations

- CellSim : Simulates chemical reactions and diffusion in 3D cells.
- connected_random_pts : Generates random points in a planar region, calculates how many connected components are created if a point is connected to all other points within a given radius.
- cross_sections_of_voronoi : Creates a random 3D voronoi tessellation and 2D image (cross-section) from it.
- ellipse : Creates an elliptical cell. It is used to model a squashed cell, and see how different it looks after blurring from a cylindrical cell.
- makebeads2 : Randomly distributes beads within a volume; created to simulate the distribution of glut4 labeled vesicles in the TIRF field of a microscope.
- oxygen : Calculates steady state partial oxygen pressure (P) as a function of distance from a blood vessel (r).
- random_paths.pl: Generates a random path (ie, diffusion of a single particle).
- vcspat : Calculates complete spatial randomness (Ripley's K function) of punctates in a volume.
- vcspat-models: Drops random objects in a volume. Objects can either be individual voxels, tiny images, or connected voxels in an image.

## Statistics

- coloc3way image1[.i2i] image2[.i2i] image3[.i2i] [coloc-image[.i2i]] : Calculates two-way and three-way binary image colocation.
- overlap3: Produces a whole set of numbers with colocalization for each of many offset directions of one image relative to the other (to compensate for slight mis-registration between the images).

## Finding events, Tracking

- exo6 : Identifies vesicle fusion events in a 2D time series. Intensity maxima are tracked, then each is analyzed to see if it meets the criteria to be a fusing vesicle.
- track : Tracks regions over time. Tracks written to stdout in rpts format. This is done by thresholding the input image. Resulting connected voxels define regions. Uses one of 3 tracking algorithms, lots of parameters can be specified to customize the tracking algorithm.
- trackinfo.pl : Analyzes output from track.
- track_pts : This programs takes a list of initial points in a 3D image and tracks them over time. It uses hill climbing or the brightest point in the neighborhood to link to in the next time point.
- paths4.pl : Reads tracking information from *.txt files, generates statistics for each track (speed, pathlength, distance, curvature). Statistics printed to stdout. There are lots of filtering and analysis options.

## Miscellaneous

- parse_macros.pl (and associated .pl files) : Produces layered ImageJ macro files with variables passed through between them. This makes it much easier to write ImageJ macro programs which call other ImageJ macro programs. All of the parameter flags available in the called program are automatically propagated up to the calling program and made available to be set.
- diffusion2 : Uses VTK to diffuse (and visualize) random points in a volume.
- find_blobs : Uses maxima and zero crossings of the 2nd derivative to find blobs and their spatial extent; works in multiple resolutions.
- thin3Ds: This program is for thinning (i.e., finding the skeleton or medial axis) 3D objects. It is modified from Lee, Kashyap, and Chu (a parallel algorithm). thin3Dp is the same but in parallel.
- channel_current2 : Calculates the current which would be produced when each object in a Ryanodine labeled image releases calcium (e.g, a "spark") which activates nearby BK channels also labeled in the image (in a different wavelength).
- morph3d : Performes morphological processing on grayscale images.
- pwarpL : pwarp warps an image into a new image based upon a set of registered points from one image to the other. A line segment can also be specified as an old "point". The new point will match the best position along this line segment.

## Shrinkwrap - a Deformable Surfaces program

We have developed and implemented a deformable surface model [Lifshitz, L. M., Fogarty, K. E., Gauch, J. M. & Moore, E. D. W. "Computer vision and graphics in fluorescence microscopy". Visualization in Biomedical Computing, SPIE vol. 1808, 521-534 (1992)], called "shrinkwrap" which has been used numerous times over the years to locate the plasma membrane of fluorescently labeled cells. The model is composed of vertices spaced along a model surface, which is initially placed within a 3D image, near its correct location. Each vertex is then allowed move in response to image intensity gradients, attempting to move to locations in the image with a higher intensity while imposing a cost based upon how much curvature its position would impose onto the surface. Our implementation calculates curvature measures (i.e., first and second derivatives of the surface) and intensity measures locally and quickly. This results in a fast algorithm, which is also easy to modify to handle additional local constraints.

Other programs related to finding the surface of a smooth muscle cell and performing related analysis:

- Corecell.pl : Takes a 3D image and an outline and keeps all data points within the specified distance of the outline (removes the "core" of the cell). It runs shrink_wrap and uses its results to do this.
- histogram2: This program creates a number of histograms of signal intensity (e.g., Calcium label) or volume as a function of distance of the object voxels from the nearest background voxel. If the "background" voxel represents the plasma membrane of a cell, the volume (or signal intensity, eg, Calcium) histogram this program produces describes, e.g., what percentage of the cell volume (or, eg, Calcium) is within 1 um of the plasma membrane. The volume histogram can be used to determine if intracellular objects are distributed randomly within the cell. Changes in the signal intensity histogram over time can indicate movement of signal (e.g., Calcium) within the cell.
- analyze_histo: Rather than apply a single threshold value to an entire image, this program produces a different threshold for each distance from the plasma membrane. It analyzes the histograms produced by histogram2 to calculate these thresholds. Different thresholds may be needed since signal intensity levels may vary as a function of distance from the plasma membrane and therefore noise levels (due to non-specific signal and/or Poisson noise) may also depend upon that distance.
- map_to_surface : This program takes a surface (which are those voxels which are nonzero) and one (or two) data images. It then maps (moves) the data voxels (nonzero voxels) to the position of the nearest surface voxel. It can be the first step in deciding (for example) whether cellular structures just below the plasma membrane are distributed randomly "on" the plasma membrane (a not infrequent question posed by physiologists).
- distance2surface_overtime.pl : Finds the distance (in pixels) of a specified point (over time) from the nearest pixel greater than 0 in a 3D image. Alternatively, it finds the distance to the boundary surface specified with the -surface option. Written to take a tracked gene in the nucleus and find its distance to the nuclear envelope over time.

## CellSim - simulating chemical reactions and diffusion in 3D cells.

CellSim [Kargacin, G., and F.S. Fay. 1991. "Ca2+ movement in smooth muscle cells studied with one- and two-dimensional diffusion model" Biophys. J. 60:1088-1100. and Zou, H.L., M. Lifshitz, R.A. Tuft, K.E. Fogarty, and J.J. Singer. 1999. "Imaging Ca2+ entering the cytoplasm through a single opening of a plasma membrane cation channel". J. Gen. Physiol. 114:575-588. ] is a program for 3-D reaction and diffusion simulation. Cellsim uses a Forward Time Centered Space fully explicit differencing scheme to solve the associated PDEs. Cellsim can run in parallel using MPI. Cellsim models multiple species of diffusible molecules and their reactions (e.g. ions and their binding targets). The simulation is saved as a series of images at specified time intervals. Images explicitly describing flow (since there can be large flow yet no net change in concentration) can also be saved.

## DAVE - Interactive quantitative visualization of multidimensional image data

The Data Visualization and Analysis Environment [Lifshitz, L. M., Collins, J. A., Moore, E. D. W. & Gauch, J., "Computer vision and graphics in fluorescence microscopy". IEEE Workshop on Biomedical Image Analysis, Seattle, 1994. pp. 166-175. doi:10.1109/BIA.1994.315854] is a real-time interactive software system for displaying and analyzing multi-wavelength, 3-D, time-series images. DAVE simultaneously displays up to three image stacks, in red, green and blue. The images can be displayed as 2-D slices or full 3-D volumes, using several different methods for volume rendering. Surfaces (polygonal) or outlines (polyline) can be simultaneously displayed with volumes. Time series can be interactively reviewed in real-time, with selected time points held on screen for visual reference. Pixels or voxels of co-localization among the three images sets can be visually highlighted, segmented, and quantified50 . DAVE calculates percentage of overlap of regions of images, and the probability that such overlap could have occurred randomly. DAVE can render voxels as an individual cubes when the location or value of an individual voxel is desired. An "intelligent" 3-D cursor automatically snaps to the brightest nearby voxel when clicked. DAVE can remotely synchronize with another instance of DAVE, thus researchers in different locations see and interact with the same exact views.

## EPR - Image Restoration and super-resolution.

Over 20 years ago the BIG pioneered the development of image deconvolution for light microscopy and the resulting EPR algorithm 25 remains one of the best. EPR uses the light microscope's point spread function (PSF) in an iterative process to move out-of-focus light back to its originating position. Unlike many image restoration algorithms it is guaranteed to converge, and when given a finely-sampled point spread function it has been shown to increase lateral image resolution beyond the diffraction limit to less than 100 nm/pixel25. When combined with structured light microscopy (i.e. Structured Light Epifluorescence EPR or SLEEPR) it can achieve resolutions better than 50 nm laterally and 250 nm axially. EPR is computationally demanding. Recently Dr. Bellvé has produced both cluster-enabled and GPU versions of EPR, the latter of which can now restore images in seconds on a single workstation, make EPR for the first time amenable for High-Througput, high-speed and intelligent acquisition modalitis.

## SignalMass - quantifying intracellular calcium microdomains.

SignalMass [ZhuGe, R., Fogarty, K., Tuft, R., Lifshitz, L., Sayar, K., and Walsh, J. "Dynamics of signaling between Ca(2+) sparks and Ca(2+)- activated K(+) channels studied with a novel image-based method for direct intracellular measurement of ryanodine receptor Ca(2+) current". J Gen Physiol 116, 845-864 (2000).] is program for the analysis of fast, time-lapse images of live cells loaded with a fluorescent calcium indicator for spatially and temporal local (microdomain) events (flashes of light) representing the release of calcium from internal stores through events called calcium "sparks" and "puffs".

During a Ca2+ "spark", free Ca2+ and Ca2+ bound to fluo-3 quickly diffuse away from the spark release site as Ca2+ continues to be discharged. To quantify the total fluorescence (and, ultimately the total Ca2+ discharged) arising from the binding of fluo-3 to the discharged Ca2+(i.e., the Ca2+ signal mass), the increase in fluo3/fluo4 fluorescence ) must be collected from a sufficiently large volume (and in 3D) to provide a measure of the total quantity of Ca2+ released.

## Kinetic Analysis of TIRF Time Series Image of Cargo and Compartment

Among key endosomal components implicated in early endosome fusion are proteins such as early endosome antigen (EEA1), adaptor protein, phosphotyrosine interaction (APPL1) (25), and Rabenosyn-5 (26). To determine the role of endosomes containing these proteins, we imaged each one simultaneously with clathrin and Tf, using multi-color Total Internal Reflection Fluorescence Microscopy (TIRF-M).

### Flux of Cargo from the Plasma Membrane Through the Endosomal Compartment System

To visualize the movement of Tf from clathrin-enriched plasma membrane regions into the endosomal system, Tf-DyLight-649 was added to cells co-transfected with TagRFP-T-clathrin and EGFP-tagged Rabenosyn-5, APPL1, or EEA1. Cells were imaged at 1 Hz (one frame/s) continuously for 15 min: 5 min before the addition, 5 min in the presence, and 5 min after the removal of Tf-DyLight-649. These complex time-lapse image sets were analyzed by systematically identifying all individual endosomal structures within the images and measuring the relative amount of Tf specifically associated with them at each time point.

### The Algorithm

1) Images of TagRFP-T-clathrin are convolved using a difference of Gaussians (DOG) filter to eliminate the diffuse signal originating from out-of-focus or auto fluorescence. 2) The positions (x,y) of local 2-D maxima (pixels with intensity > all 8 neighbors) in the filter images that exceed a threshold are identified and mapped onto the unprocessed clathrin images. 3) The mean intenisty of the 5 × 5 pixel square (25 pixels) surrounding each (x,y) position [i.e. clathrin(in)] is recorded. The mean intensity of the one-pixel-wide frame (24 pixels) surrounding the square [i.e. clathrin(out)] is subtracted thus removing the local background: [clathrin(in] - clathrin(out)]. 4) The (x,y) positions of the maximum-intensity pixels (MIPs) derived from the clathrin image are mapped onto the corresponding unprocessed Tf image and the same calculation is performed. Thus, Tf has a positive value only if the signal detected within the 5 × 5 pixel region is higher than that of the immediate surrounding region. 5) The ratios of the values obtained r(x,y,t)=[Tf(in(x,y)) − Tf(out(x,y))/clathrin(in(x,y)) − clathrin(out(x,y))] for each region then is calculated, and the mean R(t) of all regions in the image is plotted for each time point.

In the same cell, the trafficking of Tf into Rabenosyn-5–containing endosomes is measured similarly. The R(t) data are fit with a kinetic model, described in detail in the SI Appendix, “Curve Fits to Ratio Data” section, that incorporates known constants for Tf binding to the TfR, the concentration of free ligand, and the amount of Tf associated with clathrin. Tf associates rapidly and saturably with clathrin-containing regions, displaying kinetic constants consistent with binding to the TfR (SI Appendix, Table S1) (30–32). Interestingly, Tf is detected almost immediately within Rabenosyn-5–enriched regions, with virtually no time elapsing between the detection of Tf in clathrin coated regions and its detection in Rabenosyn-5 endosomes, and continues to accumulate in these endosomes as long as free Tf is present. When Tf is removed from the medium, it disappears rapidly from clathrin-enriched regions because of transfer into the endosomal pathway (SI Appendix, Table S1). Strikingly, the rate of exit of Tf from clathrin-containing membrane regions (0.0053 ± 0.002/s, mean ± SEM, n = 4), is indistinguishable from its rate of entry into Rabenosyn-5–enriched endosomes (0.0037 ± 0.002/s, mean ± SEM, n = 4) (SI Appendix, Table S1), implying that Tf internalized via clathrin-coated pits could be delivered almost directly into Rabenosyn-5–containing endosomes.

### Segmentation of vesicles from background

To generate images in which vesicles or other fluorescence structures can be analyzed without interference from diffuse fluorescence, images were then convolved with a difference of Gaussians (DOG) filter consisting of 1) a small, two dimensional, Gaussian spot with unit area (sigma = 150 nm) that acted as a vesicle matched detector, i.e. an approximation to a near-diffraction limited spot, and 2) a larger, inverted, two dimensional Gaussian (sigma = 300 nm) with negative unit area that estimated and subtracted the local background. The Gaussian smoothed images were visually thresholded (global threshold) to select for pixels belonging to objects (e.g. vesicles) and eliminate areas devoid of signal (but containing noise).

### Transferrin Trafficking

For analyzing the trafficking of transferrin, the thresholded-DOG time series image of an endocytic marker (clathrin, rabenosyn-5, APPL, EEA1) was analyzed to identify individual objects (spots, vesicles) and their central (x,y,t) positions in every time point, by finding all 2-dimensional intensity maxima (pixels brighter than their 8 neighbors) and thresholding the maxima to eliminate spurious peaks. The average local fluorescence within a 5x5 pixel region centered at each (x,y,t) position was obtained in both the running-averaged, background-subtracted endocytic marker image and the corresponding transferrin image. From this value, the average fluorescence of a 1 pixel wide frame surrounding the 5 x 5 pixel region was subtracted (e.g. Figure 1b). The fluorescence ratio of transferrin to endocytic marker was calculated for each object position, and the mean ratio and standard error for each time point was calculated (there were typically 100s of objects in each time point) and plotted (e.g. Figure 1c and 3b). The fluorescence ratio was used as an estimate of the (relative) transferrin concentration associated with each object – both signals should be roughly in some proportion to the amount of surface membrane – while also correcting for the exponential decrease in fluorescence brightness with increasing depth in TIRF. Additionally, the kinetics of trafficking of transferrin through clathrin, rabenosyn-5 and APPL1 were modeled, and the models fit to the ratio time course data to determine rates of entry/filling and exiting/emptying of transferrin for each endocytic marker.

### Curve Fits to Ratio Data

The Tf/Clathrin ratio data were fit with a simple kinetic model:

koff Tf + clathrin ↔ Tf-clathrin → ? kon kempty

The TfR is modeled as having two binding sites for Tf, one low affinity and one high affinity. TfRs are assumed to be in instantaneous, steady-state colocalization with clathrin, independent of binding Tf, therefore clathrin is a proxy for TfR. Tf/Clathrin is [Tf-Clathrin] (times a scale factor).The extracellular [Tf] is stepped from zero to a constant concentration, and the increase in the total Tf-TfR, and therefore Tf-clathrin, is proportional to 1-e-k·t where k is the sum of the on and off rates for Tf plus a rate kempty, the rate at which Tf apparently disassociates from clathrin objects. This may be because Tf in fact disassociates from a clathrin object, or a Tf-clathrin object disappears from the TIRF imaging zone (and is replaced by a clathrin-alone object, preserving the steady state). Conversely when the [Tf] is stepped back to zero, the Tf-clathrin decreases as e-k·t where kon·[Tf] is now zero since [Tf] is zero. The two Tf binding sites are assumed to be independent, so the equations for the ratio, R(t)clathrin are

R(t≤t_add)clathrin = 0 R(t_add<t≤t_wash)clathrin = A·{(1-exp(-k1·t)) + (1-exp(-k2·t))}/2 R(t>t_wash)clathrin = R(t_wash)clathrin·{exp(-k1·t) + exp(-k2·t)}/2

where A is an arbitrary constant to match the ratio data amplitude, and

k1 = k1on·[Tf] + k1off + kempty k2 = k2on·[Tf] + k2off + kempty

The equations for R(t)clathrin were fit to the ratio data using a Levenberg–Marquardt least-squares fit algorithm, with t_add, t_wash, [Tf], A, and kempty as parameters. [Tf] is either a constant (t_add< t < t_wash) or 0. Tf binding rate constants used are (Giannetti and Bjorkman, 2004)

k1on=1.0·105 M^-1 s^-1 k1off=3.2·10-3 s^-1 k2on=8.3·105 M^-1 s^-1 k2off=1.0·10-3 s^-1

The Tf-Rbsn (or Tf-APPL) ratio data were fit with a similar model, except that for the association of Tf with rabenosyn-5, we assume that the Tf must pass through the clathrin pathway to be available to taken up in rabenosyn-5 vesicles

Tf-clathrin + Rbsn → clathrin + Tf-Rbsn → ? kfill kempty where kfill = kon·[Tf-clathrin]·[Rbsn]

This yields the differential equations for the change in rabenosyn-5

d[Rbsn]/dt = kempty·[Tf-Rbsn] – kon·[Rbsn]·[Tf-clathrin] d[Tf-Rbsn]/dt = kon·[Rbsn]·[Tf-clathrin] - kempty·[Tf-Rbsn]

[Rbsn]+[Tf-Rbsn] is assumed to be constant. We used R(t)clathrin (see above) as the [Tf-clathrin] driving the association with rabenosyn-5 at time t yielding

d[Rbsn]/dt = kempty·[Tf-Rbsn] – kon·[Rbsn]·R(t)clathrin d[Tf-Rbsn]/dt = kon·[Rbsn]·R(t)clathrin – kempty·[Tf-Rbsn]

[Tf-Rbsn](t) is related to the observed ratio by an unknown scale factor, i.e.

R(t)rabenosyn = B·[Tf-Rbsn](t).

The system of ordinary differential equations was numerically integrated (for a given B, kon and kempty) using the ODE function of Scilab (http://www.scilab.org/) producing [Tf-Rbsn], and hence R(t)rabenosyn . Initial conditions were [Rbsn] =1 and [Tf-Rbsn]=0. This function was fit to the Tf-Rbsn ratio data using the L-M method, with B, kon and kempty as parameters (i.e., the ODE solver

### Programs

### Abbreviations

Tf transferrin TfR transferrin receptor Rbsn rabenosyn-5 kon, koff exponential time constants [] concentration

### References

Rabenosyn-5 defines the fate of the transferrin receptor following clathrin-mediated endocytosis.Navaroli DM, Bellvé KD, Standley C, Lifshitz LM, Cardia J, Lambright D, Leonard D, Fogarty KE, Corvera S. Proc Natl Acad Sci U S A. 2012 Feb 21;109(8):E471-80. doi: 10.1073/pnas.1115495109. Epub 2012 Jan 30.