Calculate drop delay for flow cytometry systems and methods

PRIORITY CLAIM

p-0002This patent application claims priority pursuant to 35 U.S.C. §119(e) to U.S. Provisional Patent Application 61/311,956, filed on Mar. 9, 2010, which is hereby incorporated by reference in its entirety.

BACKGROUND

p-00031. Field of the Invention

p-0004The present invention relates to systems and methods for flow cytometry. More specifically, the present invention relates to systems and methods for calculating drop delay in a flow cytometer.

p-00052. Background Art

p-0006In the field of flow cytometry, a hydrodynamic focusing technique is used to separate and align analyte particles within a stream. An oscillator perturbs the stream as it is released from a nozzle of the flow cytometer. The perturbation results in droplets containing individual analyte particles. Just prior to breaking off from the stream, the droplets can be electrically charged and an electrical field downstream from the break-off point can be used to deflect and sort the charged particle. As part of this process, it is necessary to determine the system's “drop delay” (i.e., the time between when a particle of interest is identified, to when a droplet containing the particle of interest reaches the charging location (or “break-off point”)). As such, the droplet of interest can be charged without affecting neighboring droplets.

p-0007Commercial flow cytometers typically use empirical, iterative, trial-and-error protocols to determine drop delay. For example, fluorescent beads can be run through the flow cytometer while a user (or automated program) cycles through a range of drop delay settings. When the user (or automated program) finds the drop delay setting resulting in the optimal sorting of the fluorescent beads, such drop delay is set as the operational drop delay, and is used for the sample runs. Use of fluorescent beads to determine drop delay, however, is undesirable because fluorescent beads can contaminate the flow cytometer. Further, what are needed are systems and methods to avoid the use of iterative, trial-and-error operations, which are time-consuming and can lack the necessary precision of drop delay calculation.

p-0008The systems and methods presented herein act to solve some of the problems seen in the prior art. It is noted that the systems and methods presented herein do not merely automate a previously manual process. The systems and methods presented herein have been developed through painstaking design, testing, and optimization. Challenges in achieving the necessary precision were identified and addressed in the development of these systems and methods. As such, the systems and methods presented herein provide significantly higher accuracy than prior art systems.

BRIEF SUMMARY

p-0009The invention relates to a method of calculating a drop delay in a stream of a flow cytometer. In one embodiment the method includes the steps of determining a widths plot by measuring the width of a stream based on an image of the stream and setting drop delay based on the widths plot.

p-0010In another, embodiment the method of calculating a drop delay in a stream of a flow cytometer includes identifying a laser position relative to a frame of reference; determining break-off point; and calculating drop delay by setting a strobe phase, imaging the stream, determining a widths plot by measuring the widths of the stream based on the image of the stream, determining the peaks in the widths plot, calculating phase of laser intercept, and setting drop delay in response to the peaks in the widths plot, a phase of the laser intercept relative to a droplet generation wave and a charge phase. In one embodiment the method includes smoothing the widths plot. In another embodiment the method includes calculating a derivative of the widths plot, and recording a plurality of zero-crossings of the derivative, wherein each zero-crossing is indicative of a peak in the widths plot,

p-0011In yet another embodiment the break-off point is the shortest break-off point. In still yet another embodiment the step of imaging the stream includes the steps of running a camera along a length of the stream, while taking a plurality of images of the stream; and stitching the plurality of images together to form a single image. In another embodiment the strobe phase is set to the charge phase plus 180°. In still another embodiment the drop delay is also set in response to a bias correction factor. In yet another embodiment the bias correction factor is equal to a linear function aX+b, where X is the time from nozzle to break-off and parameters “a” and “b” are determined from experimental measurement of system response. In another embodiment the smoothing of the widths plot is conducted using a discrete Fourier analysis. In one embodiment a drive voltage is increased from an operational value when the camera is at about the top of the stream, and the drive voltage is decreased to an operational value as the camera is driven down towards the break-off point. In another embodiment the calculation of drop delay is used to determine when an analyte of interest has reached the break-off point of the stream.

p-0012In another aspect the invention relates to a flow cytometer system for automatically calculating drop delay. In one embodiment the flow cytometer system includes means for determining a widths plot based on an image of the stream and means for calculating drop delay based on the widths plot. In another embodiment the flow cytometer system includes means for taking an image of a stream; means for determining a widths plot based on an image of the stream; means for determining peaks in the widths plot and means for calculating drop delay based on the determined peaks. In one embodiment the system includes means for smoothing the widths plot. In another embodiment the system includes means for calculating a derivative on the widths plot; and means for determining zero-crossing on the widths plot. In yet another embodiment the means for taking an image of the stream comprises taking a plurality of images of the stream and the cytometer system further comprises means for stitching the plurality of images of the stream to form one entire image of the stream.

p-0013In another embodiment the system includes a camera adapted to take an image of a stream; and a processor adapted to determine a widths plot based on the image of the stream, and adapted to calculate drop delay based on the widths plot. In still yet another embodiment the system includes a camera mounted on a precision-calibrated, movable camera stage, which is adapted to take a plurality of images of a stream; and a processor adapted to stitch the plurality of images of the stream to form one entire image of the stream, determine a widths plot based of the one entire image of the stream, smooth the widths plot, calculate a derivative on the smoothed widths plot, determine zero-crossing on the smoothed widths plot, and calculate drop delay based on the determined zero-crossings.

BRIEF DESCRIPTION OF THE DRAWINGS

p-0014The accompanying drawings, which are incorporated herein, form part of the specification and illustrate embodiments of flow cytometer systems and methods for calculating drop delay. Together with the description, the drawings further serve to explain the principles of, and to enable a person skilled in the relevant art(s) to make and use the systems and methods described herein. In the drawings, like reference numbers indicate identical or functionally similar elements.

p-0015<figref idrefs="DRAWINGS">FIG. 1</figref> is a schematic of a prior art flow cytometer system.

p-0016<figref idrefs="DRAWINGS">FIG. 2</figref> is a schematic view of a flow cytometer system in accordance with an embodiment of the present invention.

p-0017<figref idrefs="DRAWINGS">FIG. 3</figref> is a stitched image prepared from a plurality of images taken from a camera in a flow cytometer.

