Introduction

Cystic fibrosis (CF) is a genetic disease that affects over 70,000 patients worldwide. The primary cause of death in CF is respiratory failure. Mutations cause dysfunction of an anion channel at the surface of epithelial cells and result in osmotic imbalances. This causes mucus layers above the cells to dehydrate, inhibiting a process known as mucociliary clearance (MCC) in which hair-like structures at the surface of cells called cilia sweep away mucus. In CF, mucus dehydration causes MCC failure, which contributes to increased infection and inflammation. Nuclear imaging techniques provide a means of measuring MCC, however, they do not provide a clear picture of the underlying processes driving mucus clearance. Previous work has mathematically described the dynamics of MCC at the population level in different regions of the lung [1]. On an individual basis, however, there is visible heterogeneity within these regions and the model structure does not capture MCC dynamics for all subjects. Other work has studied MCC on a pixel-by-pixel basis in individuals, however the sampling frequency likely was too long to capture dynamic differences in the large airways [2]. The present work aims to bridge this gap by developing a mathematical model of MCC in 2D, based on average airway physiology. In doing so, we can gain valuable insight into clearance-rate limiting regions of the airways and quantify differences in MCC dynamics on a per-subject basis.

Methods

Nuclear Imaging

Study subjects inhaled two radiolabeled probes, a non-absorbable particle, Technetium-99m labeled sulfur colloid (Tc-SC) and a small molecule, Inidium-111 labeled diethylene triamine pentaacetic acid (In-DTPA) prior to 80 minutes of sequential gamma camera imaging. We focus here on Tc-SC, which provides a surrogate to measure mucus clearance and thus MCC (In-DTPA provides an additional assessment of paracellular absorption). Subjects also inhaled isotonic saline from t=10-20 min during the measurement. Healthy (n=12) and CF (n=22) subjects were included.

Image Processing

The images were processed using ImageJ (1.52v). An anatomically-based whole lung region of interest (ROI) taken from [3] was stretched and rotated to fit the right lung of a posterior transmission scan for each subject. This ROI was then translated vertically and horizontally on the posterior view of the nuclear images to align the ROI with the observed activity. The images were clipped to a bounding box around this ROI for all time points. These cropped images were divided into 8 x 16 (w x l) grids and the total Tc-SC activity in the whole lung was normalized to a value of 100 for each subject after correcting for background activity and radioactive decay of the probe.

Model Construction

To capture the branched nature of the airway tree across the lungs, we defined the large airway region (LAR) of the 8 x 16 grid. This was done by dividing the right lung of an anatomically-averaged high resolution computed tomography (HRCT) scan from [4] into an 8 x 16 grid and defining the LAR as any grid that overlapped the large airways, as shown in Figure 1. From this we constrained where mucus, and thus Tc-SC, could flow in 2D. We defined these constraints using the elevation mapping shown in Figure 1, such that mucus within the LAR flowed within the LAR towards the trachea (left edge of LAR) and mucus outside the LAR flowed toward the LAR; in terms of the elevation mapping, this meant that mucus could not flow “uphill.” Grids outside the whole lung region were fixed at an activity of 0 and mucus could not flow into or out of these grids.

[Figure 1 attached below]

The dynamics of Tc-SC clearance in each grid block were defined using mass-action kinetics. Assuming that MCC in a specific block occurred at a constant rate and flowed in all directions permitted by the flow constraints at equal rates, the concentration of activity in each block was described by the Equation 1.

[Equation 1 attached below]

Here C_k is the concentration of activity in grid k, k_k is a fitted rate coefficient for MCC leaving grid k, J_in is the set of all grids that can flow into grid i, and N is the number of grids that mucus from grid i can flow into. Rate coefficients were fit to Tc-SC retention in each block over 80 minutes on a per-subject basis by minimizing the sum of squared errors.

Model Reduction

To improve interpretability and decrease parameter uncertainty, the number of model parameters was reduced in two steps: (i) spectral clustering was used to identify contiguous 2D clusters with similar rate coefficients; and (ii) linear relationships were defined between rate coefficients of clusters that were highly correlated. In the first step, rate coefficients in the same cluster were set to be equal (i.e., a single parameter for each cluster) and in the second step, linear relationships were defined (i.e., only one rate coefficient in a highly correlated pair was independent). In each case, the dynamics describing the system were not modified, but the rate coefficients describing MCC leaving each grid were constrained to reduce the number of free parameters.

Results

The model was able to capture MCC dynamics in both CF and healthy control subjects with substantially higher accuracy than in previous work, particularly in the peripheral regions. The model from this work had a mean squared error (MSE) of 12.4±7.5 and 30.7±19.1 in the large airway and peripheral regions, respectively, versus 25.3±23.8 and 176±247 using the model from [1]. Furthermore, the number of fitted parameters could be reduce while retaining a high level of accuracy. Using spectral clustering to decrease the number of fitted parameters from 120 to 26 only increased the MSE in the large airway and peripheral regions to 13.9±7.7 and 33.6±19.8, respectively. Further reduction by fixing linear relationships between strongly correlated parameters showed that these dynamics could be captured using just 12 parameters, maintaining high accuracy (MSE of 20.3±30.8 and 50.1±32.9 in the large airway and peripheral regions, respectively).

Summary

This work describes a 2D model of MCC across the lung that improves upon existing work in its ability to capture variations between individuals using a small number of free parameters. Data is available for the same subjects that includes inhalation of a mucus hydrating agent (7% hypertonic saline), which will be used to characterize CF patient MCC improvements with hyperosmotic treatment. A further extension will explore the dynamics of paracellular junction permeability to water absorption in healthy and CF subjects, which is facilitated by the small molecule probe data from the original imaging study.

References

[1] M. R. Markovetz, T. E. Corcoran, L. W. Locke, M. M. Myerburg, J. M. Pilewski, and R. S. Parker, “A physiologically-motivated compartment-based model of the effect of inhaled hypertonic saline on mucociliary clearance and liquid transport in cystic fibrosis,” PLoS One, vol. 9, no. 11, pp. 1–13, 2014.

[2] W. D. Bennett, M. Xie, K. Zeman, H. Hurd, and S. Donaldson, “Heterogeneity of particle deposition by pixel analysis of 2D gamma scintigraphy images,” J. Aerosol Med. Pulm. Drug Deliv., vol. 28, no. 3, pp. 211–218, 2015.

[3] L. Alcoforado et al., “Anatomically Based Analysis of Radioaerosol Distribution in Pulmonary Scintigraphy: A Feasibility Study in Asthmatics,” J. Aerosol Med. Pulm. Drug Deliv., vol. 31, no. 5, pp. 298–310, Oct. 2018.

[4] E. E. Greenblatt, T. Winkler, R. S. Harris, V. J. Kelly, M. Kone, and J. Venegas, “Analysis of three-dimensional aerosol deposition in pharmacologically relevant terms: Beyond black or white rois,” J. Aerosol Med. Pulm. Drug Deliv., vol. 28, no. 2, pp. 116–129, Apr. 2015.