extfieldsofstudy:
Medicine
Computer Science

text:
Large-scale analysis of high-speed atomic force microscopy data sets using adaptive image processing

Modern high-speed atomic force microscopes generate significant quantities of data in a short amount of time. Each image in the sequence has to be processed quickly and accurately in order to obtain a true representation of the sample and its changes over time. This paper presents an automated, adaptive algorithm for the required processing of AFM images. The algorithm adaptively corrects for both common one-dimensional distortions as well as the most common two-dimensional distortions. This method uses an iterative thresholded processing algorithm for rapid and accurate separation of background and surface topography. This separation prevents artificial bias from topographic features and ensures the best possible coherence between the different images in a sequence. This method is equally applicable to all channels of AFM data, and can process images in seconds.

Introduction
Atomic force microscopes (AFMs) are a useful tool for investigating nanoscale surfaces. They have applications in physics, materials science, chemistry, biology and nanotechnology. AFMs generate detailed three-dimensional images of surfaces with nanometer and subnanometer resolution [1][2][3][4][5][6][7][8][9][10][11][12]. The raw imaging data has to be post-processed to eliminate artifacts arising from distortions inherent in the technique or a specific instrument. Sample tilt, scanner bow and other artifacts corrupt the true topography. Removal of these artifacts is requisite to obtain the true topography of a sample [13][14][15][16][17]. In most AFM studies, the goal is to generate a single image of a sample surface or surfaces of multiple samples. The resulting data sets typical distortions on model AFM data. Panels A though E show a sample with arbitrary topographic features (randomly placed and sized half spheres) and the corresponding histograms below. The first column shows the pure topography. The second column shows the addition of random noise. The third column shows the effects of an arbitrary sample tilt. The fourth column shows the effects of an additional second-order polynomial representing scanner bow. The fifth column adds arbitrary 1-D offsets in the fast-scan direction. Each artifact successively broadens the peaks in the histogram. The histogram in column A uses the vertical scale on the left. The histograms columns B through E share a common vertical scale shown on the right. All images share a common color scale shown on the right. All the histograms are plotted over this same range.
are relatively small and easy to correct by hand. In an emerging part of the field, the goal is no longer to see just the detailed surface topography; but rather, to also see how this topography changes as a function of time or treatment [18][19][20][21][22][23][24]. High-speed AFM (HS-AFM) focuses on reducing imaging time to track faster dynamics than can be observed with existing instruments [25][26][27][28][29]. Where conventional AFM generates only a few images, HS-AFM generates hundreds of images. This amount of data creates a significant increase in the amount of image processing needed to extract the true sample topography, which can no longer be performed by the user on an image-by-image basis [15,17]. In HS-AFM, there is a real need for an automated processing routine that optimizes the image processing for each image without losing essential topographic information. One of the most important considerations in automating processing is that it must not introduce processing artifacts, and it must maintain coherence between images in a dataset to allow for rapid comparisons between frames [30,31]. Automated, AFM quality-control check in semiconductors and other nanotechnologies [32][33][34] also creates large data sets, which would also benefit from reliable automated image processing. In order to automate the image processing, it is important to define a metric by which the success of the image processing can be judged, and for which the processing is optimized. In the next sections we will evaluate the source and effect of typical AFM image distortions and derive a suitable optimization metric.

Sources of inherent distortions
All AFM images will contain artifacts in the x-y-direction as well as the z-direction. This work will address only the distortions found in the z-direction. There are two fundamental classes of image distortions in the z-direction, i.e., 2-D distortions and 1-D distortions. Both of these distortions are wellknown and have been commented on since the earliest prototypes. 2-D distortions manifest over the whole image and deform 2-D surfaces into 3-D. The most common example of this is sample tilt. Tilting the sample normal relative to the scanner normal gives the entire image an apparent tilt. Piezo tube scanners (which are often used in AFM instruments) also generate inherent distortions in the image, creating an additional bow in the apparent topography [13]. Some of the early work in the field addresses the need to correct these background distortions appropriately by fitting only the background regions to a polynomial of the correct order [35]. 1-D distortions cause relative offsets from scan line to scan line. The sources of these distortions include laser-mode hopping and changes in the tip-sample interaction. These 1-D distortions cause apparent discontinuities in the topography of the sample, but do not represent an actual topographic feature. Most importantly, 1-D distortions generally appear as offsets in the z-data. Potential corrections for these types of artifacts have been addressed well in a direct fashion from starting principals by Starink and Jovin [35], and again in a less direct fashion that has some potential pitfalls with an emphasis on biological samples [36]. Figure 1 shows how these types of artifacts distort an artificial image along with corresponding histograms. In Figure 1A, half-spheres of arbitrary size were placed randomly within the image. A histogram of the z-heights in the image shows a spike at the background level, and a distribution of heights from the sphere (along the axis x-axis). Figure 1B shows the addition of a small amount of random noise added in the z-direction. This noise leads to broadening of the histogram peak. The peak now appears Gaussian with a small distribution. Figure 1C imparts a small tilt to the sample (less than 10 nm over 3 µm in either direction, which is very good for most samples). The background peak has been broadened so far as to be nearly indistinguishable. Figure 1D adds a second-order polynomial distortion to simulate scanner bow (less than 10 nm over the scan range, which is well within reason for most AFMs). The peak is no longer visible in the histogram. Finally, Figure 1E adds random vertical offsets to some of the scan lines. These offsets simulate the sorts of line skips commonly observed while imaging. As can be seen from the histograms, each successive distortion broadens the peak in the histogram. We will, therefore, test the suitability of calculating the distribution of the histogram as a metric for judging the success of automated-image-processing steps.

Basic requirements for successful image correction
The goal of the AFM image-processing algorithm described in this paper is to remove imaging artifacts and generate the best possible representation of the true sample topography. The achievement of this goal has one basic prerequisite, an area that can be used as a reference geometry. In most cases, this will be an area that can be considered as flat in the true topography of the sample. The image-processing algorithm then attempts to extract the global distortions in the image from observations of the flat reference region. These global distortions are subtracted from the complete image. In practice, this criterion is often met in AFM images, since the sample preparation generally utilizes flat substrates, such as mica or silicon wafers. In the case of semiconductor processing, the reference geometry could also be a specific area of the processed device that is flat or has a known shape. For the purpose of this work, we use the assumption that there is some part in the image that corresponds to a flat background, since this is the most common case. The algorithm automatically determines this background region and performs the necessary background correction. The algorithm uses minimization of the standard deviation of the image histogram to judge the success of the background detection and the image processing.