p-0018<figref idrefs="DRAWINGS">FIG. 4</figref> is an example widths plot measured from a stitched image, such as the image shown in <figref idrefs="DRAWINGS">FIG. 3</figref>.

p-0019<figref idrefs="DRAWINGS">FIG. 5A</figref> shows the original widths plot of <figref idrefs="DRAWINGS">FIG. 4</figref> and a filtered (smoothed) widths plot. The original widths plot has been offset by 60 to allow both plots to be shown without overlap.

p-0020<figref idrefs="DRAWINGS">FIG. 5B</figref> is a close-up of the plot of <figref idrefs="DRAWINGS">FIG. 5A</figref>, showing data near the nozzle, overlaid to show both raw and filtered (smoothed) data.

p-0021<figref idrefs="DRAWINGS">FIG. 6</figref> shows how a charge window (bottom) is at a specified phase of the oscillator cycle (top).

p-0022<figref idrefs="DRAWINGS">FIG. 7</figref> is a high-level flowchart showing a method presented herein.

p-0023<figref idrefs="DRAWINGS">FIG. 8</figref> is a flowchart showing a method of finding a laser position.

p-0024<figref idrefs="DRAWINGS">FIG. 9A</figref> is a flowchart showing a method of conducting an initialization scan.

p-0025<figref idrefs="DRAWINGS">FIG. 9B</figref> is a flowchart showing a sub-protocol of the method shown in <figref idrefs="DRAWINGS">FIG. 9A</figref>.

p-0026<figref idrefs="DRAWINGS">FIG. 9C</figref> is a flowchart showing a sub-protocol of the method of shown in <figref idrefs="DRAWINGS">FIG. 9A</figref>.

p-0027<figref idrefs="DRAWINGS">FIG. 10A</figref> is a flowchart showing a method of conducting a calibration scan.

p-0028<figref idrefs="DRAWINGS">FIG. 10B</figref> outlines the steps of determining drop delay in accordance with one embodiment of the present invention.

p-0029<figref idrefs="DRAWINGS">FIG. 10C</figref> is a sub-protocol to capture and stitch a full stream image in accordance with one embodiment of the present invention.

p-0030<figref idrefs="DRAWINGS">FIG. 10D</figref> is a sub-protocol to calculate peaks and valleys of a widths plot in accordance with one embodiment of the present invention.

p-0031<figref idrefs="DRAWINGS">FIG. 10E</figref> is a sub-protocol to calculate drop delay in accordance with one embodiment of the present invention.

DETAILED DESCRIPTION

p-0032The following detailed description of systems and methods for calculating drop delay refers to the accompanying drawings that illustrate exemplary embodiments. Unless otherwise noted, all embodiments and examples should be considered prophetic examples. Other embodiments are possible. Modifications can be made to the embodiments described herein without departing from the spirit and scope of the present invention. Therefore, the following detailed description is not meant to be limiting. Further, it would be apparent to one of skill in the art that the systems and methods described below can be implemented in many different embodiments of hardware, software, and/or firmware. Any actual hardware, software, and/or firmware described are not meant to be limiting.

p-0033<figref idrefs="DRAWINGS">FIG. 1</figref> is a schematic of a prior art flow cytometer system, as taught in U.S. Pat. No. 5,643,796, which was incorporated by reference above. <figref idrefs="DRAWINGS">FIG. 1</figref> is presented to show the basic parts of a flow cytometer system, which can be employed in the systems and methods of the present invention. <figref idrefs="DRAWINGS">FIG. 1</figref>, however, is also presented to show a system that employs an iterative, trial-and-error process to calculate drop delay, which is different from the systems and methods of the present invention.

p-0034Specifically, <figref idrefs="DRAWINGS">FIG. 1</figref> shows a droplet flow cytometer having a nozzle system <b>42</b>, which in turn includes a nozzle container <b>2</b>. Introduced within a nozzle volume <b>3</b> is sheath fluid from sheath reservoir <b>5</b> through sheath fluid port <b>4</b>. Analyte particles are introduced within the sheath fluid from substance reservoir <b>8</b> through substance introduction port <b>9</b>. The sheath fluid, together with its entrained analyte particles, are hydrodynamically focused so that the individual analyte particles line up in a single-file orientation and are emitted from nozzle exit <b>6</b> into free fall area <b>7</b>.

p-0035In order to form regular droplets, an oscillator <b>10</b> is adapted to cause pressure waves within, or actual movement of, the sheath fluid. Oscillator <b>10</b> is controlled by an alternating voltage source <b>22</b>. As a result of oscillations within the sheath fluid, stream <b>12</b> is periodically perturbed to form droplets <b>13</b>. Each droplet <b>13</b> is analyzed while in stream <b>12</b> to assess whether the droplet <b>13</b> contains an analyte of interest.

p-0036As shown in <figref idrefs="DRAWINGS">FIG. 1</figref>, substance property sensor <b>62</b> is positioned to sense emissions occurring from the analyte within stream <b>12</b> just prior to free droplet formation (i.e. at the point where the droplet breaks free from the stream). The signals sensed by substance property sensor <b>62</b> are processed by a processor, such as analysis equipment <b>15</b>. Droplets <b>13</b> can then be characterized by components including sorting equipment <b>16</b>. As shown, sorting equipment <b>16</b> creates an electrostatic field so that charged droplets are deflected, whereas uncharged droplets are not. The differential charging of one droplet as opposed to another can be achieved through the interconnection of analysis equipment <b>15</b> with electrode <b>72</b>. Since analysis equipment <b>15</b> can obtain its information from substance property sensor <b>62</b>, electrode <b>72</b> can be responsive to substance property sensor <b>62</b>. Electrode <b>72</b> can be positioned within nozzle volume <b>3</b> so as to cause a charge on the sheath fluid. Thus, the sheath fluid can be responsive to electrode <b>72</b>. The charge can be transmitted to the end of stream <b>12</b> so that at the instant when a droplet is formed, that droplet can be charged. Electrode <b>72</b> can then be quickly neutralized so that the next droplet formed does not get inadvertently charged. In this fashion, when droplets <b>13</b> pass through sorting equipment <b>16</b>, some individual droplets—namely those which do or do not contain the desired substance as determined by analysis equipment <b>15</b>—will be deflected. The deflected droplet can thus be collected in sample collection container <b>82</b>, whereas the undeflected droplets can pass to drain receptacle <b>81</b>. The appropriate timing of the charging of stream <b>12</b> through action of electrode <b>72</b> requires the accurate utilization of a droplet formation time delay (or “drop delay”) estimate; specifically, the amount of time it takes for the droplet to form after substance property sensor <b>62</b> and analysis equipment <b>15</b> have actually determined whether a droplet <b>13</b> in stream <b>12</b> contains a desired analyte.

