Modelling Airway Geometry as Stock Market Data using Bayesian Changepoint Detection

TL;DR

This paper introduces a Bayesian changepoint detection method to accurately identify airway dilatation progression in lung images, demonstrating robustness against noise and variability, with potential for disease monitoring.

Contribution

A novel probabilistic model and inference approach for detecting airway dilatation changes from longitudinal imaging data, improving accuracy and robustness over existing methods.

Findings

Detected airway dilatation with 2.5mm accuracy on simulated data

Reasonable agreement with radiologists on real IPF data

Quantified airway volume change for disease progression assessment

Abstract

Numerous lung diseases, such as idiopathic pulmonary fibrosis (IPF), exhibit dilation of the airways. Accurate measurement of dilatation enables assessment of the progression of disease. Unfortunately the combination of image noise and airway bifurcations causes high variability in the profiles of cross-sectional areas, rendering the identification of affected regions very difficult. Here we introduce a noise-robust method for automatically detecting the location of progressive airway dilatation given two profiles of the same airway acquired at different time points. We propose a probabilistic model of abrupt relative variations between profiles and perform inference via Reversible Jump Markov Chain Monte Carlo sampling. We demonstrate the efficacy of the proposed method on two datasets; (i) images of healthy airways with simulated dilatation; (ii) pairs of real images of IPF-affected airways acquired at 1 year intervals. Our model is able to detect the starting location of airway dilatation with an accuracy of 2.5mm on simulated data. The experiments on the IPF dataset display reasonable agreement with radiologists. We can compute a relative change in airway volume that may be useful for quantifying IPF disease progression. The code is available at https://github.com/quan14/Modelling_Airway_Geometry_as_Stock_Market_Data

Tables

Table 1. Table 1: Table of the image properties of voxel size, number of airways used and the time between scans. The airways were selected by a trained radiologist R1. All patients in Experiment 1 do not have IPF. Abbreviation: BL - Baseline Scan, FU - follow-up Scan, M - Months, D - Days.

Experiment	Patients	BL Voxel Size	FU Voxel Size	Airways	Time between scans
1	1	0.67, 0.67, 1.00	0.56, 0.56, 1.00	6	9M 6D
1	2	0.63, 0.63, 1.00	0.78, 0.78, 1.00	7	35M 6D
1	3	0.72, 0.72, 1.00	0.72, 0.72, 1.00	1	5M 22D
2	1	0.72, 0.72, 1.00	0.64, 0.64, 1.00	3	12M 5D
2	2	0.67, 0.67, 1.00	0.87, 0.87, 1.00	1	10M 24D

Table 2. Table 2: The percentage volume change (PVC % percent \% ) for each region of the airway.

Airway	PVC of $V_{c\to d}$	PVC of $V_{t\to d}$	PVC of $V_{c\to t}$
1	2.6 %	32.9%	1.1%
2	3.2 %	129.7%	2.2%
3	2.6 %	47.4 %	-0.3%
4	7.4 %	48.4 %	7.1%

Equations

$$a={\mathrm{argmin}}_{a\in [-5,5]}{\left|\left|\log \left(\frac{g_B(x)}{g_F(x-a)}\right)\right|\right|}_2.$$