Image term description
Throughout the rest of this paper we will use the following conventions when referring to images for the sake of clarity. First, an image is defined as an array of scan lines along the slow scan axis, where each scan line is an array of pixels in the image along the fast scan axis. (2)

Standard deviation as a metric for image flatness
In order to show that the standard deviation of an image is a suitable metric for monitoring the progress of image flattening we can use uncertainty propagation. We can start by describing the AFM image, , as the sum of two height fields, the sample topography, , and the inherent geometric distortions, : The square variance of the image can be written as: The equation simplifies to the final form because the covariance of the two fields is zero (they are uncorrelated variables for the types of distortions described in this paper). This means that ; thus, the geometric distortions will always broaden the distribution of observed heights relative to the true sample topography. This also means that as long as the processing steps only affect (flattening is only performed on truly flat regions, and not real topography), the minimization of will improve the image.

Commonly used methods for correction of inherent distortions
The goal of AFM image processing is to correct the inherent distortions mentioned above and recover an accurate representation of the undistorted surface topography. A simple, often used, method is 1-D line fitting, which is an effective way to get a rough representation of the surface by removing much of the scanner tilt and bow at the expense of some induced artifacts. Methods for the subtraction of 2-D distortions, which add less artifacts than does 1-D line fitting, are 2-D polynomial fits. 2-D polynomial fits can be performed by using either the whole figure for the fit, or only certain regions of the image determined by using thresholds. Figure 2 shows a comparison of different processing methods on a standard calibration grating as well as on a lipid bilayer of mixed composition. Figure 2A and Figure 2E show the starting data. If present, the 1-D errors must be corrected first because any attempts to correct for the 2-D distortions will be biased by these 1-D offsets. To do this, the offset caused by the 1-D distortions is removed from each line. Higher order 1-D operations, such as line fits or higher order polynomial subtraction, should not be used to correct the inherent 2-D distortions because they will destroy interline relationships in the data, and can generate false artifacts between lines [13,31,37], see Figure 2B and Figure 2F which show the results of 1-D second-order polynomial removal. This line-byline polynomial subtraction generates many artifacts in the data, for example the surface surrounding the pits appears raised along the fast scan axis in Figure 2B, and continuous levels in Figure 2F have offsets from line to line and are not perpendicular to the image plane. Figure 2C and Figure 2G show the removal of 1-D offsets followed by a second-order background polynomial removal. For the grid, there is a significant improvement in the representation of the surface. In Figure 2F, the sample topography is of the same order as the 2-D distortions. While 2-D operations are less prone to induce artifacts, performing a global 2-D polynomial fit and background subtraction leaves significant residual distortions ( Figure 2G). These distortions can be avoided by using a thresholded flattening, instead of a global flattening Figure 2D and Figure 2H (details discussed in Section Results and Discussion, "Algorithm description"). We conclude that a suitable way to process the images is to first determine and subtract the line-by-line offsets. Second, fit only the part of the image that is the flat background with a 2-D polynomial. Finally, subtract the calculated 2-D polynomial from the entire image. For an automated algorithm, the problem reduces to accurately distinguishing the flat background from the sample topography. In the rest of this paper, we describe a method for iteratively determining both the flat background, and the line-by-line offsets. Once these quantities are known, subtracting the line offsets and correcting the 2-D distortions can be performed with only two imageprocessing steps on the raw data.

Results and Discussion
1 Algorithm description Figure 3 shows a general flow diagram of the processing algorithm from data import to final output. The algorithm has three major blocks: (a) identify the background region, generate a mask of the background region and estimate the polynomial background; (b) determine 1-D background offsets from the raw data within the mask; and (c) subtract 1-D offsets from the raw data followed by a single masked background flattening. This final step ensures that a minimal number of modifications There are three main operational blocks used in the process. Each block starts with the raw data as input and generates unique outputs. The first major block applies a series of 1-D offset corrections and thresholded polynomial background subtractions to generate a masked image of the region to be flattened as well as an estimate of the polynomial background. The second block subtracts this estimated background from the raw data and runs thorough a series of median corrections and flattening steps only on the data within the masks. This sequence determines the unique 1-D offset for each line. The third block subtracts these 1-D offsets from the raw data and performs a single polynomial fitting on the data within the mask. This background is subtracted from the entire image (raw-offsets). Finally, the height of the background is shifted to a defined reference value to ensure continuity within the image sequence. All vertical scales in the figure are in nanometers.
are performed on the raw data to generate the final output. The final output is offset to consistent values and exported for later use. The results of the automated processing routine on an example lipid bilayer of mixed composition are shown as the inputs and outputs of each major block. Each block will be discussed separately in detail.
1.1 Identify the background, generate a mask, estimate the polynomial background The purpose of this section is to identify the background region in the image and estimate the initial polynomial background. In order to accomplish this, the algorithm does the following: (a) roughly corrects 1-D offsets present in the image (scars and offsets); (b) monitors the distribution of height values in the image to track the progress of each step, and rejects steps that broaden the distribution; (c) intelligently removes global tilt and scanner bow (often greater than sample topography); (d) adaptively identifies the background; (e) fits 2-D polynomials to the identified background; and (f) removes these fitted polynomials from the image.