p-0037To determine optimal drop delay, the system of <figref idrefs="DRAWINGS">FIG. 1</figref> employs an automatic droplet sampler <b>73</b> that samples the droplets sorted by sorting equipment <b>16</b>. Automatic droplet sampler <b>73</b> has a plurality of collection areas <b>75</b> situated on or supported by movable platter <b>74</b>. Through moving movable platter <b>74</b>, different containers <b>75</b> can be sequentially positioned in free fall area <b>7</b> and thus can gather droplets. Importantly, movable platter <b>74</b> can be responsive to stepper <b>76</b> so that the various containers <b>75</b> can be separately and incrementally positioned in free fall area <b>7</b>. Thus, throughout one single sampling event, a number of different samples can be taken.

p-0038In the system of <figref idrefs="DRAWINGS">FIG. 1</figref>, a droplet formation time delay adjustment <b>80</b> coordinates with automatic droplet sampler <b>73</b> to optimize the time delay. Optimization is achieved by creating an automatic substance maximization means—essentially, any component, software, or circuitry which automatically provides information to allow maximization of the substance sorted. The substance maximization means can include increment control <b>77</b>, a device which acts to automatically adjust the various time delays actually applied by analysis equipment <b>15</b> in controlling electrode <b>72</b>. Through increment control <b>77</b>, the time delay applied by analysis equipment <b>15</b> to electrode <b>72</b> can be varied in incremental amount as each different container <b>75</b> is positioned within free fall area <b>7</b>. Thus, each container <b>75</b> can represent a different time delay amount. These results can then be used to achieve an optimum (or more optimum) drop formation time delay estimate.

p-0039In such an automated configuration, a cell counter <b>79</b> yields information that allows time delay adjustment <b>80</b> to alter the actual time delay used during processing. In this fashion, different time delay amounts can be applied and tested automatically to see which yields the maximum count.

p-0040As would be appreciated by one of skill in the art, despite the fact that the system of <figref idrefs="DRAWINGS">FIG. 1</figref> is an automated system, such system proceeds through an iterative, trial-and-error protocol to determine optimal time delay.

p-0041<figref idrefs="DRAWINGS">FIG. 2</figref> is a schematic view of a flow cytometer system <b>200</b> in accordance with an embodiment of the present invention. As with other flow cytometer systems, flow cytometer <b>200</b> includes a nozzle system <b>242</b>, an oscillator <b>210</b>, and an interrogation chamber enclosure <b>209</b>.

p-0042Within interrogation chamber enclosure <b>209</b>, one or more laser sources <b>208</b> emit one or more laser beams that cross a stream <b>207</b> at laser intercept <b>250</b>. Detectors (not shown) are used to detect emissions from analyte particles in stream <b>207</b>. An analysis system (not shown) then determines whether the analyte particle crossing laser intercept <b>250</b> is an analyte of interest. Such analysis is typically based on the fluorescence, side scatter, or forward scatter of the analyte particle.

p-0043As in typical flow cytometers, oscillator <b>210</b> perturbs stream <b>207</b> to form droplets <b>213</b>. In the methods of the present invention, oscillator <b>210</b> serves as the reference frequency and reference phase for the calculations of the drop delay. Although illustrated as being connected to the bottom of nozzle system <b>242</b>, it is understood that the oscillator can be attached to other locations on the nozzle system <b>242</b>. The time between when an analyte of interest crosses laser intercept <b>250</b> and the time the analyte of interest reaches break-off point <b>214</b> is the “drop delay,” which is typically referenced in periods relative to the frequency of the oscillator <b>210</b>. Typically, the system has about 400 microseconds to determine whether the analyte particle crossing laser intercept <b>250</b> is an analyte of interest.

p-0044Unique to the systems of the present invention, flow cytometer <b>200</b> includes a precision-calibrated, movable camera stage. In one embodiment, the movable camera stage is a ball screw linear motion stage <b>201</b>. Mounted on stage <b>201</b> is a camera <b>202</b>. Typically, camera <b>202</b> is a high magnification objective camera. In operation, camera <b>202</b> is moved along the track of stage <b>201</b>, and the camera position is measured with stage position measurement capabilities, such as optical encoder <b>203</b>. Movement of camera <b>202</b> is actuated by a precision motor or drive mechanism, such as a servo motor <b>204</b>.

p-0045As further explained below, camera <b>202</b> is precisely driven down track <b>203</b> in order to obtain a plurality of high-resolution images of stream <b>207</b>. Specifically, camera <b>202</b> is incrementally driven down track <b>203</b>, taking an image at each increment. The images are then stitched together using digital image processing algorithms. In an alternative embodiment, a wide-angle lens camera can be employed to take an image (or multiple images) of the entire stream <b>207</b>. As will be fully understood below, while prior art systems use iterative, trial-and-error protocols to determine drop delay, the systems and methods of the present invention calculate drop delay based on the stitched image of stream <b>207</b>. In essence, the methods presented herein process the stitched image to identify the number of droplets, peaks, and/or periods between laser intercept <b>250</b> and break-off point <b>214</b>. After identifying the number of droplets, peaks, and/or periods between these two points, the system can accurately predict the drop delay by compensating for variables that affect drop delay. As such, the systems and methods of the present invention do not require fluorescent beads, or empirical, iterative protocols.

p-0046<figref idrefs="DRAWINGS">FIG. 7</figref> is a high-level flowchart showing a method presented herein. The method of <figref idrefs="DRAWINGS">FIG. 7</figref> includes four main sub-methods, namely: a method <b>800</b> for precisely finding the laser position of the flow cytometer; a method <b>900</b> of conducting an initialization scan; a method <b>1000</b> of conducting a calibration scan, in which the drop delay is calculated; and a method <b>1100</b> for maintaining the drop delay. Method <b>1100</b> for maintaining drop delay is known to those skilled in the art, is not pertinent to the present invention, and is therefore not discussed in any detail in this specification.