$$\mathbf{y}=\{\begin{array}{l}(y_1,\dots ,y_{\tau })\sim F({\theta }_1)\\ (y_{\tau +1},\dots ,y_n)\sim F({\theta }_2)\end{array},$$

$$p(\mathbf{y}|\tau ,\theta ,M)=\prod _{l=1}^{k+1}F(y_{{\tau }_{l-1}:{\tau }_l}|{\theta }_l),$$

$$\mu$$

$${\sigma }^2$$

$$\nu$$

$$M$$

$$\begin{split}\alpha=\textup{min}\Bigg{\{}1,\frac{p(\textbf{y}|{\tau_{new}},{\theta_{new}},{M_{new}})}{p(\textbf{y}|\tau,\theta,M)}\frac{p({\tau_{new}})p({\theta_{new}})p({M_{new}})}{p(\tau)p(\theta)p(M)}\\ \frac{q(\tau|{\tau_{new}})}{q({\tau_{new}}|\tau)}\frac{q(\theta|{\theta_{new}})}{q({\theta_{new}}|\theta)}\frac{q(M|{M_{new}})}{q({M_{new}}|M)}\,|J_{M\rightarrow M+1}|\Bigg{\}}\end{split}$$

$$\begin{split}\alpha=\textup{min}\Bigg{\{}1,\frac{p(\textbf{y}|\tau{new},{\theta_{new}},{M_{new}})}{p(\textbf{y}|\tau,\theta,M)}\frac{p(\tau_{new})p({\theta_{new}})p({M_{new}})}{p(\tau)p(\theta)p(M)}\\ \frac{q(\tau|\tau_{new})}{q(\tau_{new}|\tau)}\frac{q(\theta|{\theta_{new}})}{q({\theta_{new}}|\theta)}\frac{q(M|{M_{new}})}{q({M_{new}}|M)}\,|J_{M+1\rightarrow M}|\Bigg{\}}\end{split}$$

$$\mathbf{l}=\frac{M}{1+e^{-k(x-\alpha )}}.$$

$$J(K)=\sum _{r=0}^{K-1}\sum _{i=k_r}^{k_{r+1}-1}\Delta (y_i,y_{k_r:k_{r+1}-1})+\beta K,$$

$$\Delta (y_i,y_{k_r:k_{r-1}-1})=y_i-\mathrm{mean}(y_{k_r:k_{r+1}-1})$$

Full text

1 Introduction

In fibrosing lung disease, contraction of the lung interstitium pulls on airway walls, and causes a dilatation. Airway dilatation is calculated using crude visual scores have been shown to be powerful predictors of outcome in idiopathic pulmonary fibrosis (IPF) [8]. Computed Tomography (CT) coupled with image analysis algorithms to measure the airway lumen can provide an excellent tool to quantitatively assess these dilatation. CT gives a very high image contrast between air within the airway lumen and the surrounding airway wall. However, airway lumen measuring algorithms are subjected to sources of noise. Variations in voxel sizes [1] and CT reconstruction algorithms [6] have been shown to affect the precision of the cross-sectional area measurements of the airways. In addition, there can be variations in airway cross-sectional area due to normal biological processes such as points of airway bifurcation where the airway transiently dilates [12, 3]. These sources of noise makes post-processing of measurements difficult to locate dilatation on an airway track. Identifying the point of dilatation and subsequently calculating the volume of abnormal airways in the lung potentially provides a sensitive means for quantifying subtle worsening of disease in IPF across longitudinal scans.

In this work, we introduce a variant of Bayesian Changepoint Detection (BCPD) models to automatically identify the location of airway dilatation in the presence of strong measurement noise, given two measurements of the same airway acquired at different time points. BCPD model, typically used in DNA sequencing [13] and stock market data analysis [10], aims to capture abrupt variations in the underlying distributions of a given signal or a time series. The method processes a series of airway cross-sectional area changes between baseline (first) and follow-up (second) CT scans, and generates the posterior distribution over multiple possible points of abnormal variations, whose the mode is taken as the final prediction. We test the efficacy of the method on (1) CT images of healthy airways with simulated dilatation, and (2) pairs of real images of IPF-affected airways acquired approximately 1 year intervals. For the simulated dataset, we measured the detection accuracy with respect to the commonly used naive thresholding and maximum likelihood based methods [9]. For the longitudinal IPF dataset, we compare the predictions of our model to measurements from radiologists based on two different protocols.

2 Method

In this section, we introduce our method for quantifying progression of IPF on a patients serial CT scans. The method proceeds as follows. Firstly, we fit a tubular shape model to the airways in both baseline and follow-up CT scans, and acquire estimates of the cross-sectional areas along them (Sec. 2.1). We then treat the difference in the cross-sectional areas as a 1D signal, and employ the proposed Bayesian changepoint model to estimate the posterior distribution over locations of abrupt airway dilatation (Sec. 2.2 and 2.3). Lastly, we post-process this posterior distribution to determine the region of dilatation (Sec. 2.4).

2.1 Airway Pre-processing

In this work, for each airway track, we acquire a series of cross-sectional area measurements using the method proposed by Quan et al [12] (Fig. 1). Following a semi-manual segmentation of the airway, the method computes the airway centreline and the corresponding normal planes at contiguous intervals, each of which is then used to estimate the cross-sectional area. The final output is a 1D function of cross-sectional area along the airways for baseline $f_{B}(x)$ and the follow-up $f_{F}(x)$ scans.

The next task was to align the signals on both baseline $f_{B}(x)$ and follow-up $f_{F}(x)$ scans. To this end, we resample the signal to 1mm using cubic interpolation. We considered the first 50 points on both signals $g_{B},g_{F}$ from the start of the carina. To register the two airways together, we apply the transformation $f_{F}(x-a)$ where

[TABLE]

The longest of the two signal were truncated from the right hand side such that both signals were of the same length. For the rest of the methodology, we will only consider the series difference defined as $\boldsymbol{y}=\log(f_{F})-\log(f_{B})$.

2.2 Bayesian Changepoint Model

The progression of fibrosing lung diseases can manifest itself as dilatation of the internal airways. Thus, we hypothesise that at the start of dilatation, the series $\boldsymbol{y}$ undergoes an abrupt variation, which we refer to as a changepoint. More formally, given signal $\boldsymbol{y}=(y_{1},\dots,y_{n})$ of length $n$, we define a changepoint $\tau$ as the location where there exists a change in parameters $\theta$ in the underlying distribution $F$. In other words, at changepoint $1<\tau<n$, the observations $\boldsymbol{y}$ can be separated at $\tau$ such that:

[TABLE]

where $\theta_{1}\neq\theta_{2}$. This definition can be naturally extended to the scenario with $M$ changepoints; we denote the changepoint location vector by $\tau=(\tau_{1},\dots,\tau_{k})$, with parameters $\theta=(\theta_{1},\dots,\theta_{k+1})$ for each respective segment. For ease of notation, we also denote $\tau_{0}=1$ and $\tau_{k+1}=n$. Assuming statistical independence between segments, the likelihood factorises as:

[TABLE]

where $y_{\tau_{l-1}:\tau_{l}}=(y_{\tau_{l-1}},\dots,y_{\tau_{l}})$ We also specify prior distributions on the the number of changepoints $p(M;\delta)$, the locations of the changepoints $p(\tau|M;\gamma)$, and the parameters of the corresponding segments $p(\theta|M;\beta)$ where $\beta$, $\gamma$ and $\delta$ represent the hyper-parameters.

Given the likelihood and the prior distributions above, we would like to estimate the posterior distributions over the number of changepoints $p(M|\boldsymbol{y})$, the locations of changepoints $p(\tau|\boldsymbol{y})$ and the parameters of the respective segments $p(\theta|\boldsymbol{y})$. Commonly in Bayesian changepoint analysis, likelihoods with a conjugate prior are chosen to allow calculation of analytical posteriors for changepoints and parameters, enabling faster convergence and better accuracy [4]. However, in our work, we wish to relax this conjugacy assumption in order to allow flexibility in the choices of prior and likelihoods. In the next section, we describe the inference scheme in such highly general setting, viable for arbitrary forms of the likelihood distribution $F$ and the prior. In this paper, $F$ was defined as the Student t distribution with parameters $\theta=(\nu,\mu,\sigma^{2})$ with degrees of freedom $\nu$, mean $\mu$ and variance $\sigma^{2}$. We chose t-distribution for its robustness to outliers [11] caused by noise from the area measurements.

2.3 Reversible Jump MCMC for Posterior Inference

Posterior inference with our model possesses two challenges. Firstly, without the conjugacy assumption, computing the posterior distributions is intractable. Secondly, the dimensionality of the posterior distribution over the changepoints $\tau$ is given by $M$ and varies during inference. To combat the first problem, we use the Metropolis-Hasting (MH) algorithm [2], a variant of Markov Chain Monte Carlo (MCMC) methods that can sample from the posterior, with or without conjugacy. Given that the number of changepoints $M$ is known, MH can be used to sample from the posterior distributions over the changepoints $\tau$ and segment parameters $\theta$. To address the second problem of varying posterior dimensionality $M$, we extend the above sampling scheme to the Reversible Jump MCMC framework [5]. Taken all together, the method is capable of traversing the full posterior distributions for $M,\tau,\theta$ and we refer to this as Reversible Jump Metropolis Hasting (RJMH) algorithm.

Overview of RJMH:

The RJMH proceeds by randomly executing one of four possible moves at each iteration (see Algorithm 1 for the pseudo-code):

The parameters of every segment is updated. 2. 2.

One of the changepoints is selected randomly and updated, and the parameters of its neighbouring two segments are updated. 3. 3.

Randomly choose a location in the time series and add it as a changepoint. The parameters of the resultant two new segments are also sampled. 4. 4.

One of the changepoints is randomly removed and the parameters of the new segment are sampled

We denote these moves as: $\curlyvee_{\theta},\curlyvee_{\tau},\curlyvee_{M\rightarrow M+1},\curlyvee_{M+1\rightarrow M}$ respectively. We also define the maximum number of changepoints $k_{max}$ and at the boundary cases for $k$, we impose restrictions such that $\curlyvee_{M+1\rightarrow M}$ is not possible when $k=0$ and $\curlyvee_{M\rightarrow M+1}$ is not possible when $k=k_{max}$. Each move updates the appropriate subset of parameters $\theta,\tau$ by sampling from the corresponding proposal distributions $q(\theta_{new}|\theta_{old})$ and $q(\tau_{new}|\tau_{old})$, and is only executed if it passes the associated acceptance criteria $\alpha$.

We can subdivide these four moves into two sets: fixed $M$ (MH Steps) and changing $M$ (RJ Steps) i.e. resampling parameters and moving changepoints does not affect the changepoint model $M$ whereas adding or deleting a changepoint changes the value of $M$.

MH Steps:

For $\curlyvee_{\theta}$, we set $q(\theta_{new}|\theta_{old})=(\mu_{old},\sigma_{old}^{2},\nu_{old})+$N$(0,\epsilon)$. This step resamples parameters of each segment by proposing Gaussian perturbations around the current values of parameters for all segments.

For $\curlyvee_{\tau}$, we set $q(\tau_{new}|\tau_{old})=\tau_{old}+(-1)^{b}$Poi$(\lambda)$ where $b\sim$Binary$[0,1]$. This step selects a changepoint $\tau$ at random and shifts it with a Poisson perturbation. The segments neighbouring this new changepoint location have parameters $\theta$ resampled as in move $\curlyvee_{\theta}$ using the current segment parameters.

Note that each proposal in the MH steps will be symmetrical and therefore, cancel out with the reverse proposal term (See Algorithm 1).

RJ Steps:

For the changing $M$ moves, the algorithm moves in dimensionality and we need to introduce alternative proposal terms for the new segments produced as there are no current segment parameters to base our proposals on.

For $\curlyvee_{M\rightarrow M+1}$ i.e. $M_{new}=M+1$, we proposed random new changepoints over our data, $\tau_{new}\sim$U$[1,n-1]$. The proposed $\tau_{new}$ split an existing segment into a new left segment $\theta_{l}=(\mu_{l},\sigma^{2}_{l},\nu_{l})$ and new right segment $\theta_{r}=(\mu_{r},\sigma^{2}_{r},\nu_{r})$. Our proposal for $\mu_{i},\sigma^{2}_{i}$ are defined by a Gaussian perturbation on empirical values of the respective $i=l,r$ segments (Fig. 2). The proposal for $\nu_{i}$, is Gaussian perturbation of the previous update $\nu_{old}$. Due to dependence of the $\nu_{i}$ proposal, a Jacobian term is introduced $|J_{M\rightarrow M+1}|=2$.

For $\curlyvee_{M+1\rightarrow M}$ i.e. $M_{new}=M-1$, we remove a changepoint $\tau_{new}$. As before, proposals for $\mu_{new}$, $\sigma^{2}_{new}$ are defined using empirical values of the segments and Gaussian perturbation (Fig. 2). The proposal for $\nu_{new}$ is the mean of the previous $\nu_{l},\nu_{r}$. The move introduces Jacobian term $|J_{M+1\rightarrow M}|=0.5$.

Initialization:

The overall algorithm structure can be seen in the pseudo-code Algorithm 1. With the model defined, we mentioned the priors used and address some common implementation techniques. The priors used are given as follows:

[TABLE]

The hyper-parameters for $\mu,\sigma^{2},\nu$ were chosen to be non-informative and within plausible ranges. In order to reduce the number of changepoints we could detect (in the case of acquisition noise from the CT), we set the expectation for $M$ to be sufficiently low. In terms of implementation, we follow a standard procedure of setting a burn-in for the number of iterations to ignore any potential issues with intialisation at approximately 25% of the total iteration count. To remove any chance of auto-correlation, we thin the number of samples by only storing the $5^{\textup{th}}$ iteration, after the burn-in period. Additionally, we choose the number of iterations, $N$, to be sufficiently large in order to achieve convergence.

2.4 Locating airway dilatation

For airways affected by IPF, dilatation starts at the distal point and progresses in the proximal direction [7]. Therefore, we topologically can assume that each affected airway undergoes a single changepoint from which dilatation starts. To locate such unique changepoint, we consider the posterior probability of the changepoint $p(\tau|\boldsymbol{y})$, and perform the following post-processing steps.

On each airway track, the proximal region is surrounded by cartilage [14]. As the airway track loses cartilage support, the geometry of the airway changes and results in a changepoint. Since such biological changepoint is independent of the disease state and occurs prior to dilatation, we eliminate it by discounting the most proximal peak in the posterior distribution $p(\tau|\boldsymbol{y})$. We then selected highest peak on the modified posterior $p(\tau|\boldsymbol{y})$ as the final estimation for the starting point of dilatation.

3 Evaluation & Results

We evaluated our proposed method with two experiments. Firstly, we used a set of healthy airways, in which we had simulated dilatation to quantify the accuracy of our method and compared the results with conventional changepoint methodologies. Secondly, we assessed the clinical utility of our method on a dataset of airways affected by IPF. We compared our labels with those from two experienced thoracic radiologists. The properties of the images used in both experiments are displayed on Tab. 1.

3.1 Disease simulation

To quantitatively assess accuracy, a ground truth is required. To this end, we applied our point detection algorithm on augmented healthy airway series to simulate the airway dilatation caused by IPF. After obtaining written informed consent from 3 patients, a trained radiologist (R1) selected 14 pairs of healthy airways in both baseline and follow-up scans. The image properties are displayed on Tab. 1. They were acquired from different scanners and used different reconstruction kernels. The airways were pre-processed as described in Sec. 2.1 to produce a function of area change along the length of the airway. We interpret this output signal as the overall noise caused from normal biological changes and acquisition.

We modelled the change in dilatation with a logistic function;

[TABLE]

We interpreted $M$ as the magnitude of dilatation and $\alpha$ as the point of dilatation. The parameters $\alpha$ are set such that the dilatation starts 10-40mm from the distal point in increments of 5mm. In addition, we set $M$ to range from 0.3-3 in increments of 0.25. Finally we set $k=0.5$mm*-1*, in order to create an abnormal increase in cross sectional area. To simulate the dilatation on the airway; the logistic function was added to the area change of the healthy airway, as shown on Fig. 3. We applied every permutation of $M$ and $\alpha$ on each of the 14 healthy airways. After the signal augmentation, we applied our proposed Bayesian changepoint detection algorithm.

The proposed method was compared against two conventional methods. First, we use a basic thresholding method. For a given signal, we applied a 5mm moving mean average and thresholded the point at which the signal reached above the upper quartile from the right hand side. Secondly, we implemented the method based on Lavielle [9], in summary we consider $K$ changepoint and these changepoints $y_{i}$, minimize the function:

[TABLE]

where $\beta$ is modified such that the function finds less than $K$ changepoints. For our paper, different $\Delta$ were evaluated and we found

[TABLE]

gives the most accurate results. To replicate the post processing of our proposed method, we consider $K=2$ possible changepoints. This takes into account the changepoint caused by the support cartilage. Finally, we set a minimum distance of 20mm. Once we acquired the changepoints, the most peripheral point was chosen as the point of dilatation. The implementation was performed through Matlab inbuilt function; findchangepts111https://www.mathworks.com/help/signal/ref/findchangepts.html last accessed on ..

Fig. 4 shows the accuracy for each individual method as a heatmap. The metric used to quantify accuracy was the displacement between the ground truth and the changepoint given by the proposed method. A positive displacement (mm) corresponds to an overestimation of the ground truth towards the distal point. Each entry on the heatmap corresponds to the median average of all displacements over 14 airway pairs.

When the magnitude of dilation is larger than $M>0.75$, our proposed method achieves consistently higher accuracy than Lavielle [9]. We note that the accuracy gain in the peripheral regions of the airways at $\alpha=$ 10-30mm from the distal point are the most clinically relevant in IPF as lung parenchymal damage begins in the lung periphery and progresses proximally [7]. Furthermore on the same peripheral regions $\alpha=$ 10-30mm, the baseline method showed systematic bias in accuracy towards the central airways. This was due to the baseline method being influenced by outliers from the longer expanses of normal airway regions. The proposed method uses the t-distribution as the likelihood thus making it robust to possible outliers within the data [11]. On the other hand, the proposed method suffers from poor accuracy below magnitudes of dilatation $M=0.75$. However, in physical terms a dilatation of $M=0.75$ corresponds to a percentage increase in cross sectional area of $e^{0.75}-1\approx 112\%$. This is within the range of normal biological change of the airways [3].

3.2 Application to Airways affected by IPF

We applied our method to airways affected by IPF. The purpose was to compare our measurement to the labels provided by a radiologist and compute the volume change of the identified diseased airway regions. For our dataset, we acquired 4 airway pairs from 2 patients after obtained a waiver for consent from the local Research Ethics Committee. All airways were judged by the radiologist R1 to be dilated as a consequence of IPF on baseline and to have visually worsened on follow-up imaging. Image properties are displayed in Tab. 1.

We compared the performance of our method against two trained thoracic subspecialist radiologists. Two radiologists R1, R2 identified the point of at which a given airway was seen to demonstrate increased dilatation on the follow-up CT scan. To assessed the reproducibility of manual labeling, each radiologist labelled the same airway twice through two different protocols. In the first method, the radiologists interrogated axial CT images. Using 2 separate workstations and the airway centreline, the radiologists identified the point on the centreline (on the follow-up scan) where the airway demonstrated definitive worsened dilatation. For the second method, the radiologist compared the aligned reconstructed cross-sectional planes on baseline and follow-up scans. The radiologist then selected the slice where the airway had worsened when evaluated against the baseline scan. An example of both protocols are displayed on Fig. 5.

Fig. 6 compares the predictions of our method with the labels from radiologists obtained in two different protocols. The results indicate that the predictions for Airway 1, 3 and 4 are within the range of the radiologists’ labels. In the case of Airway 2, although our method based on the maximum peak overestimates with respect to the radiologists’ predictions, the posterior distribution contains another equally probable peak that underestimates the radiologists’ labels (see the second highest peak at 70mm), potentially indicating a more proximal point of dilatation. To test this, we delineated the boundary of the lumen on the reconstructed cross sectional slices at baseline in the neighbouhood of this peak, and Fig. 7 shows the initial few slices (62-64mm). We observed that, when the delineated boundary from baseline was superimposed on the follow-up scan, the boundary is contained inside of the follow-up lumen with several lumen pixels are consistently outside from the boundary. Thus, this result indicates that the starting point of dilatation is more proximal than the labels from the radiologists.

To demonstrate the clinical utility of locating the starting point of airway dilatation, we compared longitudinal airway volume changes in diseased and healthy regions of each of the airway tracks. To find the volume of the airway, we considered the aligned signals, $f_{F}$,$f_{B}$. These signals are measurements of area against the airway arc length. Thus volume can be computed via the area under the curve. For our work, we used the trapezium rule to find the volume. Three volumetric regions were considered: (i) the entire airway $V_{c\rightarrow d}$ i.e. from the carina, $c$ to the distal point, $d$, (ii) carina to the dilatation point $t$, denoted as $V_{c\rightarrow t}$. (iii) the dilatation point to the distal point, denoted by $V_{t\rightarrow d}$. Note that $V_{c\rightarrow d}=V_{c\rightarrow t}\cup V_{t\rightarrow d}$ and $V_{c\rightarrow t}$ does not overlap with $V_{t\rightarrow d}$. Tab. 2, shows the results of the percentage volume change. The volume change in $V_{t\rightarrow d}$ had greater sensitivity for selecting progressive airway dilatation in IPF than the volume change in the entire airway $V_{c\rightarrow d}$.

4 Discussion & Conclusion

In this paper, we modelled changes in area along airway tracks as a time series with the purpose of detecting abnormal dilatation caused by IPF using Bayesian changepoint detection. Experiments on simulated data show that our model is able to detect the starting location of airway dilatation with superior accuracy than the relevant baseline methods. The results on the IPF longitudinal dataset display reasonable agreement with radiologists, while in one case indicating a more plausible location of dilatation, potentially missed by the experts.

Identifying changepoints and thereby calculating a change in airway dilatation over time could become a sensitive measure of IPF aggravation. This would form an important secondary endpoint in drug trials, where our measurements could indicate whether a new drug ameliorates disease progression better than existing medications. In the future, we hope to evaluate such utility of our proposed method on a larger cohort of IPF subjects. Furthermore, our method is applicable to other lung airway diseases, characterised by dilatation (such as cystic fibrosis) or other geometrical deformations, and such extensions remain valuable future work.

5 Acknowledgements

Kin Quan is supported by the EPSRC-funded UCL Centre for Doctoral Training in Medical Imaging (EP/L016478/1) and the Department of Health NIHR-funded Biomedical Research Centre at University College London Hospitals. Ryutaro Tanno is supported by Microsoft Research Scholarship. Micheal Duong is funded by the Centre for Doctoral Training in Financial Computing and Analytics. Arjun Nair is funded by National Institute for Health Research Biomedical Research Centre. Mark Jones and Christopher Brereton are supported by National Institute for Health Research. Joseph Jacob is a recipient of Wellcome Trust Clinical Research Career Development Fellowship 209553/Z/17/Z.

An adapted version of this manuscript was accepted to The 10th International Workshop on Machine Learning in Medical Imaging (MLMI 2019). In conjunction with MICCAI 2019, Shenzhen, China. It is cited as follows: Quan K. et al. (2019) Modelling Airway Geometry as Stock Market Data Using Bayesian Changepoint Detection. In: Suk HI., Liu M., Yan P., Lian C. (eds) Machine Learning in Medical Imaging. MLMI 2019. Lecture Notes in Computer Science, vol 11861. Springer, Cham. DOI: https://doi.org/10.1007/978-3-030-32692-0_40.