Scar identification and median correction:
The first process is to import the data from the binary format and convert it into relevant units. The standard deviation of the image is calculated as a reference point. Figure 4A shows the data as captured. Sample tilt dominates the image.
A basic scar identification and correction step is used to minimize the influence of the scars on the correction of 1-D artifacts. Scars are identified by detecting sections within a line L n that are greater than surrounding lines by more than two times the standard deviation of the difference in medians.
( Equation 5). Then each data point identified within a scar is replaced with the median of five neighboring data points in the last line before the scar: diff(X) calculates differences between adjacent elements of X. If X is an array, then diff(X) returns an array, one element shorter than X, of differences between adjacent elements: [X 2 − X 1 , Second, a weighted median-difference correction corrects artifacts from the median correction caused by differences in topography in the image. This process is only performed on lines not previously identified as scars. The median difference correction calculates the median of the difference between two consecutive scan lines excluding outliers greater than two times the standard deviation of the differences between scan lines. This outlier exclusion reduces influence from edges within a scan line. The median difference value is subtracted from all the values within the second scan line, see Equation 7.
Third, an additional iterative refinement minimizes residual offsets generated from changes in topography. This subtracts a running average of the median difference between the scan line and up to the preceding k scan lines within a defined threshold. This does two things: one, it further excludes edge effects along topographic contours; and, two, it maintains continuity over long distances in the image and prevents the build up of highfrequency errors from scan line to scan line. This process runs iteratively in alternating directions. Currently, the algorithm uses four iterations with a maximum span, k, equal to 15 lines.
These steps handle discontinuities in the original image that do not represent actual topography. In cases where there is significant topographic variation parallel to the fast scan axis, i.e., a calibration structure, these corrections can actually introduce significant artifacts into the data and should not be used. In these cases, these 1-D corrections can be turned off in the algorithm. In Figure 4B, the initial median correction improves the image by removing both the 1-D artifacts as well as most of the tilt in the vertical direction. It is important to note, that in all these operations, there is always only a constant offset subtracted from all points in a line, and not a linear fit or higher polynomial. Such higher-order 1-D corrections can create artifacts that cannot be undone at a later stage.

2-D Tilt removal:
The third process is a basic polynomial surface fit of either first or second order and subsequent removal to facilitate the thresholding process by removing enough distortion to make it possible to separate topographic features from the background. Primarily, this removes sample tilt, which has a significant contribution to the variation in the observed heights. In cases of very flat samples where the changes of the topography are on the same order or smaller than the sample tilt or scanner bow, removal of a second-order surface generally improves the separation of the topography from the background. The algorithm does a separate first-order and second-order removal process and calculates the standard deviation of both results. The standard deviations are compared to the standard deviation after the median correction. The image with the lowest standard deviation progresses to the next step. Figure 4C shows the results of the first, unthresholded background correction. A second-order polynomial removal improved the image the most. After this step, the two distinct levels in the system manifest in both the image and the corresponding histogram. Some residual tilt remains in the image, as the bottom is somewhat higher than the top.

Thresholded flattening:
The fourth process is an adaptive image thresholding algorithm and background removal. In order to calculate the appropriate threshold for the 2-D background removal, the algorithm identifies discrete peaks in the histogram corresponding to different topographic levels in the image. In the absence of discrete levels, only one general peak is used. The histogram is fit to a summation of m ≤ 8 Gaussian functions centered at each peak, z i , where i ranges from 1 to m.
In the case that the number of levels is greater than m, only m peaks with the greatest area will be fit. From the Gaussian fit, the standard deviations, σ i , and areas covered, A i , are calculated. Nominally, the threshold for the 2-D surface fitting is z j ± nσ j , where n is a user-defined range and j is the peak used for the 2-D surface fit. (The algorithm default is to use the highest peak found in the histogram. Based on prior knowledge of the sample, the user can override this default and specify which peak to fit.) Since this span can include contributions from other peaks, the actual span can be narrowed to the intersection of the jth peak with the next closest peak. All the data within the threshold is fit to a polynomial surface of a given order (up to fifth). The extrapolated background is subtracted from the whole image.

Improvement check:
In order to test if a particular step has improved the image, the peak detection and Gaussian fit is repeated on a histogram of the corrected data. While the thresholded flattening only operates on a single level, the improvement check monitors the progress of all the fit peaks. This allows the algorithm to follow the progress of a sample with many levels. The criteria for this conditional check are: 1. The total area of the fits must not decrease by more than a small tolerance with respect to the image before the processing step. 2. A weighted average, , of the standard deviations is computed both before and after. This term must not increase by more than a small threshold.
If either of these conditions fails, the algorithm rejects the process and passes the previous image to the next step. In the case of the initial background removal the whole image is used. After the initial background removal, only regions identified as independent levels are considered. The user defines the number of levels present in the system. In the case of the lipids or grids, both independent levels are considered; whereas, in the case of a random sample on a substrate, only the substrate would be considered. Figure 4D shows the results of the first thresholded processing. This results in a significant improvement of the overall image flatness, shown in both the image and histogram. The algorithm based the thresholding on the higher level since it had the greatest area.

Determine 1-D background offsets
The purpose of this section is to accurately determine the 1-D offsets that exist from line to line. While the operations used above in Section 1.1 do a reasonable job, they are still prone to influence from sample topography. In order to minimize the influence of topography the algorithm performs its corrections only on the masked background region determined in Section 1.1. First, the estimated background from Section 1.1 is subtracted from the raw data in order to generate a reasonable starting point for the first masked median correction ( Figure 5A). The masked median correction subtracts the median of all the data within the mask on a given line, from that line ( Figure 5B).

(9)
The second step is an additional polynomial flattening of the masked data in order to ensure a flat background ( Figure 5C). This step is followed by an additional masked median correction ( Figure 5D). Following each step, a check is performed to ensure that the step decreased the standard deviation of the data Figure 5: This figure shows the structure used to determine the 1-D offsets present in the image. First, the estimated background from the prior block is subtracted from the raw image. Second, a median correction is performed on the data within the mask. Third, an additional polynomial flattening on the masked region reduces the influence from the topography. Fourth, a second masked median correction refines the calculated offsets. After each step, except the initial background subtraction, an improvement check on the masked region verifies that the step improved the image. The final offsets are determined and passed through to the next section. The vertical scale of all images and horizontal scales of all histograms are in nanometers.
within the mask. Finally, the cumulative offsets are calculated for each line. In the case where a line has no pixels contained in the mask, the algorithm calculates the offset from the final output of Section 1.1 instead. These final outputs are passed on to the next section.

Subtract offsets, masked flattening, reference shift data
This final section of the algorithm subtracts the offsets calculated in Section 1.2 from the raw data ( Figure 6A). Next, the data within the mask calculated in Section 1.1 is fit to a single 2-D polynomial. The resulting polynomial is subtracted from all the data ( Figure 6B). A final histogram is computed and fit with Gaussians. The difference between the defined background reference height and the fitted center of the background peak in the final histogram is added to all the data. This final data is exported for future use.
2 Processing high-speed dynamic images in biological systems

HS-specific artifacts
Operating AFMs at high speeds has the potential to generate new distortions not normally observed at conventional speeds. Such artifacts include structural scanner resonances coupling into the sample topography. These sorts of spurious resonances can be observed in piezo tube scanners as a turn-around ripple in the images (damped oscillations of a given frequency). This stems from excitations in the lateral direction coupling into the vertical motion. Because the scanner used in these experiments has a relatively low structural resonance, this turn-around can be excited to a substantial extent at even moderate scan speeds. Given that the amplitude of the turnaround ripple in tube scanners like the one used in these experiments can easily be greater than the topography of the system, the turn-around ripple must be dealt with in some fashion. The best way is to avoid it altogether through either input shaping [38][39][40][41] of the drive signals or through electrical damping of the resonances [42]. For our experiments, we use a self optimizing method that determines the scanner resonances and compensates them with an input The fully automated processing using the lower level for the background calculation is shown in the middle row. In the bottom row, an additional vertical median correction has been used. This corrects some of the residual artifacts from the turnaround ripple of the scanner, which are not corrected with the resonance compensator described by Burns et al. [38]. Each image is a 4 µm image with 256 lines captured at 78 lines/second while tapping in fluid. These imaging parameters give an approximate acquisition of 3 seconds per image. Each frame in this figure represents every third image in the overall sequence, approximately 9 seconds apart. (Full image sequence as Figure S1 in Supporting Information File 1 or in movie format as Supporting Information File 2). In the unthresholded sequence, many artifacts remain, making it difficult to track the progression. In the fully automated sequence, the automated algorithm does a much better job of correcting the data and outputting consistent images for comparison.
shaper [38]. Using this resonance-compensator system, the amplitude of the turn-around ripple in the image can be nearly eliminated in a fashion that is sample independent and completely transparent to the rest of the experiment. While this method works very well, some residual, subnanometer distortions can remain.

Vertical median correction
While model-based filters do a very good job of reducing the amplitude of the turn-around ripple, they may not completely eliminate it. This is especially problematic on very flat samples with nanometer scale topography, where a subnanometer distortion is a significant part of the overall topography, such as in the images of lipid membranes shown in Figure 7. This residual turn-around ripple can be seen on the right-hand side of the images in the middle row. The turn-around ripple appears on the right side of the image because retrace images are being shown. This residual error can be corrected by using an additional vertical median correction on the masked data after the final background subtraction in Section 1.3 and before the final offset. This vertical median correction is exactly the same as Equation 9, but it runs in the perpendicular direction. This correction is only applicable to HS-AFM data that exhibits some turn-around ripple and should not be used otherwise. Figure 7 demonstrates the utility of the automated processing algorithm on a dynamic model biological system. The sample consists of two phase-segregated lipids, DLPC (blue) and DPPC (gold). The image sequence shows the rapid degradation of the DPPC domains following rapid heating from interactions with the AFM tip. The sequence shows every third image. The upper sequence shows the results of standard unthresholded flattening and the lower sequence shows the results of the fully automated processing. Both sequences use the same polynomial background orders. The fully automated processing does a much better job of identifying the background and flattening it correctly. This makes it much easier to follow the changes in the system. Moreover, this automated flattening process allowed for a significant improvement in the speed of processing without sacrificing the quality of the results. The images shown in Figure 7 and Figure S1 in Supporting Information File 1 are only a few sample images from a much longer image sequence many hundreds of images long. The entire sequence was processed in only a few hours, representing a significant time saving compared to hand-correcting each image.

Conclusion
We have developed an automated, adaptive image-processing algorithm for high-speed AFM image sequences. This is achieved by identifying the background region and determining line-by-line offsets using an iterative process. The output of this iterative process is used to perform a single line-by-line offset correction, followed by a single 2-D polynomial-background removal step. This order of operations ensures that a minimal amount of manipulation is performed on the raw data to generate the final output. During the iterative background-identification process, the algorithm uses the standard deviation of one or more peaks in the height histogram of the image as a metric for determining the accuracy of the thresholded background identification. A similar metric is used to determine the accuracy of the line-by-line offset correction. The algorithm is specifically tailored to image sequences, ensuring consistent processing, offsets and contrast settings between frames in AFM movies in a short time period. On a standard PC, each channel takes about 10 seconds to process. This algorithm can be used on AFM data from most modes. We have demonstrated its applicability to data acquired in tapping mode in air, tapping mode in fluid, and quantitative nanomechanical mapping (QNM) in fluid. Finally, this algorithm can also be applied to strongly stepped samples, such as atomic layers. Figure S1 in Supporting Information File 1 shows the correction of samples with four and eight distinct levels. In both cases, the error is less than the noise in the image. With the release of modern HS-AFM systems by commercial manufacturers, this sort of automated processing will provide a significant benefit to this emerging research area in surface science and related fields. The algorithm and an associated graphical user interface are available at http://lbni.epfl.ch under the software section.

Experimental Lipid preparation
Small unilamelar vesicle mixtures of 1,2-dilauroyl-sn-glycero-3-phosphocholine (DLPC) and 1,2-dipalmitoyl-sn-glycero-3phosphocholine (DPPC) were prepared by sonication. Both lipid types were purchased from Avanti Polar Lipids Incorporated (Alabaster, AL, USA). Lipid powders were mixed prior to vesicle formation at a nominal molar ratio of 1:2, DLPC/DPPC. Vesicle solutions at 1 mg/mL were formed by transferring an appropriate mass of lipid into glass vials and dissolved with chloroform. The chloroform was evaporated off with dry nitrogen gas, leaving a thin film on the glass vial. The film was hydrated with Milli-Q water (Millipore, Billerica, MA, USA), generating large multilaminar vesicles (LMVs). The LMVs were then sonicated with a probe sonicator (BioLogics Inc, Manassas, VA, USA) to generate small unilaminar vesicles (SUVs). The SUVs were centrifuged to remove metal particles left from the probe sonicator. A 35 µL amount of the lipid preparation was warmed to 37 °C and deposited onto freshly cleaved mica surfaces, forming bilayers through vesicle fusion. Surfaces were allowed to incubate for at least a half hour in a humid environment at room temperature.

AFM imaging in air
Images were captured on a Multimode system with an E-scanner (Bruker Nano: Santa Barbara, CA, USA). Standard TESPA (Bruker AFM Probes: Camarillo, CA, USA) tapping cantilevers were used. The cantilever was driven at 348.44 kHz with an amplitude of 20.14 mV. Images at a size of 14.25 µm were captured at 512 × 512 pixels with a line rate of 1.5 lines/ second. The data shown in the paper are crops from the center of the acquired image.

Quantitative Nanomechanical Mapping (QNM) -AFM imaging in fluid
Images were captured on a Multimode system with an E-scanner (Bruker Nano: Santa Barbara, CA, USA). A standard DNP-A (Bruker AFM Probes: Camarillo, CA, USA) cantilever was used with a spring constant of 0.40 N/m. Images at 5 µm were captured at 512 × 512 pixels with a line rate of 1 Hz. Manual control of the QNM parameters was used to minimize the applied force on the sample and the QNM drive amplitude. Each channel's acquired limits were minimized to limit bit quantization in the DSP.

HS-AFM imaging in fluid
Images were captured on a modified Multimode system with an E-scanner (Bruker Nano: Santa Barbara, CA, USA). A customized small-lever head allowed for the use of small cantilevers (SCL-Sensor.Tech., Vienna, Austria). The cantilever had a resonance frequency in fluid of 266.49 kHz, a spring constant of 0.54 N/m and a Q value of 2.68. Square areas of 4 µm were scanned at 78 lines/second at 256 × 256 pixels giving an approximate image acquisition time of 3 seconds/ image. The fast-scan drive signal was passed through a custom filter designed to minimize the excitation of the tube scanner resonance [38]. The internal PID feedback of the Nanoscope 5 controller was bypassed with a Labview controlled FPGA based PID with a loop rate of 575 kHz.

Supporting Information File 2
Image sequence movie The movie shows the entire image sequence of the fully corrected data, with vertical median correction, from


extfieldsofstudy:
Mathematics
Computer Science
Medicine

text:
WaVPeak: picking NMR peaks through wavelet-based smoothing and volume-based filtering

Motivation: Nuclear magnetic resonance (NMR) has been widely used as a powerful tool to determine the 3D structures of proteins in vivo. However, the post-spectra processing stage of NMR structure determination usually involves a tremendous amount of time and expert knowledge, which includes peak picking, chemical shift assignment and structure calculation steps. Detecting accurate peaks from the NMR spectra is a prerequisite for all following steps, and thus remains a key problem in automatic NMR structure determination. Results: We introduce WaVPeak, a fully automatic peak detection method. WaVPeak first smoothes the given NMR spectrum by wavelets. The peaks are then identified as the local maxima. The false positive peaks are filtered out efficiently by considering the volume of the peaks. WaVPeak has two major advantages over the state-of-the-art peak-picking methods. First, through wavelet-based smoothing, WaVPeak does not eliminate any data point in the spectra. Therefore, WaVPeak is able to detect weak peaks that are embedded in the noise level. NMR spectroscopists need the most help isolating these weak peaks. Second, WaVPeak estimates the volume of the peaks to filter the false positives. This is more reliable than intensity-based filters that are widely used in existing methods. We evaluate the performance of WaVPeak on the benchmark set proposed by PICKY (Alipanahi et al., 2009), one of the most accurate methods in the literature. The dataset comprises 32 2D and 3D spectra from eight different proteins. Experimental results demonstrate that WaVPeak achieves an average of 96%, 91%, 88%, 76% and 85% recall on 15N-HSQC, HNCO, HNCA, HNCACB and CBCA(CO)NH, respectively. When the same number of peaks are considered, WaVPeak significantly outperforms PICKY. Availability: WaVPeak is an open source program. The source code and two test spectra of WaVPeak are available at http://faculty.kaust.edu.sa/sites/xingao/Pages/Publications.aspx. The online server is under construction. Contact: statliuzhi@xmu.edu.cn; ahmed.abbas@kaust.edu.sa; majing@ust.hk; xin.gao@kaust.edu.sa

Among the existing automatic methods, AUTOPSY (Koradi et al., 1998) and PICKY (Alipanahi et al., 2009) are the most accurate. Both of these methods attempt to estimate the noise level of the given spectrum. The noise is assumed to be white Gaussian noise and is estimated within small regions of the spectrum. All the data points that have lower intensities than the estimated noise level are eliminated, which eliminates most parts of the spectrum. The remaining spectrum looks like a set of disconnected components. Instead of dealing with all the components altogether, both AUTOPSY and PICKY extract separate components. AUTOPSY extracts components by a simple flood fill algorithm, whereas PICKY further subdivides weakly connected components into smaller ones. After the components are formed, both methods conduct peak picking within each component. AUTOPSY first identifies strong, obvious peaks from components. Each of these components is then decomposed as the outer product of 1D lineshapes. The lineshapes are then clustered and used to detect overlapping peaks from the remaining components. Different from AUTOPSY, which uses symmetry properties of peaks, PICKY directly applies SVD on each of components. Since the components formed by PICKY are much smaller and simpler than the ones formed by AUTOPSY, rank-one SVD seems to be powerful enough to identify peaks from the components. However, both AUTOPSY and PICKY have a quite high false positive rate. Therefore, a refinement step is performed in each method to reduce the number of false peaks. In AUTOPSY, the integration, symmetrization and filtering modules are applied. In PICKY, the peaks are first sorted according to the intensities. A certain number of the strongest peaks (usually 1.2K, where K is the expected number of peaks) is kept. The peaks from different NMR spectra, which share common atoms, are used to cross-eliminate false positive peaks.
Although AUTOPSY and PICKY have demonstrated impressive accuracy on different benchmark spectra, they both have two bottlenecks. The first is that all the data points with low intensities are eliminated in both methods. However, a number of true peaks actually have low intensities due to various reasons, including the sensitivity of the NMR spectrometer, the strength of the magnetic field, the characteristics of the target proteins and the local dynamics of the spin systems. Unfortunately, spectroscopists need the most help to deal with these weak peaks. An ideal peak-picking method should therefore be able to identify the weak peaks instead of disregarding them.
The second bottleneck is the way the two methods eliminate the false positive peaks. Apparently, a brute force method that selects all the local maxima as peaks should have higher sensitivity than any other method. The issue is that since the real spectra are very noisy, the brute force method will generate a huge number of peaks, making it almost impossible for users to identify the true ones. On the other hand, all peak-picking methods try to smooth the spectra such that there is a good tradeoff between the sensitivity and the specificity. For example, PICKY first eliminates all the weak data points and then smoothes the components by rank-one SVD. The problem now is how to rank the peak candidates in an order such that the true peaks are at the top. The existing methods use either intensitybased ranking or cross-references from other spectra. However, the intensity of a single data point is not informative enough to distinguish a true peak from a false one. Cross-references from other less reliable spectra, on the other hand, may eliminate some true peaks.
In this article, we introduce WaVPeak, a fully automatic peakpicking method that is based on wavelet-based smoothing and volume-based filtering and overcomes the two bottlenecks of the existing methods. We first propose the use of the Daubechies 3 wavelet for the peak-picking problem. We suggest and demonstrate that this wavelet is the most suitable one for this problem because it can smooth the spectrum, sharpen the peak shapes and maintain the peak locations. In this way, no data point will be eliminated. Furthermore, the weak peaks that are embedded in the noise will be recovered and will become visible. Then, a brute force method is applied to isolate all the local maxima in the wavelet-smoothed spectrum. However, the number of local maxima in the waveletsmoothed spectrum is much smaller than the number in the original spectrum. The peak candidates are then sorted according to the estimated volume of the peak shapes, which has a much stronger distinguishable power than intensity-based sorting has.

METHODS
WaVPeak consists of three main steps, i.e. wavelet-based smoothing, brute force peak picking and volume-based filtering. Given any spectrum, it is first smoothed by wavelets. A brute force algorithm is then applied to identify all the local maxima as the initial peak candidates. The initial peaks are ranked according to their estimated volume. The details of the wavelet-based smoothing and volume-based filtering are discussed in this section.

Wavelet-based smoothing
Wavelets are mathematical functions that cut up data into different frequency components. Subsequently, each component is studied with a resolution matched to its scale. Wavelets have advantages over traditional Fourier methods in analyzing physical situations in which the signal contains discontinuities and sharp spikes. Interestingly, wavelets were developed independently in the fields of mathematics, quantum physics, electrical engineering and seismic geology.
Despite such applications of wavelets in NMR signal processing, waveletbased automatic peak picking has not been the focus of research attention. The most relevant work is Dancea and Güntert, 2005. Dancea and Güntert combined different wavelet procedures together in a consensus manner to de-noise the 15 N-NOESY spectrum in order to generate a 3D structure for the Sud protein from Wolinella succinogenes. They considered the consensus results from three wavelet-base functions, i.e. Symmlet, Daubechies (Daubechies, 1992) and Coiflet. Although the structure was accurately determined, the proposed multi-stage de-noising technique was specific for one particular spectrum and could not be easily generalized.
De-noising or smoothing by wavelets is an indispensable step for automatic NMR peak picking due to the various sources of noise in the spectra. It is an important step as the de-noising is carried out without smoothing the sharp structures. The result is a cleaned-up signal that still captures important details. The de-noising consists of three steps: (1) transforming the spectra to the wavelet domain using some appropriate 2D wavelets (more discussion below); (2) applying thresholding methods (e.g. hard or soft thresholding), i.e.
setting all coefficients to zero that are less than a particular threshold (we adopt soft thresholding later in this article); and (3) inverse-transforming the thresholded coefficients to reconstruct the data set in the signal domain.
We now elaborate more on 2D wavelets in Step 1. First, 2D wavelets can be constructed from 1D wavelets (see Mallat, 1989). Given a 1D scale function, φ(x), and its corresponding wavelets, ψ(x), we define the 2D scale function as and the corresponding 2D wavelets as The critical issue then becomes how to choose the 1D wavelets φ(x) and ψ(x) as different wavelets serve different purposes. Since our objective is to identify peaks, the ideal wavelet function should resemble the shapes of true peaks, i.e. a cone shape in our NMR peak-picking problem. We have explored various wavelet families with different parameters on MATLAB, including Daubechies, Symlets, Coiflets and Biorthogonal. It turns out that Daubechies 3, Symlets 3 and Biorthogonal 2.4 significantly outperform others on the peak-picking accuracy, when evaluated with real NMR spectra. The accuracy of these three wavelets is very similar with Daubechies 3 having slightly higher sensitivity. Therefore, in WaVPeak, we use 2D Daubechies 3 as the default wavelet. Figure 1 illustrates the scaling function, φ, and the wavelet function, ψ, of the 1D Daubechies 3 wavelet. Clearly, both functions produce a cone shape, which looks like a peak shape after smoothing. This can be further confirmed from Figure 2, which shows the original spectrum of the 15 N-HSQC of the VRAR protein with wavelet-based smoothing using Daubechies 3. Clearly, the original spectrum is much noisier than the wavelet-smoothed spectrum. Furthermore, the peak shapes are smoothed and sharpened in the smoothed spectrum as well.

Volume-based filtering
After the initial peaks are identified by the brute force algorithm on the wavelet-smoothed spectrum, we need to identify the true peaks from the relatively large set of initial peaks. The idea is to rank the initial peaks according to a certain criterion, so that all the true peaks are ranked as highly as possible. Traditional methods use the intensity of the peak location as the filtering criterion. However, the intensity of one single data point can be badly biased by various sources of noise in the NMR spectrum; intensity is thus not a reliable criterion.
Here, we propose the use of the volume of the peak instead of the intensity. However, the number of data points that belong to a certain peak shape is not apparent. We therefore need to select a region around the peak location in The original spectrum The wavelet-smoothed spectrum order to estimate the volume. A trivial solution is to set a predefined window size, such as three. This region includes all the data points that differ by at most one in all the dimensions. To estimate the volume of the peak shape in this region, we assume that the region of the entire peak shape is at least larger than this small region with a window size of three. The volume can thus be estimated as the sum of the volume of the 'pillars' (in this case, nine), each of which corresponds to a part of the peak shape over a grid. The volume of a pillar can then be simply estimated as the area of the grid multiplied by the intensity of the data points corresponding to that grid. Note that in a given spectrum, the area of the grid is a fixed value. Therefore, we can use the sum of all the intensities in this region to estimate the volume of the peak.
However, in real spectra, the span of a peak shape over different dimensions can have different numbers of data points. Figure 3 shows an example of a peak shape on a 2D spectrum. This peak shape spans over three and five data points on the two dimensions, respectively. Thus, taking a fixed region of size 3×3 does not give a good estimation of the volume of the peak. In order to solve this problem, we propose a self-adapted, window-based volume estimation.
As shown in Figure 3, the ideal size of the window is 3×5, which includes all the data points that differ by at most one in one dimension and by at most two in the other dimension. From Figure 3, it is easy to show that: where r x and r y are the levels of neighbors being considered for dimensions x and y, respectively, and I p , I x , I y are the intensities for the peak, the direct WaVPeak Fig. 3. A peak shape with different spans over the x-and y-dimensions.

Input:
The wavelet-smoothed NMR spectrum in d-dimensional space; Input: n initial peak candidates (p 11 ,...,p 1d ),..., (p n1 ,...,p nd ). Input: The length of the protein, L. Output: The re-ranked list of the initial peaks. Set level size to be l 1 = l 2 = ··· = l d = 1. for t = 1,...,n do Calculate the volume as end for Sort the initial peaks according to the volume V and take the top L.
Calculate the average of the intensity difference between the peak and its direct neighbor for each dimension, end for Sort the initial peaks according to the volume V and return the top L.
neighbor of this peak on x-dimension, and the direct neighbor of this peak on the y-dimension, respectively. Consequently, we have To adjust the size of the region in which we calculate the volume, we set the level size that corresponds to the dimension with the smaller r value to be one. Without loss of generality, assume r x < r y . The level size for the y-dimension is set to [

Ip−Ix
Ip−Iy ]. The generalization to higher dimensional spectra is straightforward. Algorithm 1 gives the pseudocode of the selfadapted, window-based filtering algorithm for an arbitrary spectrum.

Evaluation criteria
In order to evaluate a peak-picking method objectively, we adopt three criteria, i.e. recall, precision and F-score. Let TP denote the number of true peaks that are discovered by the method, FP the number of false peaks returned by the method and FN the number of true peaks that are not discovered by the method. Recall is defined as TP/(TP +FN), where TP +FN is the total number of ideal peaks. Similarly, precision is defined as TP/(TP +FP), where TP +FP is the total number of peaks returned by the method.
Recall measures the ability of the method to discover the true peaks, whereas precision measures the ability to reject false peaks. In the peakpicking problem, recall is more important than precision because a false peak can be possibly eliminated in the following NMR data analysis process.
We further apply the F-score to measure the tradeoff between recall and precision, where the F-score is defined as the harmonic mean of recall and precision, i.e. 2·Recall ·Precision/(Recall +Precision).

Data set
To fairly evaluate the performance of WaVPeak, we test WaVPeak on the spectra set proposed in Alipanahi et al. (2009), which is one of the largest benchmark sets on the peak-picking problem. The set contains 32 spectra extracted from 15 N-HSQC, HNCO, HNCA, HNCACB and CBCA(CO)NH from eight proteins, i.e. TM1112, YST0336, RP3384, ATC1776, CASKIN, HACS1, VRAR and COILIN.

Performance on the benchmark set
WaVPeak is applied on each spectrum of the benchmark set to automatically pick peaks. Although different spectra in the data set share the common 15 N and 1 H atoms, the cross-referencing is not applied in our experiments so that we may objectively test the peak-picking accuracy of the proposed method.
The Daubechies 3 wavelet is employed to smooth the original spectrum. Figure 2 shows one example of a wavelet-smoothed 15 N-HSQC spectrum of the VRAR protein. The smoothed spectrum is much smoother than the original spectrum. Furthermore, the property of the Daubechies 3 wavelet ensures that the peak shapes are more obvious in the smoothed spectrum. A brute force algorithm is then used to select all the local maxima in the smoothed spectrum to build the initial peak list. Note that due to the smoothness of the wavelet-smoothed spectrum, the number of local maxima becomes an order of magnitude smaller than the number in the original spectrum.
All initial peaks are then ranked according to the volume-based estimation. In the first round, the default window size (differing by at most one step size in each dimension) is used to define the neighbors. Then, the top N peaks are extracted, where N is the length of the target protein, which is given as the input. The window size is adjusted according to Algorithm 1. The initial peaks are then re-ranked according to the estimated volume in the updated neighborhood. Finally, the top M peaks are selected as the prediction results, where M is either given by the users or is set to 1.2K as the default (where K is the expected number of peaks in this spectrum).
We compared the performance of the publicly available version of PICKY with WaVPeak, with a default noise cutoff threshold of 5 according to Alipanahi et al. (2009). Table 1 lists the performance of PICKY and WaVPeak on the benchmark set when the top 1.2K peaks are considered. It can be seen that when the same number of peaks is considered, WaVPeak has a consistently better recall than PICKY on four of the five types of spectra. It is well acknowledged in the NMR community that 15 N-HSQC is the most reliable and clean spectrum among the five, whereas the HNCACB and CBCA(CO)NH spectra are the noisier. WaVPeak has an average improvement of 3%, 2%, 16% and 13% on recall over PICKY on 15 N-HSQC, HNCO, HNCACB and CBCA(CO)NH, respectively. On the other hand, the precision of WaVPeak is comparable to that of PICKY on 15 N-HSQC, HNCO and HNCA. The improvement in precision is significant for CBCA(CO)NH. However, PICKY has a much higher precision than WaVPeak on HNCACB. Since the number of top peaks is the same for both methods, this higher precision of PICKY seems to conflict with the lower recall. This is due to the fact that PICKY detects fewer than 1.2K peaks in the three HNCACB spectra. Figure 4a-e show the receiver operating characteristic (ROC) curves of PICKY and WaVPeak when different numbers of top peaks are considered. Therefore, WaVPeak is consistently and significantly more sensitive than PICKY, and the overall performance (F-score in Table 1 and area under curves (AUCs) in the Figure  4a-e is also better than that of PICKY.
It is clear from Table 1 that PICKY has a higher recall than WaVPeak on only 3 of the 32 spectra. It should be noted that on eight spectra, WaVPeak is able to identify significantly more true peaks than PICKY identifies (the improvement in recall is over 15%). These spectra are the 15 N-HSQC spectrum of CASKIN, the HNCACB spectra of CASKIN, VRAR and COILIN, and the CBCA(CO)NH spectra of RP3384, CASKIN, VRAR and COILIN. With the exception of the CBCA(CO)NH spectrum of RP3384, PICKY's noise level elimination step overestimates the noise level. The hardthreshold-based elimination step then removes many true peaks. Some of these true peaks still have relatively strong intensities, but these intensities are lower than the estimated noise level. The others are weak peaks that are completely embedded in the noise, which are impossible for any hardthreshold-based method to discover. Therefore, PICKY has very low recall on these seven spectra, but quite high precision due to the small number of total peaks returned. For CBCA(CO)NH of RP3384, PICKY returns a large number of peaks. However, the recall is still 20% lower than that of WaVPeak. This further confirms that the wavelet-smoothed spectrum can discover the baseline peak shapes of the weak peaks that are embedded in the noise level, and the volume-based filtering can select the peaks, though they have very low intensities.
In Alipanahi et al. (2009), the cross-reference step is used to refine the peak lists of different spectra that share common atoms. This, apparently, is an efficient way to remove the false positives in practice. However, this assumes that a set of spectra are available. Our focus is to discover a more powerful and general method for peak picking for any given spectrum. Therefore, WaVPeak does not have a default cross-referencing step. It should be noted that cross-referencing can be easily implemented and applied as a post-processing step for any peak-picking method. When compared with the accuracy of PICKY after cross-referencing, WaVPeak without crossreferencing still has comparable F-scores (87% versus 87%, 82% versus 83%, 81% versus 80%, 72% versus 70% and 78% versus 78% for the five types of spectra, respectively). Table 1 demonstrates the performance of WaVPeak when the top 1.2K peaks are considered. It is also necessary to see how the accuracy changes when different number of peaks are considered. When the accuracy of WaVPeak and PICKY is compared using different numbers of peaks, the results are very similar to those presented in Table 1. We further evaluate the contribution of wavelet-based smoothing and volume-based filtering. To measure the contribution of wavelet-based smoothing, we compare the smoothing module of WaVPeak plus the traditional intensity-based filtering by PICKY, which includes intensitybased filtering as a default. To evaluate the contribution of volume-based filtering, we compare the smoothing module of WaVPeak with different filtering methods, including intensity-based filtering and self-adapted, window-based filtering. Figure 4f shows the recall curve of the three methods when different numbers of peaks are considered in the CBCA(CO)NH spectrum of the ATC1776 protein. The expected number of peaks is 180. Therefore, when <180 top peaks are considered, the recall for all methods is relatively low. The three methods have the same recall when <80 top peaks are taken into account. When we consider >80 peaks, wavelet-based smoothing plus the self-adapted, window-based filtering, i.e. WaVPeak, consistently has highest recall, whereas wavelet-based smoothing plus intensity-based filtering has the second highest recall. Note that the list of initial peaks for the waveletsmoothed spectrum is the same for the two different filtering methods. The filtering method only re-ranks the peaks according to different measures. This clearly demonstrates that volume has much stronger power to distinguish peaks than intensity has. When using intensity-based filtering, wavelet-based smoothing has slightly better recall than PICKY. This demonstrates the power of the wavelet to discover the weak peaks that are discarded in PICKY's noise level elimination. By combining wavelet-based smoothing and self-adapted windows, WaVPeak has significantly higher recall than PICKY has. When triple the expected number of peaks are considered, WaVPeak has 95% recall, whereas the recall of PICKY is 89%.

Implementation details
One other advantage of WaVPeak is that it can be easily implemented. WaVPeak does not have as many steps as the other methods. WaVPeak is also an open source program. The MATLAB version of the source code with two testing spectra are available at http://faculty.kaust.edu.sa/sites/xingao/Pages/Publications.aspx.
We are currently developing a web server and making WaVPeak a plug-in for interactive NMR data processing tools, such as SPARKY (Goddard and Kneller, 2007).

DISCUSSION
WaVPeak is proposed as a general peak-picking approach for any NMR spectrum. Therefore, WaVPeak does not use the properties of the specific spectra or the specific amino acids. Such properties, of course, will be useful to improve the accuracy of WaVPeak. One efficient way to reduce the number of false positive peaks further is to take the consensus of the peak lists from different peak-picking methods, such as PICKY and WaVPeak. This can significantly increase the precision (data not shown) while maintaining comparable recall. We have made WaVPeak an open source program so that expert knowledge can be easily encoded. Furthermore, we will make WaVPeak a plug-in to SPARKY or NMRView so that users can do interactive peak picking. WaVPeak estimates the volume within self-adapted rectangular regions. The rectangular regions, however, may not be ideal, especially for peaks with a lot of overlap. In such a case, the volume of one peak may be added to that of another overlapped peak. A possible solution is to use a V-N plot to identify true peaks. For each of the initial peak candidates, we can start with the smallest diamond-shaped neighborhood region (differing by at most one step size in at most one dimension), and gradually increase the size of the region. This will generate a series of alternating diamondshaped and rectangular regions. Consequently, we will have a series of estimated volumes for these regions. We can then plot the volume in relations to each neighborhood, which will be referred to as the V-N plot. The V-N plot (or its smoothed version) may reveal a lot of information about the patterns of the true peaks. For instance, if two peaks overlap a lot with each other or are even completely overlapping, it is possible that we only observe one peak in the wavelet-smoothed spectrum. However, in the V-Nplot, we can see that the AUC is almost twice the area of the other curves. If a local maximum is a fake peak, it is very unlikely that it has a regular and smooth peak shape. In the V-N plot, the curve will not be smooth either. The gradient change of the curve might be informative in identifying false peaks.
In evaluating the performance on the benchmark set in Section 3, the top M peaks are selected as the prediction results, where M is  either given by the user or is set to 1.2K as the default. In reality, it would be desirable to have a fully automatic data-driven method to identify true peaks with high recall rates among the candidate peaks obtained after wavelet-based smoothing. This is, in fact, a multiple testing problem, which has received much attention in the statistical literature since the seminal article by Benjamini and Hochberg (1995). Work is currently underway to clarify how it can help in resolve our peak-picking problem.

CONCLUSION
In this article, we introduced WaVPeak, an automatic peak-picking method that is based on wavelet-based smoothing and volume-based filtering. Experimental results demonstrate that the wavelet-based smoothing outperforms existing decomposition techniques for this task, and volume-based filtering outperforms traditional intensitybased filtering. The combination of these two ideas results in a novel method that is significantly more accurate than the state-of-the-art peak-picking methods.
grateful to Ming Li for making PICKY publicly available. We thank Virginia Unkefer for proof-reading the manuscript.