p-0047<figref idrefs="DRAWINGS">FIG. 8</figref> is a flowchart showing method <b>800</b> for precisely finding the laser position of the flow cytometer. Typically, upon putting together a flow cytometer, a user does not care where in physical space the laser is located. The user typically only cares that the laser is properly aligned with the stream. For the methods of calculating drop delay in accordance with this invention, however, the exact location of the laser's intercept with the stream is important because the location is a variable in calculating drop delay. With reference to flow cytometer <b>200</b> of <figref idrefs="DRAWINGS">FIG. 2</figref>, method <b>800</b> begins at step <b>802</b> by first moving camera <b>202</b> to the top of travel (i.e., where stream <b>207</b> exits nozzle <b>242</b>). At the top of travel, the precise location of camera <b>202</b> is known. The location of camera <b>202</b> at the top of travel serves as the reference point for all ensuing calculations. In step <b>804</b>, the state of the laser shutters is saved such that only the laser of interest in turned on. In step <b>806</b>, the camera's <b>202</b> brightness, contrast, and/or exposure are set for laser viewing. In step <b>808</b>, a snapshot of the laser is taken, and the image is processed into a binary image. In step <b>810</b>, the centroid of the laser is measured as a function of image pixels. With the measurement of the laser centroid taken, the exact position of the intercept of the laser with the stream, relative to the camera stage, is known.

p-0048<figref idrefs="DRAWINGS">FIG. 9A</figref> is a flowchart showing method <b>900</b> of conducting an initialization scan. Method <b>900</b> includes a step <b>902</b> for scanning the oscillator <b>10</b> through a range of frequencies to determine the shortest droplet break-off; a step <b>904</b> to confirm that the laser spot location and size are valid (or plausible); and a step <b>906</b> to calibrate the pixel size in microns per pixel.

p-0049<figref idrefs="DRAWINGS">FIG. 9B</figref> is a flowchart showing sub-protocol step <b>902</b> of method <b>900</b> shown in <figref idrefs="DRAWINGS">FIG. 9A</figref>. Sub-protocol step <b>902</b> can be programmed into a processor for determining the shortest droplet break-off of a stream exiting the flow cytometer nozzle. In general what is done in this sub-protocol is the oscillator <b>210</b> (<figref idrefs="DRAWINGS">FIG. 2</figref>) is set at given frequency and the stream <b>207</b> is scanned by the camera <b>202</b> to determine where the drops <b>213</b> break off from the stream <b>214</b>. Once this location <b>214</b> is found, the oscillator frequency is changed by a given increment and the stream is scanned again by the camera. This process is repeated until the break-off point stops moving toward the nozzle and begins to move away from the nozzle. The frequency which causes the point of closest position of break-off to the nozzle is the optimal frequency. The oscillator frequency is then held at this optimal frequency and the stream is scanned again at finer resolution until the point <b>214</b> (<figref idrefs="DRAWINGS">FIG. 3</figref>) at which break-off just begins is determined.

p-0050In more detail in step <b>911</b>, the program sets upper and lower oscillator drive frequency limits for the oscillator frequency scan. In step <b>913</b>, the camera is moved below the nozzle position. In step <b>915</b>, the scan step size is initially set to 1200 Hz. In alternative embodiment, the scan step size can be set to an alternative frequency. In step <b>917</b>, a quantile list of the position of first drop separation is initialized in memory so as to record the subscans of the shortest 70% of the course scan and the shortest 75% of the next finer scan. In step <b>919</b>, the program asks if all scans are done. If they are, the program ends. Otherwise, the program proceeds to step <b>921</b>, where the program scans the stream with the camera over each frequency in the oscillator frequency range and records the break-off position at each frequency increment. This series of steps generates a series of measurements which forms a distribution about the point of closest approach of the break-off point to the nozzle. It is this distribution that will be divided into quantiles and the upper and lower quantiles discarded.

p-0051Step <b>921</b> is performed by the protocol outlined in steps <b>921</b>A-<b>921</b>G. Specifically, while the frequency is less than the upper limit <b>921</b>A, the scans continue. At each set oscillator frequency, the camera is moved until the first detached drop is visible, in step <b>921</b>B. Then, in step <b>921</b>C, the average distance (termed a drop wavelength) between detached drops is calculated by averaging the distance between centers of all visible drops. In step <b>921</b>D, the position of the last attached drop is calculated by subtracting one average drop wavelength from the center of the first detached drop. In step <b>921</b>E, the break-off position is calculated by subtracting the center of mass of the last attached drop from the current camera motor position. In step <b>921</b>F, the oscillator frequency and break-off position are recorded. In step <b>921</b>G, the oscillator frequency is incremented, and the protocol returns to step <b>921</b>A.

p-0052Sub-protocol step <b>921</b> then moves to step <b>923</b>, where the scan's break-off positions are sorted and the ones above the specified quantile in the distribution of break-off points are saved. In step <b>925</b>, the scan step size is incremented by one half. In step <b>927</b>, the quantile index into memory is incremented. In step <b>929</b>, the upper and lower frequency limits are set to re-scan the saved list using a finer step size.

p-0053<figref idrefs="DRAWINGS">FIG. 9C</figref> is a flowchart showing sub-protocol step <b>906</b> of the method of shown in <figref idrefs="DRAWINGS">FIG. 9A</figref>. In sub-protocol <b>906</b>, the microns per pixel of image is measured and calibrated. In step <b>930</b>, the microns per pixel are set to a default value (based on optical design nominal values). Steps <b>932</b>, <b>934</b>, and <b>936</b> are then repeated three times. In step <b>932</b>, the position of the camera is set so the top drop is in the lower 25% of the screen. In step <b>934</b>, the camera is moved down a known distance (in motor position coordinates). In one embodiment, the distance moved corresponds to approximately 50% of the display screen (in image rows). In step <b>936</b>, the new position on screen of top drop is measured and the original position is subtracted to determine the number of image rows moved in response to commanded motor move (in microns). After repeating steps <b>932</b>, <b>934</b>, and <b>936</b>, microns per pixel is set to equal the median of the three measured values.

