Assessing motile cilia coverage and beat frequency in mammalian in vitro cell culture tissues

Cilia density, distribution and beating frequency are important properties of airway epithelial tissues. These parameters are critical in diagnosing primary ciliary dyskinesia and examining in vitro models, including those derived from induced pluripotent stem cells. Video microscopy can be used to characterize these parameters, but most tools available at the moment are limited in the type of information they can provide, usually only describing the ciliary beat frequency of very small areas, while requiring human intervention and training for their use. We propose a novel and open-source method to fully characterize cilia beating frequency and motile cilia coverage in an automated fashion without user intervention. We demonstrate the ability to differentiate between different coverage densities, identifying even small patches of cilia in a larger field of view, and to fully characterize the cilia beating frequency of all moving areas. We also show that the method can be used to combine multiple fields of view to better describe a sample without relying on small pre-selected regions of interest. This is released with a simple graphical user interface for file handling, enabling a full analysis of individual fields of view in a few minutes on a typical personal computer.

closely resembles in vivo conditions and drives differentiation towards a mucociliary phenotype [1][2][3]. The tissue can be started from primary cells or cell lines. Moreover, in the last decade, multiple protocols have demonstrated the ability to derive hAECs from induced pluripotent stem cell (iPSCs) [4,5].
In vitro models of human airway epithelium are an effective tool to study mucociliary clearance (MCC) and disease pathogenesis in the respiratory tract [6]. However, depending on the cell origin, culture medium and growth factors employed, these cultured tissues may differ for the amount and type of cells (i.e. ciliated, basal, secretory cells). Therefore, number, distribution and ciliary length, and physical properties of the airway surface liquid (ASL) are often highly diverse among the samples, leading to a variability in cilia beating frequency (CBF), coordination and alignment [7,8]. Moreover, the ability of motile cilia to generate an effective MCC can be altered in the presence of genetic diseases such as cystic fibrosis (CF), where the ASL hydration level decreases drastically, or primary ciliary dyskinesia (PCD), where cilia have defective function [9]. Hence, the need to measure coverage in motile cilia and CBF in in vitro mucociliated tissues in order to quantify physical differences between different cells and different culture protocols and be able to correlate them to the MCC efficiency.
While improvements have been made in recent years, there are still a limited set of tools for the analysis of microscopy videos of ciliated cells, with the standard being the commercial software Sisson-Ammons Video Analysis (Ammons Engineering, http://www.ammonsengineering.com/SAVA/sava.html) [10]. Even though new alternatives have appeared in recent times [11], this method typically only measures the CBF and mostly in user selected areas. This usually requires human intervention in either the areas selection or the frequency range filtering, a potential source of bias in the analysis. In addition, and unlike the method presented in this manuscript, SAVA software does not evaluate how much of the visualized area is covered with beating cilia, and thus only allows for CBF comparison of small preselected areas. Alternative techniques have been developed over the years, spanning from methods that make use of cheaper cell phone cameras for image acquisition [12], to others attempting to eliminate or automate the region of interest (ROI) selection [13], or trying to determine spatial distribution of CBF [14]. However, these techniques are still limited to CBF measurements only, and mostly limited to small areas of the epithelium or only able to run on video of cilia imaged from the side. A different method must thus be developed to better quantify the differences in beating cilia coverage and frequency, particularly in heterogeneous samples that can require the acquisition of a high number of fields of view (FOV) to obtain a representative characterization, and for samples that can present the full range from no visible moving cilia to full coverage in motile cilia. Identifying regions with cilia beating is also a prerequisite for other analysis, for example multiDDM [15] aiming to measure the spatial coherence of ciliary dynamics.
The approach presented here can, with no input from the user, accurately detect areas containing beating cilia, calculate the CBF distribution across the entire FOV and, given a sufficient number of FOV for a particular sample, provide a description of the average percentage of tissue covered in motile cilia, and aggregate to a CBF distribution over an entire sample of FOV from matching conditions, as illustrated in figure 1.

Ciliated epithelia
The algorithm developed here has been tested on videos collected from in vitro airway epithelium derived from both human primary cells and iPSC. Commercially available human Bronchial Epithelia reconstituted in vitro (MucilAir) were purchased from Epithelix Sàrl and maintained following the protocol provided by the company. Other samples have been derived from human iPSC by Laetitia Pinte and Marta Vila Gonzalez at the Wellcome-MRC Cambridge Stem Cell Institute. Briefly, human iPSCs were cultured as previously described in [16]. For differentiation towards airway epithelia, protocols detailed in both [17,18] were used with minor modifications. In both cases, epithelial tissues are reconstituted on a microporous filter in an plastic insert of 6.5 mm diameter and imaged to assess cilia properties at the ALI stage.

Microscopy set-up
For these measurements, bright field high-speed videos are typically recorded at least at approximately 80 frames per second (fps) and saved in an uncompressed video format. The algorithm should work with any file type that can be opened with MATLAB'S own VideoReader but we would recommend recording without compression to avoid any possible artefacts.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 230185 Our set-up consists of a Nikon Ti-E inverted microscope equipped with a 60 ×-water immersion Nikon objective (NA 1.2), using a CMOS camera (model no. GS3-U3-23S6M-C; Point Grey Research/ FLIR Integrated Imaging Solutions, resulting in 1 pixel ¼ 0:098 m. Each FOV ≈ 187 μm by 117 μm). Samples of ciliated cells are imaged in a custom-made chamber which regulates temperature, CO 2 and humidity continuously, maintaining values of 37°C, 5% and 100%, respectively. The insert is placed in a glass-bottom Petri dish and the semipermeable membrane is kept on a thin layer of culture medium during the experiments. The objective used was a CFI Plan Apochromat VC 60XC WI with a working distance of 0.29 (0.31-0.28) mm. The thickness of the bottom of the Petri dish is the same a cover glass slide (0.17 mm) and the membrane is placed directly on the glass. With this set-up, it was possible to clearly focus over the entire height of the epithelium and cilia layers.

Software development
One of the basic parameters used to evaluate ciliary function and MCC efficiency in airway epithelial tissues is the CBF [9,19]. The majority of algorithms employed for CBF measurement are based on fast Fourier transform (FFT) [10,20]. Briefly, CBF can be extracted from the frequency of light intensity fluctuation by Fourier transforming the pixel intensity over time, therefore measuring the ciliary cycles. A common limitation of these methods is that to determine the relevant frequency of the pixels across the FOV, it is required to distinguish between areas where cilia are beating and areas without cilia movement. Only the first group should contribute to the final CBF measurement. This might seem like a simple task, but it is in fact hard to develop a robust and automated method working reliably in the general range of possible conditions (very low to high coverage). It is therefore common for these techniques to require human intervention in the selection of the areas to analyse, creating a potential source of bias and limiting the number of FOV that can be analysed in a given time.
For this purpose, an algorithm was developed relying on the differences in frequency and intensity in the power spectrum density (PSD) of areas with and without beating cilia, over different length scales, to first determine the areas where movement is present, followed by a Fourier power spectrum analysis to determine the CBF of the ciliated areas. The first approach identifies areas where beating cilia are likely to be located by observing how the power spectrum changes between regions with and without movement. Areas with no movement are typically characterized by a Fourier power spectrum with higher amplitudes in the lowest frequencies, mostly representing the background noise inherent to the camera and acquisition settings. On the other hand, areas with high cilia coverage will present higher and sharper peaks for relatively bigger frequencies.
The development and video analysis were performed using MATLAB R2022b (MathWorks Inc, Natick, MA, USA), running on a Linux workstation (AMD Ryzen 3900X, 32 Gb RAM, Nvidia RTX 3090 24 Gb).   Figure 1. The method presented in this paper allows us to determine the areas of an ALI culture where cilia are beating, in conditions of high (a), medium (b) and low (c) coverages, and even a single cilia patch (d ). The frequencies of the pixels in the areas where movement is detected are organized into a histogram to visualize the CBF distribution. The characteristic frequency for the FOV is calculated as the median of the distribution (red line) and the scatter of the frequencies (in orange) has been quantified by measuring the corresponding 25th and 75th percentiles. Each sample measurements is described above the histogram as median (25th percentile-75th percentile). Each FOV is ≈187 μm by 117 μm.

Initial image normalization
As the illumination quality may vary across the different videos analysed, and in order to achieve consistency in the grey-scale range among them, each video is first normalized. The normalized intensity at height i, width j and time t (I n (i, j, t)), is obtained by subtracting the minimum value measured in the entire video (min (I)) to the original intensity signal I(i, j, t), and dividing the result by the intensity range This operation sets the value of all pixels in the video in a range between 0 and 1. Background intensity is then decreased by calculating the average intensity for each pixel over time, and subtracting it from each pixel obtaining I n . This removes most of the background and normalizes the framestack to values between −1 and 1.

Pixel Fourier transform analysis
The frame-stack auto correlation is then calculated, and this processed frame-stack is subsequently used in the FFT analysis. Pixel PSD P(i, j, f ) is estimated for each position (i, j ) through FFT over time, applying a Hann window to the signal I n . This consists in multiplying the signal by a sinusoidal function that touches zero at both ends eliminating all discontinuity at the edge of the signal. The resulting spectrum is filtered through a band pass filter (2-30 Hz).

Box average function
Areas with beating cilia within the FOV are determined by observing how each pixel power spectrum changes within square 'boxes' of increasing size, centred in each pixel. Fifteen different box sizes between 5 and 33 pixel are first determined Then, a 'box average function' is applied to the original spectrum for each box size. That is, the spectrum for a given pixel is calculated as the average of the spectra of the B size surrounding pixels, for the different box sizes. A visualization of the key steps in the process is described in figure 2.
Formally, let the input be a three-dimensional array, P(i, j, f ), with dimensions H × W × F, where each layer in f corresponds to the amplitude of the PSD for all pixels in a FOV with i height and j width for a given frequency f. Let the output be a smoothed version of this stack, A s . The box average function is applied by dividing each slice in the stack into overlapping blocks, or 'boxes' of size B size × B size . The output for each pixel A s (i, j, f ) is then calculated by taking the average amplitude of all the pixels within the box centred in (i, j ), for each layer f in the stack.

Movement map determination
In the resulting three-dimensional matrix, each element A s (i, j, f ) represents the average power spectrum amplitude at the frequency f evaluated in a box of size s centred in the pixel (i, j ). To determine the movement map for the FOV, we consider how the box-averaged power spectra contained in A s changes when increasing box size by comparing the slope of their linear fit. Pixels far from any area with active movement present a PSD similar to that of their close neighbours. Thus, by averaging the spectra on different box-sizes, the curves should be as similar as their linear fit slopes. Pixels within an area with beating cilia display variations in the PSD when changing box size, and thus variations in their corresponding linear fit.
First, the three-dimensional array A s is simplified by fitting each spectrum with a linear least square fit and storing the slope of the fit (see figure 3 for a visual reference of the overall sequence of steps). Let the input be the previously computed three-dimensional matrix, A s , and let the output be a matrix of the slopes of the best linear fit for the amplitudes along the f axis of A s , indicated as the matrix a s (i, j ). Let F be the length of f. In this case, for each s we apply a linear least square fit to each element of the matrix A s along its third dimension, using f as the independent variable. a s ði, jÞ ¼  Note that with this linear fit (figure 3a), we were not trying to achieve a good interpolation of the data; the purpose of the fit was to identify in a numerically rapid and robust fashion if there is an outstanding low-frequency peak that rises above background or noise. This information is obtained in a particular pixel (i, j ) for each box sizes, resulting in 15 different two-dimensional matrices a s , s = 1, …, 15. The matrices are stacked together is a three-dimensional matrix M(i, j, s) of size (H × W × 15), i.e. an array of 'box slopes' (see figure 3a for a visual representation of this step, for the case of a single pixel).
A second linear least-squares fit is then applied over the box size axis s, again storing the value of the calculated slopes, resulting in a final matrix describing the profile changes over the different box sizes for each pixel spectrum in the entire FOV, m(i, j ). Pixels in a background region, that is, away from any movement area, should have relatively stable spectra profiles over the different box sizes, thus low m values. Pixels within a high movement area display slightly higher variations with different box sizes due to small local variations in cilia beating frequency. On the other hand, pixels on the edge between a high movement area and a background area should present a decreasing or increasing slope profile, depending on whether they are within the movement area or just outside of it (figure 3b,c). The resulting matrix m is again box-averaged with a five pixel box size to remove noise, and a binary threshold is applied, so that any absolute value under 0.6 × 10 −7 is considered background. This value was determined experimentally as the most flexible regarding FOV with high, medium and low coverage, as well as unfocused FOV, in order to present the least amount of artefacts (figure 3d). The obtained movement map is later used to calculate the cilia coverage as the percentage of area containing movement within the FOV.

Cilia beating frequency
To determine the distribution of frequencies within the FOV, the original pixel PSD map P is firstly box averaged with a three pixel box size to remove noise. The characteristic frequency for each pixel is then determined by calculating the local maxima of the PSD for each box, and selecting the one presenting the highest amplitude (figure 4a), generating a frequency map of the FOV (figure 4b). The resulting matrix is again filtered to remove any artefact resulting from the box average operations before, by removing an edge of pixels on the outside of the FOV correspondent to the size of the largest box used, and is then combined with the previously determined movement map creating a frequency map for the areas containing movement (figure 4c). The characteristic frequency for the FOV is calculated as the median frequency measured, and the corresponding 25th and 75th percentiles are also evaluated to quantify the scatter of the frequencies in the FOV (figure 4d).
In this way, it is possible to determine the beating cilia coverage for FOV with different coverage, from low to high, with minimal error, even when only small patches of cilia are detected.

Detecting focussing issues and debris
As already discussed, in vitro cultured hAEC tissues can significantly differ from each other for coverage, CBF, and ASL properties depending on the cell type and the culture technique. Moreover, even following a well-defined imaging protocol, the collected videos may vary for the focusing choices or the presence of blurred elements and debris. For this reason, we implemented the algorithm with a debris detection and a bad focus warning. Both these sections are based on each pixel intensity power spectrum and standard deviation data.

Bad focus detection
The cells lining the tract of the human airways going from the trachea to respiratory bronchioles form a pseudostrafied epithelium. Here, even though all cells make contact with the basal membrane creating a single cell layer, the nuclei are disposed on different levels, causing the illusion of cellular stratification. In some samples, this effect is more pronounced, making it challenging to have all the ciliated cells entirely in focus in the same focal plane. Moreover, the presence of mucus may hinder the correct focusing on cilia royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 230185 when this layer is extremely viscous, leading to a globally blurred image (figure 5). It is possible to attempt to identify a blurred image by its STD distribution. For each video, the STD is computed as I n ði, j, tÞ À hI n ði, j, tÞi t Á 2 s , obtaining a H × W matrix in which each element assumes the value of the pixel standard deviation over time. The computed STD values are normalized to a range between 0 and 1, and the matrix is box averaged with a box size of 10 pixel. The STD distribution in the FOV is then fitted with a log-logistic  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 230185 regression, and the resulting mean (μ) and spread (σ) are calculated. With the STD histogram being quite sharp when the video is blurred, it is possible to implement a bad focus detection section based on the spread of the distribution. If σ < 0.05, a warning message is triggered, and a flag is added into the output file, so that the output provided can be further verified for accuracy.

Debris detection
It is not uncommon for videos of biological samples to present unwanted debris crossing the FOV. These debris can be detached cells, mucus vesicles or dust particles. While the analysis method described here appears to be robust enough that the presence of small debris should not significantly affect the final result, it still takes note of cases where debris are detectable. The PSD of the pixels crossed by debris when transported across the FOV can present more noise, resulting in less sharp peaks. This can be quantified by measuring the prominence of the highest peak in the PSD (i.e. how much the peak stands out from the surrounding baseline of the signal due to its height and its position compared to other peaks). If a local maxima is part of a range of peaks (i.e. noisy PSD representing background), it is more likely to have a smaller prominence compared to an isolated, sharp one, which commonly describes cilia beating. This information can be collected in a ratio map R(i, j ), where every entry represents the ratio between the prominence and the amplitude of the highest local maxima (both a direct output from the CBF calculation in the previous section). Once computed, the matrix R is box averaged with a 10 pixel box size. The entries in R with values close to 1 are more likely to represent areas with beating cilia, due to the presence of a sharp peak in the PSD. On the other hand, if the value of the ratio is small, the pixel presents a noisy spectrum and it is likely to represent an area with debris moving or background. In addition, the presence of debris can result in an increase in the STD for the pixels in that specific region, due to the variation in the intensity signal compared to the average. The STD map σ video (i, j ) is calculated and box averaged as described above.
Since the pixels potentially affected by the debris movement are characterized by a high STD signal (moving high contrast feature) and a small ratio value (noisy PSD), a debris mask can be finally calculated as Entries true in the mask would represent pixels possibly crossed by a debris while the video has been captured. If the total number of pixels tagged in the debris mask is higher than 2% of the total pixel count, the FOV is tagged as potentially containing debris ( figure 6).

Sample data aggregation
To fully characterize a sample while avoiding any potential bias, it is important to collect several FOV across the entire tissue and aggregate their information. This algorithm includes an aggregation function that can calculate the frequency distribution and average CBF for the entire sample. This is done by pooling together the pixel information relative to each video as one larger FOV describing the whole sample. royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 230185

Graphical user interface
To make it easier for the user to run the analysis on multiple samples and FOV, it is useful to work through a graphical user interface. This allows the user to input the folder where the data is stored, the output folder where the results should be written to, as well as filename based sample identification and renaming (figure 7). This input data are saved as a configuration file in the output folder so that in the case of any unexpected shutdown the analysis can continue without the need to re-input all the information.

Detecting differences in cilia beating frequency and coverage
It is important to be able to clearly distinguish small differences in both beating cilia coverage and their corresponding frequency. Figure 1 shows that this method can easily identify the moving areas in a FOV, from fully covered tissue to very small patches of beating cilia. The algorithm also provide a way to show differences in the CBF distribution for the different videos, in this case displaying a relative decrease in CBF following the decrease in ciliated areas.

Sample characterization from multiple FOV
In order to measure CBF and coverage across the entire sample, our typical protocol is to collect 20 FOV for each insert, in predefined un-biased positions as described in figure 8a. These are selected to include information from both the middle and on the edge areas of the sample, while avoiding any potential royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 230185 selection bias. To estimate a final characterization of a sample, the frequency and movement maps obtained for each FOV are pooled together (figure 8b). As described above, the motile cilia percentage is calculated as the accumulated percentage of pixels moving across all the FOV, that is, the total number of pixels classified as moving accumulated over all FOV, divided by the total number of pixels in all FOV for a given sample. The CBF is calculated as the median frequency of all 20 FOV, with the corresponding 25th and 75th percentiles (figure 8c).

Run time and operation
Run time may vary between different systems depending on their hardware configuration. The current version of the algorithm makes use of MATLAB'S GPU acceleration when available. We tested the average run time for the system described in the previous sections, both with GPU acceleration and as CPU only. A FOV movie was first loaded into memory, followed by five sequential repeated analyses. While the use of the GPU can drastically decrease the analysis time for each FOV (from about 300 s to about 75 s in this configuration), it is still possible to fully characterize a sample from 20 FOV in about 100 min while only using the CPU (compared with 25 min while using the GPU), highlighting the quick turnaround time for a full sample analysis.

Current limitations
While the method described can be a powerful tool in the analysis of ciliated samples, it does present some limitations. First, it can only detect motile beating cilia, and while this is an important parameter to quantify, other assays may be needed to fully characterize a sample regarding total cilia coverage when these are not motile. In addition, as mentioned before, the focusing can have significant effects on the performance of this method. It is up to the user to make sure that a sample is correctly focused, however some challenges may appear when imaging samples presenting more three-dimensional features with large shifts in focusing planes.

Conclusion
In this manuscript, we presented a new powerful method for the characterization of biophysical properties of airway ciliated epithelium, able to reliably measure CBF and coverage in motile beating cilia in an automatic fashion. The algorithm we developed is able to independently isolate areas with beating cilia from the background, estimate the distribution of characteristic frequencies in the video, and combine information from multiple collected FOV to characterize the entire sample. The ability to automatically detect motion areas removes a source of human bias in the ROI selection, at the same time decreasing the knowledge barrier required to make use of the software, and increasing the throughput of FOV that can be analysed in a short amount of time. In addition, and unlike most alternatives so far, the ability to quantify and compare the areas of the tissue that presents moving  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 230185 cilia allows for a more thorough characterization of a sample, not only providing information about the frequency at which cilia beat in each specific region, but also measuring how they are distributed across the sample and how this affects their movement. This method also does not require sideways cilia imaging, thus removing further sample processing steps, and is capable of isolating the movement of even small patches of cilia from the background.