Our proposed method has four stages: 1) Image Calibration and Wave Extraction, 2) Signal Masking, 3) Global Analysis, and 4) Local Analysis (Figure 1). Raw images captured by the CCD camera are passed to the Image Calibration and Wave Extraction stage. An optional intermediate signal masking process could be applied to speed up the analysis and to get more robust measurements. The results (with or without masking) are then analyzed either by a global or a local approach, leading to the final outputs: propagation velocity, wave decay, and wave amplitude.
Image calibration and wave extraction
The purpose of this stage is to calibrate the data and isolate the Ca2+ that is released to the media using the raw images as input. The outputs are the basal conditions, statistical information of the set, and a sequence of images showing only the Ca2+ above the basal values. This change of concentration is associated with the passing wave, and hence it is important to accurately calibrate the whole set of frames, and to reduce any possible source of errors. To accomplish this, the user has to identify the frame where the wave propagation begins, and our software processes the frames according to the following steps:
Image registration
Since the Ca2+ wave is induced by a mechanical stimulation of a cell, minor motion may be introduced to the sample, and hence physical locations do not match exactly in all the frames. To correct that problem an automatic rigid alignment of the frames was implemented. We estimate the displacement between frames computing the Fourier Transform of each image, and calculating the difference of phase angles with respect to a reference frame. Since the images are captured with different filters, they do not have the same intensity mapping. A wavelet based pre-filtering is applied to the input data so that to obtain images with comparable intensities and reduced noise. This pre-filtering process consists on the following steps: first, the image is decomposed into wavelet layers, using the starlet transform [6] (formerly known as A Trous wavelets transform). Second, all wavelet layers are discarded, except for that corresponding to structures of 1px or 2px (user defined). Finally, this remaining layer is max-min rescaled to a normalized range.
Figure 2 shows an example of this procedure. Figure a) shows the difference between a reference frame captured at 340 nm and a target frame at 380 nm. In this case the misregistration is about 3px in the vertical axis. Figure b) shows the difference between the rescaled layers of 1px of the same images. Figure c) is the new difference of frames, after registration with the Fourier phase correlation.
Calcium concentration calculation
The calcium concentration can be obtained after computing the ratio between the 340 and 380 nm filtered images. However, to estimate precisely such concentration an additional calibration function must be applied [7]. Since the 340:380 ratio shows an approximately linear relation with calcium concentration and we were focused in characterizing the wave propagation (and not absolute calcium concentrations), we performed our data calculations from the computed ratio without such calibration.
Marker decay correction
The fluorescent marker suffers a natural decay over time, following closely an exponential function [8]. If the decay time were similar for 340 and 380 nm, computing the ratio of both filtered images would correct that distortion. However, we empirically found that each wavelength shows a different decay time, generating a bias on the computed ratio. To correct this effect, we analyze the total flux of the computed ratio images and fit an exponential function to the frames captured before the mechanical stimulus. This function is extrapolated to the following frames. Since the fluorescence does not generally return to basal conditions in the captured sequence, data points are not fitted at the end, thus ensuring a better fit. Figure 3 illustrates an example of this procedure.
Basal conditions calculation
To detect the Ca2+ wave, we have to determine the basal conditions of the experiment. This is done analyzing the ratio images (with decay correction) prior to the mechanical stimulus. They are averaged to increase signal to noise ratio. To improve the results, a sigma reject algorithm is used to prevent noise or artifacts in the averaging process (pixel values outside a 2 sigma range from the mean value are discarded).
Wave estimation
If we assume that the Ca2+ wave consists only of molecules liberated to the media, without other sources, subtraction of the basal conditions image yields the desired result. Furthermore, since no negative values are expected, except from noise, these values are truncated to work only with positive data.
Noise reduction
We used a median filter to reduce noise [9], while preserving sharp and small-scale structures. The filter kernel should be smaller or equal to 5 × 5 pixels to achieve those goals.
Statistical data output
In this last step, statistical data of the wave is saved as additional images. This includes the mean, maximum, standard deviation values and time of the maximum for each pixel. Such information will be used for further analysis in later processing steps.
Signal masking
Once we have a series of images containing the Ca2+ that is generated, we select an area of interest based on an intensity threshold. This speeds up the processes and avoids spurious data. We use images containing the statistical information from the previous step (e.g. the maximum value of the series and the average value for a given pixel) to choose the threshold as one standard deviation over the background intensity level. Although this value works for most of the cases, fine-tuning may be needed.
The result is a binary image, which is often contaminated with spurious small-scale structures associated with noise or artifacts, and with small signal voids in the excited cells. We solved these problems applying morphological closure and aperture binary operators [9]. The size of the operators (typically discs of 5 to 9 pixels in diameter) and order of execution yields different results, but running times are short, thus evaluation by visual inspection is affordable, almost in real time, testing a few configurations. This inspection is based mainly on preserving the shapes of the outer cells, without including background zones and excluding isolated sources outside the main connected group of cells.
Additionally, it is useful to restrict the maximum radius to bed of the samples to be inspected. Signal far away from the stimulus may corrupt calculations as the signal to noise ratio decreases with the distance.
Global analysis
The purpose of the global analysis is to simplify the data, assuming radial symmetry of the propagation wave. To summarize the process, a pixel is defined as the propagation center and then sample regions are automatically defined on the images. Radial distance and calcium intensity are measured in those samples, and a single output image is generated. Each pixel on the output image displays the average intensity as a function of radial distance and time. This image will be further analyzed to measure the wave velocity and to fit decay models. The process is performed in four steps:
Propagation source determination
The user needs to define by visual inspection the frame in which the mechanical stimulus was performed and calcium begins to be released. The physical location of the stimulus is found automatically from this image. To do that, the image is first filtered with gray scale morphological operations [9] (erosion followed by dilation, both 3 × 3 pixels) so that to reduce noise and the result is linearly rescaled into a 0.0 to 1.0 range (normalized range). Pixels with intensity values below 0.2 are discarded (pixel values are replaced with zeroes). The center of mass is calculated by successive approximations or iterations. Each iteration consists of: (1) narrowing the region of calculation by half on each dimension and centered on the previously calculated center of mass. (2) Updating the center of mass. This process is done until at least a rectangle of 20 × 20 pixels is defined. This approach returns the center of mass of the biggest and brightest structure in the image, which should correspond to the center of the wave of interest (Figure 4). The process ensures that secondary calcium releases or noise do not affect the calculations. If this algorithm fails, the operator may manually correct those coordinates.
Sample generation
Since we are interested in radially symmetric wave propagations, we use concentric one pixel-width ring samples. Pixel values are averaged inside each ring. Each averaged ring is mapped as a new pixel in an output image, creating a graph of radial distance and time (Figure 5).
Linear global fit
Performing a linear regression fit to the simplified data can be done to quickly extract information. We characterize the wave propagation using three parameters: time of maximum intensity, time of maximum increase of intensity, and maximum intensity, each one as a factor of the radial distance from the mechanical stimulus.
-
Time of Maximum Intensity: A linear regression of the radial distance versus time of the maximum intensity gives the mean propagation velocity of the maximum calcium concentration.
-
Time of Maximum Increase of Intensity: Instead of looking for the maximum intensity in the image, we look at the maximum positive gradient. Since the calcium release to the medium is a fast process, individual pixels present a quick and large increase of intensity. Thus, a linear regression of the radial distance versus time of the maximum gradient gives another good measurement of the wave-front velocity.
-
Maximum Intensity: Far from the mechanical stimulus, the maximum intensity that is reached can be well fitted by a linear decay function. Extrapolation of this function gives a measure of the maximum propagation radius.
Non-linear global fit
Our focus is now the time evolution of the calcium concentration for a given radial distance. For this, we used the CMinPack libraries (University of Chicago, 1999) [10], which consists of a Levenberg-Marquardt algorithm to fit non-linear models to data. To test the procedure we propose exponential, Eq. 1, and quadratic, Eq. 2, decays. Both models can be expressed in terms of three parameters, the wave maximum intensity (“A”), the time of activation (“t0”), and a decay rate (“τ”, tau).where H(t-t0) is the Heaviside function.
(1)
(2)
We analyze these fitted variables as function of the radial distance (Figure 6) to achieve a deeper understanding of the wave propagation. In particular, a linear regression of the radial distance versus measured time of activation gives another measurement of the mean wave propagation speed.
Local analysis
A non-linear fit can be performed to individual pixels, following the same procedure as in the Global Analysis. We use the Levenberg-Marquardt algorithm to find our three parameters pixel by pixel. This requires thousands of minimization operations. Masking the signal, as in section 2.2, greatly improves the efficiency of this procedure, since we analyze only those pixels that are part of the wave. The result gives a localized map of the variables that characterize the calcium concentration behavior in time.
One filter wave propagation measurement
To calculate the calcium concentration, consecutive filtered images are acquired at two different wavelengths. If we are interested only in measuring the wave-front velocity, calculations can be done using only the 380 nm set. By using only one filter we are saving the time of the second exposure and that required to mechanically change the filters, allowing a higher sampling rate (from 0.8 to 2.8 frames per second in our setup). This yields to a more reliable estimation of velocity, since it is derived from the slope of the maximum gradient linear regression, in the global analysis.
The characterization procedures remain the same but without computing ratios.
The propagation velocity can be computed from both, the time at the maximum pixel intensity or the maximum increase in pixel values. As will be shown, the mean velocity of the calcium wave-front estimated from the maximum gradient is highly correlated to that found using two filters. Maximum intensity in the 380 nm images is not well correlated to wave maximum intensities found with two filters, thus estimations of the mean propagation velocity cannot be done from the time of maximum intensity.
Implementation details
All stages were implemented in C++ language, as external modules for the image processing platform PixInsight (Pleiades Astrophoto S. L., Valencia, Spain). Raw images were stored in an 8-bit TIFF format. Loaded images were transformed to 32-bit floating-point samples, and outputs were stored as 32-bit FITS files to maintain numerical accuracy.
Experimental setup
Primary cultures of oviductal ciliated cells
Epithelial cell cultures from rat oviduct were obtained as in [11]. Female adult rats Sprague Dawley (250-300 g) were used. Animals were kept under artificial illumination and controlled temperature (23-25°C) from 6:00 to 21:00 hours. Water and pelleted food were supplied ad libitum. Day 1 of the rat estrous cycle was assessed by the presence of a vaginal plug. On day 1, between 11:00 and 12:00 hours, animals were sacrificed by cervical dislocation after ketamine/Xilacine anesthesia i.p.
The genital tract was removed and oviducts were placed into Hanks’ solution (Sigma-Aldrich) at pH 7.4. The distal portion of the oviduct (i.e. the fimbria) was mechanically removed with fine forceps. The ciliated epithelium was extracted and transferred to a fresh culture medium. Epithelial explants were placed on 0.1% gelatin-treated coverslips of sealed chambers. The Rose Chambers were perfused with NHS medium containing 10% heat-inactivated horse serum. NHS had the following composition in mM: NaCl 137; KCl 5.09; Na2HPO4 1.14; KH2PO4 0.18; MgCl2 0.923; CaCl2 0.91; NaHCO3 4.07; glucose 21.5 and glutamine 0.2; supplemented with 1.0% vitamins, 1.0% essential amino acids, 1.0% non-essential amino acids, 1.0% pyruvate and antibiotics (0.2 mg/ml neomycin and 0.12 mg/ml penicillin), pH 7.2-7.4 at 37°C. Primary cultures were incubated at 37°C for 5-7 days. When this cell monolayer showed spontaneous ciliary activity, the chambers were opened and the cultures were washed three times with fresh Hanks’ solution before calcium measurements were performed.
These protocols were approved by the Animal Care and Ethics Committee of Pontificia Universidad Catolica de Chile.
Measurement of calcium waves
From manual estimations with a few samples, we observed some differences on the calcium wave propagation when the ciliated epithelium samples were exposed to different hormones [12]. Thus, we decided to perform three kinds of experiments: cultures as described in the previous section, and cultures exposed either to leptin or adiponectin.
Intracellular calcium levels were determined using a spectrofluorometric technique in primary cultures of oviductal epithelial cells as described in [4], with Fura 2-AM (acetoximetil ester, Invitrogen Molecular Probes, Inc. OR 97402-9132, USA). Intracellular calcium wave propagation velocity was measured in ciliated cell cultures previously incubated with different concentrations of leptin (Recombinant Rat leptin, R&D Systems, Inc., Minneapolis, MN 55413 USA), and adiponectin (Acrp30, R&D System, Inc. Minneapolis, MN 55413 USA) for a period of 48 hours.
Calcium waves were generated with a mechanical stimulus, which consisted in a borosilicate glass micropipette (Harvard Apparatus Lts, UK) for patch-clamp, mounted on a micromanipulator. The micropipette descended gradually to touch the surface of a cell, generating an increase of calcium, which spreads to neighboring cells. Video microscopy images were obtained in a fluorescence phase microscope (Olympus Bx 51W1L, Olympus Optical, Tokyo, Japan) with water immersion lens coupled to an imaging detection system (QImaging Retiga 13001 CCD, cooled monochromatic digital camera 12-bit, QImaging Burnaby, Canada). Images were captured with Metafluor software (Universal Imaging Corp., PA, USA). To provide a sufficient number of frames for proper basal conditions estimation, at least 10 seconds were captured before the mechanical stimulus at each experiment. Experiments where stopped when the operator could not see any significant intensity change in the ratio image displayed by the Metafluor software.