p-0054<figref idrefs="DRAWINGS">FIG. 10A</figref> is a flowchart showing a method of conducting a calibration scan. Next, in general, the position of the break-off point is affected by not only the frequency of the oscillator but also the amplitude of the driving frequency. As a result a series of changes in oscillator amplitude are made and the effect of the changes on the repeatability of the image determined. Similarly the strobing of the image is also varied and the effect of the strobe changes on the repeatability of the image determined. These steps together are termed drop delay maintenance. Step <b>1010</b> (strobe scan to optimize picture for drop delay maintenance) and step <b>1020</b> (amplitude scan to measure stream response for drop delay maintenance) are steps related to the maintenance for the set drop delay (step <b>1100</b>), are known in the art, and are therefore not discussed in any detail herein.

p-0055Step <b>1030</b> provides an exemplary algorithm for determining drop delay. <figref idrefs="DRAWINGS">FIG. 10B</figref> outlines the steps of determining drop delay in accordance with one embodiment. In step <b>1040</b>, the strobe phase (the position of the application of the camera strobe relative to the peaks in the drop wavelength) is set to be the charge phase (the position at which occurs the application of charge to the stream relative to the peaks in the drop wavelength) plus 180°. In practice, the charge phase is typically set to be slightly offset from the reference phase (i.e., the phase of the oscillator drive frequency). In step <b>1050</b>, a movable camera is used to capture a plurality of images of the stream, and a processing unit is used to stitch the plurality of images to form a full stream image.

p-0056<figref idrefs="DRAWINGS">FIG. 10C</figref> is a sub-protocol of step <b>1050</b>, outlining the steps to capture and stitch a full stream image in accordance with one embodiment of the present invention. <figref idrefs="DRAWINGS">FIG. 3</figref> shows an exemplary stitched image, which was prepared from a plurality of images taken from a camera in a flow cytometer in accordance with the present invention. With reference to flow cytometer <b>200</b> of <figref idrefs="DRAWINGS">FIG. 2</figref>, camera <b>202</b> is moved to a home position (i.e., the top of stream <b>207</b>), step <b>1051</b>, where the first image is taken, step <b>1052</b>. In step <b>1053</b>, the stored background image is subtracted for the position. In step <b>1054</b>, camera <b>202</b> is moved down stage <b>201</b> approximately one screen-full, but leaving a specified amount of overlap (e.g., 10% of an image). A subsequent image is then taken. This process is repeated (step <b>1056</b>) until camera <b>202</b> images only droplets <b>213</b> (i.e., camera <b>202</b> has moved past break-off point <b>214</b> and is at the bottom of motor travel).

p-0057Images near the top of stream <b>207</b> can have poorly defined width signatures, which can complicate the extraction of peak information. For this reason, when operating with a low drive amplitude setting of the oscillator, the drop drive voltage setting in one embodiment is increased while taking the first several images. Increasing the drop drive voltage has the effect of increasing the amplitude of the stream perturbations, but does not shift the droplet position. As such, the drop count remains consistent. However, the drive voltage should be decreased to its operational setting before moving the camera too far down the stream because the drive voltage setting will change the absolute break-off of the droplets.

p-0058In one embodiment, there is provided a means for taking a plurality of images of a flow cytometer stream by performing the following program, written in C++ programming language:

p-0059<tables id="TABLE-US-00001" num="00001"><table frame="none" colsep="0" rowsep="0"><tgroup align="left" colsep="0" rowsep="0" cols="2"><colspec colname="offset" colwidth="14pt" align="left" /><colspec colname="1" colwidth="203pt" align="left" /><thead><row><entry /><entry namest="offset" nameend="1" align="center" rowsep="1" /></row></thead><tbody valign="top"><row><entry /><entry>Calibrate micronsperpixel</entry></row><row><entry /><entry>//micronsperpixel is a measurement of the number of microns of</entry></row><row><entry /><entry>stream //length represented by each image pixel</entry></row><row><entry /><entry>Move camera to home</entry></row><row><entry /><entry>do</entry></row><row><entry /><entry> if near nozzle</entry></row><row><entry /><entry>bump up drop drive amplitude</entry></row><row><entry /><entry> else</entry></row><row><entry /><entry>  set drop drive amplitude to operational</entry></row><row><entry /><entry> end if</entry></row><row><entry /><entry> take (horizontally) centered picture</entry></row><row><entry /><entry> if only drops in picture break</entry></row><row><entry /><entry> command motor to move to curpos + (1−overlap)*1280*</entry></row><row><entry /><entry> micronsperpixel</entry></row><row><entry /><entry>while</entry></row><row><entry /><entry namest="offset" nameend="1" align="center" rowsep="1" /></row></tbody></tgroup></table></tables>

p-0060It is understood that the above computer program, and any other computer program presented herein, is merely exemplary. Other logic, in the same or other programming language, can be used to achieve the same purpose.

p-0061The method then moves on to step <b>1057</b>, where the images are stitched by overlapping portions to the position with the highest correlated pixel value, as further discussed below. The plurality of images is stitched together using digital image processing algorithms. Typically, the images are poorly aligned. Vertical alignment is a function of the precision of the motor movement, the consistency of the camera stage (camera pointing) and the estimated value for microns per pixel. The horizontal alignment is a function of motor travel and stream vertical alignment, and tends to be off by several tens of pixels.

p-0062To find the best image match, each of the images gathered in the first phase is fitted to the previous image immediately above to find the best alignment in pixel space. To accomplish the fit, the bottom image is moved with respect to the upper image to find the best location (in top image pixel space) in terms of fit. The images are sampled in a correlation window that is common to both images. In one embodiment, there is provided a means for stitching a plurality of images of a flow cytometer stream by performing the following program, written in C++ programming language, wherein in the algorithm below: (dr) is ½ the width of the row search space, and (dc) is ½ the column search space for the image correlation; the values of rowoffset and rowdepth are chosen such that the image can be move downward and leftward without the correlation window of the second image leaving the first; the function ordered-pixels transforms a rectangular domain of an image (defined by [rowbegin, rowend, colbegin, colend]) into a vector by unwrapping each row scanline consecutively.

p-0063<tables id="TABLE-US-00002" num="00002"><table frame="none" colsep="0" rowsep="0"><tgroup align="left" colsep="0" rowsep="0" cols="1"><colspec colname="1" colwidth="217pt" align="left" /><thead><row><entry namest="1" nameend="1" align="center" rowsep="1" /></row></thead><tbody valign="top"><row><entry>startrow = 1280*(1−overlap)</entry></row><row><entry>// Bottom image pixels</entry></row><row><entry>vb = orderedpixels([rowoffset, rowoffset+rowdepth, width/3, 2*width/3])</entry></row><row><entry>for row = startrow − dr to startrow + dr</entry></row><row><entry>for col = −dc to dc</entry></row><row><entry> // Top image pixels</entry></row><row><entry> vt = orderedpixels([startrow + rowoffset + row,</entry></row><row><entry>           startrow + rowoffset + row + rowdepth,</entry></row><row><entry> width/3+col, 2*width.3+col)</entry></row><row><entry> c=( vt, vb)\(| vt || vb |)</entry></row><row><entry>end for</entry></row><row><entry>end for</entry></row><row><entry>return row and col with max c.</entry></row><row><entry namest="1" nameend="1" align="center" rowsep="1" /></row></tbody></tgroup></table></tables>

p-0064It is understood that the above computer program, and any other computer program presented herein, is merely exemplary. Other logic, in the same or other programming language, can be used to achieve the same purpose.

p-0065The images are often mismatched upon first processing of the images. For example, with a 400×1280 pixel image, and a 10% overlap, it is predicted that the bottom image will fit at pixel offset g=(0, 1152). For all pixels within the box that stays fixed to the second image space, the pixels of the two images are treated as two large vectors (vt and vb), with entries defined by the pixel colors (unrolled by row and column). The correlation factor is determined by computing (and normalizing) the dot product of the two vectors, c=(vt, vb)\(|vt∥vb|). This factor (0≦c≦1) determines how well the two images align for a given offset g. The fit is largest when the two images are exactly the same.

p-0066The factor g is moved within the given correlation window and the best offset is used as a stitched offset for these two images. This process continues down the stream until all images are correlated to the one above. As the process continues, the stitching offsets are kept relative to the first image so as to place the entire sequence into a single image (by simply copying pixels from relative to top image space). Where images overlap, the two sets of pixels are averaged, which allows visual inspection of the image to determine the quality of the stitching.

p-0067Due to the periodic nature of the stream image, an initial guess must be good enough so that it is within a half of a drop wavelength of the correct pixel offset. Otherwise a given waveform might potentially be matched to the wave before or after the correct one, thus either inserting or deleting a period from the system.

p-0068After the image is stitched together, as shown in <figref idrefs="DRAWINGS">FIG. 3</figref>, the width of the stream at each line in the image is extracted into a widths plot (<figref idrefs="DRAWINGS">FIG. 4</figref>) that is used to detect the peaks and valleys of the stream width image. <figref idrefs="DRAWINGS">FIG. 10D</figref> is a sub-protocol to calculate peaks and valleys of a widths plot in accordance with one embodiment of the present invention.

p-0069<figref idrefs="DRAWINGS">FIG. 4</figref> is an example widths plot measured from a stitched image, such as the image shown in <figref idrefs="DRAWINGS">FIG. 3</figref>. Assuming that the speed of the particles is closely tied to the wavespeed, the count of these peaks allows for an estimate of the amount of time a particle takes to traverse from interrogation (or laser intercept <b>250</b>) to break-off point <b>214</b>, and thus allows for the calculation of drop delay.

p-0070In one embodiment, the stitched image is converted to black and white via standard processing (histogram equalization, Gaussian smoothing, and row-by-row mean value thresholding). The resulting black and white image, where white represents the background population and black the stream population, is then converted into a stream widths plot. In one embodiment, there is provided a means for converting a stitched image into a black and white image by performing the following program, which is written in C++ programming language:

p-0071<tables id="TABLE-US-00003" num="00003"><table frame="none" colsep="0" rowsep="0"><tgroup align="left" colsep="0" rowsep="0" cols="2"><colspec colname="offset" colwidth="14pt" align="left" /><colspec colname="1" colwidth="203pt" align="left" /><thead><row><entry /><entry namest="offset" nameend="1" align="center" rowsep="1" /></row></thead><tbody valign="top"><row><entry /><entry>For each row</entry></row><row><entry /><entry> Width(row) = 0</entry></row><row><entry /><entry> // scan from left edge until we hit stream pixels</entry></row><row><entry /><entry>While pixel(row,col) is WHITE and column &lt; width</entry></row><row><entry /><entry>  Next column</entry></row><row><entry /><entry> End while white pixels</entry></row><row><entry /><entry> // now count all contiguous black pixels as part of stream width</entry></row><row><entry /><entry> While pixel(row,col) is BLACK and column &lt; width</entry></row><row><entry /><entry>  Increment width(row)</entry></row><row><entry /><entry> End while black pixels</entry></row><row><entry /><entry>End for</entry></row><row><entry /><entry namest="offset" nameend="1" align="center" rowsep="1" /></row></tbody></tgroup></table></tables>

p-0072It is understood that the above computer program, and any other computer program presented herein, is merely exemplary. Other logic, in the same or other programming language, can be used to achieve the same purpose.

p-0073Since the widths plot is taken from the monochrome image, it is fairly noisy. In particular, near the nozzle where the amplitude of the oscillator on the stream is low, the signal-to-noise ratio of the widths plot is not very good. From row to row, the threshold algorithm results in a few pixels being assigned to the wrong category (either false positive or false negative) resulting in a very noisy signal. Trying to find peaks and valleys on this signal would be very difficult. As such, in one embodiment, a Fourier analysis is conducted to filter the widths plot in order to extract the relevant signal characteristics (peaks and valleys of the main signal). Specifically, a Fourier analysis for frequency domain filtering is conducted to filter (or smooth) the widths plot. Filtering in the Fourier domain allows for the retention of selected frequencies and rejection of others via Fourier transform of the spatial domain signal, application of a simple windowing function to the frequency domain signal, followed by an inverse Fourier transform back to the spatial domain. Such Fourier analysis is outlined in <figref idrefs="DRAWINGS">FIG. 10D</figref>.

p-0074Preprocessing can be useful to effectively filter a discrete signal using a windowing function. In step <b>1081</b>, the widths plot is trimmed to remove the nozzle area of the widths plot (i.e., the left edge of the widths plot). Efficient implementations of the Fast Discrete Fourier Transform operate on data whose length is an even power of two. In step <b>1082</b>, the original data is padded from its original length with zeros to the next higher power of two length. Sudden transition from the signal values to a padding value would introduce step function to the data if care is not taken, and would result in Gibbs ringing of the transformed data. In order to pad the data with as smooth of a transition as possible, the data is extended from the left edge since the measured data at row zero converges (to the left) to a DC offset equal to the mean width of the stream. In one embodiment, there is provided a means for padding the widths plot by performing the following program, which is written in C++ programming language:

p-0075<tables id="TABLE-US-00004" num="00004"><table frame="none" colsep="0" rowsep="0"><tgroup align="left" colsep="0" rowsep="0" cols="2"><colspec colname="offset" colwidth="14pt" align="left" /><colspec colname="1" colwidth="203pt" align="left" /><thead><row><entry /><entry namest="offset" nameend="1" align="center" rowsep="1" /></row></thead><tbody valign="top"><row><entry /><entry>Calculate stream mean width</entry></row><row><entry /><entry>Calculate NN = next power of two length, from widths.size( )</entry></row><row><entry /><entry>Pad = NN − widths.size( )</entry></row><row><entry /><entry>For i=0 to pad</entry></row><row><entry /><entry> paddedWidths[i] = mean width</entry></row><row><entry /><entry>end for</entry></row><row><entry /><entry>for i=pad to NN</entry></row><row><entry /><entry> paddedWidths[i] = widths[i−pad]</entry></row><row><entry /><entry>end for</entry></row><row><entry /><entry namest="offset" nameend="1" align="center" rowsep="1" /></row></tbody></tgroup></table></tables>

p-0076It is understood that the above computer program, and any other computer program presented herein, is merely exemplary. Other logic, in the same or other programming language, can be used to achieve the same purpose.

p-0077In step <b>1083</b>, a discrete Fourier transform is performed on the widths plot. While the flowchart of <figref idrefs="DRAWINGS">FIG. 10D</figref> illustrates the use of a Fourier transform and Hann filter analysis, alternative filtering techniques can be designed to achieve the same objectives. Windowing is then performed to filter the data. The simplest filter window is the rectangular window—zero for stopped frequencies, one for passed frequencies. This window produces ringing in the inverse transformed signal, though, and can make it difficult to resolve either a) signals whose relative strength is very different, or b) signals whose frequencies are too close together. Instead of a rectangular window, an alternative embodiment uses a commonly used moderate windowing function called a Hann window. In step <b>1084</b> a Hann filter window is created and centered at the droplet wavelength. The Hann window (function provided below) reduces spectral leakage compared to the rectangular window, but still has good sensitivity for low-bandwidth signals like the widths plot:

p-0078<maths id="MATH-US-00001" num="00001"><math overflow="scroll"><mrow><mrow><mi>w</mi><mo>⁡</mo><mrow><mo>(</mo><mi>n</mi><mo>)</mo></mrow></mrow><mo>=</mo><mrow><mn>0.5</mn><mo>⁢</mo><mrow><mo>(</mo><mrow><mn>1</mn><mo>-</mo><mrow><mi>cos</mi><mo>⁡</mo><mrow><mo>(</mo><mfrac><mrow><mn>2</mn><mo>⁢</mo><mi>π</mi><mo>⁢</mo><mstyle><mspace width="0.3em" height="0.3ex" /></mstyle><mo>⁢</mo><mi>n</mi></mrow><mrow><mi>N</mi><mo>-</mo><mn>1</mn></mrow></mfrac><mo>)</mo></mrow></mrow></mrow><mo>)</mo></mrow></mrow></mrow></math></maths>

p-0079The center frequency of the Hann window is chosen to match the frequency of the first peak on the Fourier transformed data. The width of the Hann window can be chosen fairly arbitrarily to encompass most of the first peak information for a given nozzle size (e.g., 70 micron nozzles), while excluding the secondary and higher order peaks in the frequency data. The left edge and right edge cutoff values of the filter are symmetric around the peak and cutoff frequencies whose magnitudes are less than about 10% of the magnitude of the peak. The width of the filter window can be parameterized to deal with other nozzle and stream information.

p-0080In step <b>1085</b>, the data is filtered by multiplying the filter window by the frequency domain padded data. Then, in step <b>1086</b>, the frequency domain data is transformed back to the spatial domain with an inverse-transform function. <figref idrefs="DRAWINGS">FIG. 5A</figref> shows the original widths plot of <figref idrefs="DRAWINGS">FIG. 4</figref> and a filtered (smoothed) widths plot. The original widths plot has been offset vertically to allow both plots to be shown without overlap. <figref idrefs="DRAWINGS">FIG. 5B</figref> is a close-up of the plot of <figref idrefs="DRAWINGS">FIG. 5A</figref> showing data near the nozzle, overlaid to show both raw and filtered (smoothed) data. Overall, the filtered data has peaks and valleys that track very closely to the original data, but which have smooth inflections allowing the location of the peak values to be calculated.

p-0081After filtering the peak data, the derivative of the filtered data (by simple difference) is computed, in step <b>1087</b>, and the zero crossings are identified, in step <b>1088</b>. The peaks are identified by selecting the zero-crossings that are decreasing left to right. At the nozzle, an extrapolation is performed because the last waveform or two can be ill formed due to the presence of the nozzle. In one embodiment, there is provided a means for finding peaks by performing the following program, which is written in C++ programming language:

p-0082<tables id="TABLE-US-00005" num="00005"><table frame="none" colsep="0" rowsep="0"><tgroup align="left" colsep="0" rowsep="0" cols="2"><colspec colname="offset" colwidth="21pt" align="left" /><colspec colname="1" colwidth="196pt" align="left" /><thead><row><entry /><entry namest="offset" nameend="1" align="center" rowsep="1" /></row></thead><tbody valign="top"><row><entry /><entry>//let sw be the smoothed width vector.</entry></row><row><entry /><entry>resize(ds, length(sw)−1)</entry></row><row><entry /><entry>for i = 1 to length(ds)</entry></row><row><entry /><entry>ds = sw(i) − sw(i−1)</entry></row><row><entry /><entry>end for</entry></row><row><entry /><entry>for i = 2 to length(ds)</entry></row><row><entry /><entry> if (ds(i−1) >= 0 and ds(i) < 0) push_back(peaks, ds(i))</entry></row><row><entry /><entry>end for</entry></row><row><entry /><entry namest="offset" nameend="1" align="center" rowsep="1" /></row></tbody></tgroup></table></tables>

p-0083It is understood that the above computer program, and any other computer program presented herein, is merely exemplary. Other logic, in the same or other programming language, can be used to achieve the same purpose.

p-0084<figref idrefs="DRAWINGS">FIG. 10E</figref> is a sub-protocol to calculate drop delay in accordance with one embodiment of the present invention. The method of calculating drop delay is outlined as follows. In step <b>1092</b>, the phase of the laser intercept relative to the oscillator frequency is calculated by measuring the position of the widths plot peak before the intercept and the widths plot peak after the intercept. The phase is then set to equal the portion of the wavelength between the peaks that the laser intercept as measured from the previous peak. In step <b>1094</b>, the system is interrogated to identify the phase offset of the droplet stream charge pulse. The charge phase is typically offset from the reference drive frequency to account for the ramp-up time of the charge pulse. In step <b>1096</b>, the drop delay is calculated by setting the drop delay to equal the number of peaks in the filtered widths data, plus the phase of the laser intercept relative to the droplet generation wave, minus the charge phase, plus a bias component. The protocols of steps <b>1092</b>, <b>1094</b>, and <b>1096</b> will be further understood based on the description provided below.

p-0085In one embodiment, a method for calculating drop delay accounts for phase variables such as the phase of the oscillator, charge electrode, strobe system, camera images, and/or laser. Specifically, to obtain the fractional part of drop delay, various phase variables within the system are reconciled. These variables can be ignored when using beads to experimentally determine drop delay, as in prior art systems, but should be accounted for in the predictive methods of the present invention.

p-0086The appropriate charge phase is selected by moving the charge phase window until clean sort streams are obtained. This process effectively aligns the offset between the drop drive clock and charge to optimize droplet charge. A key realization is that the droplet charge pulse can occur only at a discreet set of times, namely once per drop drive period. <figref idrefs="DRAWINGS">FIG. 6</figref> shows how a charge window (bottom) is at a specified phase of the oscillator cycle (top).

p-0087In practice, the basic life cycle for an event is: 1) a timestamp is assigned at the laser intercept; 2) drop delay is added; and 3) the selected droplet is charged. Specifically, at the moment of laser interrogation, a high precision time-stamp is added to the event. The timestamp is passed into software, which estimates the eventual drop that the particle will end up in. The formula for this is: drop=(intercept time+drop delay) div (drop drive period), where div is integer division (throw away remainder). By choosing drop delay, a mapping is created between a moving window in time (exactly one period long) and a drop number. Each subsequent window in time will map to the subsequent drops.

p-0088The effect of changing the fractional part of drop delay is to slide the window in time of events that are predicted to land in a drop. When the calculation of the fractional part is wrong, the system effectively predicts the “pre drop” incorrectly, and events that fall within these times map to more than one droplet, or to a droplet other than the desired.

p-0089With respect to the laser phase, a very simple starting point for illustrating the effect of laser phase is to consider the case when the laser intercepts at the valley between two peaks, for a picture taken at zero strobe phase relative to drop drive period. A basic assumption (not always strictly true) is that such a region centered around a droplet peak eventually contracts (as it moves downstream) to form a droplet (i.e., the “pre-drop”).

p-0090Consider a hypothetical particle passing the laser at an instant. Its timestamp would mark the beginning of a time interval that is desired to map into a particular drop. Particles at the top of the “pre-drop” would pass the laser at one period later, so it is desired to arrange the drop delay such that the entire interval in time maps (under the div operation) to the same drop. This arrangement is the essence of the fractional part of drop delay. The timestamp on such a particle would have a timestamp that is a perfect integer multiple of the drop drive period. Under the div operation, all times that follow within the period would lose their remainder and map to the same drop as this first particle.

p-0091The primary introduction of fractional drop delay is to account for the case when the laser is not exactly at the beginning of a “pre drop” at phase zero. For example, if the laser were instead, at the center of the bulge in the stream. All particles would hit this laser at ½ of a period earlier than they would if the laser were centered on a trough between bulges of the stream. To account for this translation in laser position, in one embodiment, the shift is accounted for by artificially adding ½ A of a period to the ADC time: <br />drop=({intercept time+½ period}+drop delay)div(drop drive period)

p-0092Practically, however, it is not possible to change the ADC time stamp given to an event. Instead, the associativity of addition is used to absorb this fractional offset into the drop delay: <br />drop=(intercept time+{½ period+drop delay})div(drop drive period)

p-0093As such, the delay has a fractional component. If the laser is at some other fractional offset (besides ½), the same procedure is followed. Of course, within a given period it is very unlikely that more than one particle will pass the laser. However, the procedure above still applies since the count of the last attached droplet is incremented by one each period.

p-0094With respect to the charge phase, since a droplet will break off cleanly at some (unknown) phase during the period, it is possible (depending on charge phase) that the actual drop that is charged for sorting can be the one that is still connected in the image at phase zero, or the next higher (for later charge phases). To enable the image processing to clearly detect which is the case, the image is taken at exactly the charge phase. In this way it can be detected which drop is to be selected for charging. Taking the image at a later phase effectively moves the laser up the “pre drop” by the corresponding fractional phase. To make up for this effective move, the algorithm effectively moves the laser back down by subtracting the fractional movement, and the final drop delay estimate is: <br />drop delay={fractional laser offset+integral drop delay−fractional charge phase},<br /> where the image is strobed at exactly the charge phase

p-0095A final accommodation can be made for an experimentally observed bias. For alternative nozzle sizes, a bias can exist. Such bias can be experimentally determined and can be deduced to an algebraic function in order to account for the bias (bias=aX+b). For example, the bias for a 70 micron nozzle, for example, is set as a function of: <br />−0.015465<i>X+</i>0.891189, where <i>X=</i>0.1(break-offMicrons−charge/360*lastWaveLenMicrons)/avg_velocity.

p-0096This variable is basically the time from nozzle to break-off in 10 s of microns (roughly a period at 100 k Hz).

p-0097The breadth and scope of the present invention should not be limited by any of the above-described exemplary embodiments, but should be defined only in accordance with the following claims and their equivalents.