Spectral graph theory efficiently characterizes ventilation heterogeneity in lung airway networks

This paper introduces a linear operator for the purposes of quantifying the spectral properties of transport within resistive trees, such as airflow in lung airway networks. The operator, which we call the Maury matrix, acts only on the terminal nodes of the tree and is equivalent to the adjacency matrix of a complete graph summarizing the relationships between all pairs of terminal nodes. We show that the eigenmodes of the Maury operator have a direct physical interpretation as the relaxation, or resistive, modes of the network. We apply these findings to both idealized and image-based models of ventilation in lung airway trees and show that the spectral properties of the Maury matrix characterize the flow asymmetry in these networks more concisely than the Laplacian modes, and that eigenvector centrality in the Maury spectrum is closely related to the phenomenon of ventilation heterogeneity caused by airway narrowing or obstruction. This method has applications in dimensionality reduction in simulations of lung mechanics, as well as for characterization of models of the airway tree derived from medical images.

This paper introduces a linear operator for the purposes of quantifying the spectral properties of transport within resistive trees, such as airflow in lung airway networks. The operator, which we call the Maury matrix, acts only on the terminal nodes of the tree and is equivalent to the adjacency matrix of a complete graph summarizing the relationships between all pairs of terminal nodes. We show that the eigenmodes of the Maury operator have a direct physical interpretation as the relaxation, or resistive, modes of the network. We apply these findings to both idealized and image-based models of ventilation in lung airway trees and show that the spectral properties of the Maury matrix characterize the flow asymmetry in these networks more concisely than the Laplacian modes, and that eigenvector centrality in the Maury spectrum is closely related to the phenomenon of ventilation heterogeneity caused by airway narrowing or obstruction. This method has applications in dimensionality reduction in simulations of lung mechanics, as well as for characterization of models of the airway tree derived from medical images.

Background
In healthy human lungs, the airways form a bifurcating tree where, on average, around the first 16 generations of airways are purely conductive and serve to transport gas from the mouth to the alveolar region where the majority of gas exchange takes place. The conducting airways terminate in approximately 30 000 respiratory units (or acini) [1], where diffusion becomes the dominant transport mechanism. Resistive flow in the conducting airways, as well as tissue compliance, determine the ventilation mechanics and hence require a high-dimensional mathematical model to be represented accurately [2][3][4]. Efficient dimensionality reduction algorithms are needed to improve the computation time and tractability of using such models in inverse problems in the future.
Furthermore, understanding the relationship between structure and function in the lungs is crucial to identifying lung disease physiology and its progression. One particular marker of progression in obstructive lung conditions is 'ventilation heterogeneity' (VH), or the uneven distribution of fresh gas in the lungs, which can result from obstructions in the airways [5][6][7]. The methods presented here provide a new approach to characterize the resistance structure and its effect on VH. This will be useful for visualization and characterization of complex airway models based on patient computed tomography (CT) images, lung casts or micro-CT scans of excised lungs.
Graph-theoretic perspectives have proven useful in characterizing, classifying and understanding transport in various examples of biological networks [8][9][10]. The approach we take in this paper uses concepts from spectral graph theory, where a linear operator describing a particular physical property or characteristic of the network is decomposed into its eigenspectum. This has proved useful for dimensionality reduction and visualization of large datasets in various areas of science [11] and has a number of realworld applications [12] including diffusive transport networks [13], but has not previously been applied to the airways. Some standard operators, such as the adjacency and Laplacian operators, are generic to all networks and can be adapted to characterize different processes or uncover important motifs in the network [14]. However, for lung airway networks and other resistive trees, these operators are not necessarily the optimal representation of the physical processes of interest on the network, as we show here.
In this paper, we demonstrate that the resistive properties of realistic tree networks are characterized by a tree-specific operator that we call the Maury matrix, first introduced in [15]. This matrix provides a complete description of the linear resistance relations on the network. In §2, we formulate the flow problem on an airway network in terms of the conductance Laplacian and the Maury matrix, and show how each can be approximated via spectral reduction. We show how VH is evaluated for a set of realistic airway network models (four of which are based on CT imaging). In §3, we compare spectral approximations of the Laplacian and Maury operators, and show in particular how the reduced Maury operator efficiently captures spatial patterns of VH, giving a new method for dimensionality reduction in these systems.

Methods
Unlike several other examples of biological transport networks, mammalian airways develop robustly into tree networks [1], such that they contain no cycles and airways can be defined hierarchically by parent airways branching into children airways. Therefore, we model a tree network N ¼ {V, E} consisting of a set of nodes (or vertices), V ¼ {v i }, and a set of branches (or edges) representing airways, We define each branch's orientation to point from its more proximal node v i to the more distal node v i 0 based on its graph distance from the root node v 1 .
The set of terminal nodes is denoted T , V and we order the indices i so that the terminal nodes are listed last T ¼ {v jVjÀjT jþ1 , . . . , v jVj }, where |.| indicates the number of objects in the set. The tree topology means that every branch e j has a unique distal node v i 0 so we index the branches such that j = i 0 − 1. In the conducting airways of the lung, the root node represents the entrance from the upper airway, which is connected to the trachea that we label as branch e 1 . See figure 1a for an overview of this notation. The distribution of flow through the airway network in the lung is primarily dependent on the interplay between resistance and compliance of the airways and tissue. In this paper, we focus on new ways to model and analyse the resistance of the airway network, and its implications for the distribution of gas flow.

The conductance Laplacian
The network incidence matrix N [ Z jVjÂjEj maps the nodes of the network to their associated branches. Its entries are 8 < : (2:1) We model flow through the conducting airways of the lung by the linear resistance equation where q [ R jEj is the vector of branch fluxes, r [ R jEj is the vector of branch resistances and P [ R jVj is the vector of node pressures. The incidence matrix N is a discrete boundary operator, mapping 1-chains (branch quantities) to 0-chains (node quantities), and is the discrete analogue of the divergence operator in vector calculus [16]. Similarly, N T is analogous to the gradient operator and performs the reverse mapping. Applying N to (2.2) gives where L ¼ N diag(r) À1 N T is the network conductance Laplacian. The node flux Q ¼ Nq for each node v i is equal to the sum of fluxes into node i from the branches connected to it. Assuming incompressibility, these are zero except on any nodes that are sources or sinks. In the lung, gas can only enter or exit the conducting airways via the upper airway or the transitional bronchioles (which feed the acini), so henceforth, we assume that the root v 1 and terminal nodes v i [ T are the only sources T e 1 r 1 + r 3 + r 6 r 1 + r 2 + r 4 r 1 + r 2 + r 5 r 1 + r 3 + r 7 r 1 +  royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 17: 20200253 or sinks. Therefore, the vectors of node and branch quantities can be written as (2:4) Here, the separate blocks of P and Q distinguish the quantities defined on the root node v 1 , internal nodes V int ¼ {v 2 , . . . , v jVjÀjT j }, and terminal nodes T ¼ {v jVjÀjT jþ1 , . . . , v jVj }. Likewise, q and r are separated into quantities defined on the root branch e 1 , the internal branches E int ¼ {e 2 , . . . , e jVjÀjT jÀ1 } and terminal branches E term ¼ {e jVjÀjT j , . . . , e jVjÀ1 }. In block form, the incidence and Laplacian operators are Here, a = (1, 0, …, 0) T and I [ R jT jÂjEtermj is the identity matrix.
denotes the all-one vector of size M, and so e jVj is the vector identifying all nodes. Thus, e jVj is a zero-eigenvalue mode of L, and multiplying (2.3) by e T jVj gives the condition for global mass conservation Q 1 ¼ À P vi[T Q i . In fully connected networks, as considered here, this is the only zero-eigenvalue mode and so the conductance Laplacian L is rank jVj À 1. It follows that jT j þ 1 boundary conditions are required to define the jVj þ jT j unknowns in (2.3) (jVj node pressures and jT j terminal fluxes).
In the context of lungs, we are generally interested in pressure boundary conditions on P 1 and P term . To remove the arbitrary pressure constant, we first reformulate in terms of the pressure drop relative to the entry node v 1 such that DP ¼ P 1 e jVj À P, and so ΔP = (0, ΔP int , ΔP term ) T . Thus, given ΔP term and inserting (2.4) into (2.3), the system of equations to be solved is and Equation (2.6) can be solved efficiently for ΔP int by numerical linear algebra methods due to the sparsity of the operator L int . In general, for connected networks, L int is invertible and so substituting (2.6) into (2.7) gives the terminal fluxes in terms of the pressure boundary conditions In practice, we solve this system of equations by implementing the 'SparseLU' method from the C++ library Eigen [17], based on LU factorization as detailed in electronic supplementary material, text S1.1.

Associated spectra of the conductance Laplacian
The solution to (2.6) can also be expressed in terms of the eigenbasis of the operator L int . In general, the internal Laplacian L int has jVj À jT j À 1 non-zero eigenvalues that we label l 1 l 2 . . . l jVjÀjT jÀ1 with corresponding eigenvectorŝ u 1 , …,û jVjÀjT jÀ1 , where the hat indicates a normalized vector such thatû T kûk 0 ¼ d kk 0 and δ is the Kronecker delta. Substituting the eigen-decomposition of L int into (2.6) and solving for ΔP int gives Substituting this solution into (2.7) gives the corresponding solution for Q term . Dimensionality reduction can in principle be achieved using spectral methods by reconstructing the system with a subset of the eigenmodes. We can measure the relative dominance of each mode in the decomposition in (2.9) by the size of the coefficient jp k j ¼ jû T k BG term DP term =l k j. However, the relative weight of these modes is dependent on the boundary condition ΔP term . A more general measure of how well a subset of the eigenmodes captures the behaviour of the whole operator, independent of boundary conditions, is to consider how closely these modes approximate the inverse of the operator, L int , which we measure by the normalized distance δ L [11] given by where ||·|| F indicates the Frobenius norm and M [ {1, . . . , jVj À jT j À 1} is the number of modes used in the approximation. Evidently, the modes with the smallest eigenvalues contribute the most to this reconstruction, however depending on the boundary conditions used they may not be the most dominant in the reconstruction of the solutions ΔP int and Q term .

The Maury matrix
An alternative representation of the resistive airway tree network is presented in [15] for the case of a symmetric dyadic tree. Here, we generalize this operator to all connected trees and name it the Maury matrix in reference to [15]. To derive this operator, we begin by noting that connected subtrees of the airway network can be defined due to the hierarchical structure of trees. We define the matrix S ¼ (s 1 , . . . , s jEj ) where each column s j maps the branch e j to a subset of nodes V j , {v 2 , . . . , v jVj } descended from it such that The first column of S is s 1 ¼ e jVjÀ1 as all nodes except v 1 are descended from the root branch. The column vectors s j satisfy the relation . , e jn } are the n child branches of e j , and i j is the jth column of the identity matrix I [ R (jVjÀ1)Â(jVjÀ1) . This relation follows from noting that each sub-tree E j is a union of its child subtrees and its connecting node, such that V j ¼ {v jþ1 , V j1 , . . . , V jn }. It follows that SN T ¼ (e jVjÀ1 , À I) and therefore where the rows of T [ R jT jÂjEj are comprised of the final jT j rows of the matrix S. Substituting this into (2.2) and multiplying both sides by ÀS diag(r) gives royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 17: 20200253 3 Comparing (2.13) to (2.3), we see that S T diag(r)S ; L À1 0 where L À1 0 is the conductance Laplacian L with the first row and column removed (or the reduced Laplacian, as features in the matrix-tree theorem [18]). Finally, the bottom jT j rows of (2.13) directly relate the terminal pressure drops to the fluxes where R ¼ T diag(r)T T is the aforementioned Maury operator. Comparing with (2.8) gives The entries of the terminal node map T ij are 1 if the terminal node v 1þjVj int þi is descended from branch e j and zero otherwise (see (2.11)). We denote the set of subset of terminal nodes descended from e j as T j , T (figure 1a). The Maury matrix is symmetric and in index notation its components are [19] between the node v 1 and node v iLCA such that v iLCA is the lowest common ancestor (LCA) node of the terminal nodes v i and v i 0 . A direct solve of (2.14) is an inefficient method for solving the system of equations (compared to (2.6) and (2.7)) because the matrix R has, in general, jT j 2 non-zero entries and so cannot benefit from any algorithms optimized for the factorization of sparse matrices. Nonetheless, this does not inhibit computation of the spectra of R because we can decompose R into sparse, non-square matrices. Details of the numerical methods used are given in electronic supplementary material, text S1.2.

Spectral properties of the Maury matrix
Diagonalizing the Maury matrix R into its normalized eigenbasis and substituting into (2.14) gives for each eigenvalue μ k and associated eigenvectorv k . The matrix R has the form of an adjacency matrix for a complete network of resistors (with self-loops) relating the terminal nodes (as shown figure 1b). Therefore, as shown in (2.17), the modes act as resistors in parallel with resistance μ k , each subject to the pressure drop DP T term Áv k . The matrix R is non-negative, regular and symmetric and so by the Perron-Frobenius theorem its largest eigenvalue mode μ 1 is unique and the eigenvalues of the Maury matrix are all positive. Thus, we order the indices by m 1 ! m 2 ! . . . ! m jT j (the opposite ordering to modes of the Laplacian).
From (2.17), the solution for the terminal node fluxes is and so an approximation to Q (approx) term is given by limiting the sum in (2.18) to some subset of the mode spectrum. The magnitude of the coefficient jq k j ¼ jDP T termvk =m k j quantifies the relative importance of each term in the expansion of (2.18). As for the internal Laplacian operator, we measure the accuracy of such an approximation for arbitrary boundary conditions by the convergence of the inverse of R such that where M is the number of (smallest eigenvalue) modes used in the expansion. We will compare this with δ L (M) for specific networks to compare the general efficiency of dimensionality reduction for these two operators.

Dimensionality reduction in simulations of ventilation in the lungs
In models of ventilation in the lungs, the effects of lung tissue compliance generally dominate the mechanics compared to airway resistance, except in cases where there is airway narrowing or blockage [20]. Consider the simplest model of ventilation where each terminal node v i [ T is connected to an elastic bag of volume V i (t) and elastance κ with all bags subject to the same pleural pressure P pl (t) driving the breathing motion at time t.
The terminal fluxes are related to the bag volumes via Q term (t) ¼ À _ V(t) and the terminal pressures are P term ¼ [P 1 þ P pl (t)]e jT j þ kV(t). Then (2.14) becomes (2:20) We are interested in characterizing how the distribution of ventilation is affected by the airway resistance, and so have assumed that the lung unit compliances are all equal to κ. Diagonalizing R into its eigenbasis in (2.20) gives Thus each mode of R is characterized by an independent ODE with a resistance term μ k and elastance κ driven by a pressure term proportional to e T jT jv k . Solving (2.21) for a k (t) gives (2:22) The first term in the solution is the elastic contribution from all units, while the second term is the resistive contribution. The sum in the resistive term can be approximated by the largest eigenvalue modes with a cut-off of the order μ k /(κτ) ∼ O(1) (where τ is the typical breath timescale), since modes satisfying μ k /(κτ) ≪ 1, will be dominated by lung unit compliance (the first term in (2.22)).
In the simulations used in this paper, we take the pleural pressure to be sinusoidal such that P pl = P pl0 + P s sin(2πt/τ). Substituting into (2.22), we see that the periodic solution (when t ≫ μ 1 /κ) is (2:23) We define the ventilation of each terminal unit as , where the maximum and minimum are taken over a single breath cycle 0 ≤ t < τ, which gives : (2:24) Inefficient ventilation of the lungs can occur if there are blockages or narrowing of the airway lumen, which can be caused by numerous pathophysiological factors, resulting in increased resistance of those airways and inhomogeneous delivery of gas to the lung units, commonly known as ventilation heterogeneity (VH). To characterize this in simulations, we royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 17: 20200253 measure ventilation heterogeneity via the coefficient of variation of the ventilation σ V across all terminal units, (2:25) Note that the assumption of a sinusoidal pressure profile is used in this paper as it provides analytical solutions, however the above can be extended to generic pressure profiles simply by using numerical integration to estimate the value of the integral in (2.22). In many cases, it would be desirable to set the pressure boundary condition to achieve a particular flow rate at the mouth. This results in a system of equations to be solved for the weights a k (t) and P pl (t), which again can be approximated by reducing to a subset of the largest eigenvalue modes. Therefore, the above decomposition gives a method to approximate the ventilation efficiency of an airway network given κ and τ by truncating the sums in (2.24) to M , jT j modes of the Maury matrix. While this method is specifically derived for sinusoidal breathing, it provides a useful qualitative measure of relative efficiency that can be used for comparing different networks and is easily computed using a subset of the full spectrum of the Maury operator.

Airway models
In this paper, we use idealized models of airway geometry and models derived from CT images, which we now briefly describe.

Asymmetric Weibel branching
First, we adopt the Weibel-type model used in [21,22] to systematically test the effect of asymmetry in resistance on the spectral properties of the Maury and Laplace operators. This idealized model describes a dyadic tree where each airway branches into major and minor children such that where a and l are the radii and lengths, respectively, of the parent and (major and minor) child branches as labelled by the subscripts. The parameter 0 ≤ A < 1 quantifies the asymmetry in branching (A = (1 − 2r) in terms of the notation in [21]). The cube root results in the two daughter airways having the same combined effective Poiseuille resistance and volume as the symmetric (A = 0) case. In this model, we assume a lengthto-diameter ratio of 3 for all airways and a fixed number of divisions N between the root branch and the terminal branches (resulting in N + 1 generations of airways).

Horsfield model 1
In real lungs, branches terminate in a wide range of generations.
To characterize this effect we use the Horsfield model 1 [23] derived from cast measurements of an adult male, with airways isotropically scaled to 60% of their original volume to a lung with approximately 3 l functional residual capacity (approximately that of the average adult). The model is terminated at Horsfield order 4 (as defined in [23]) because this corresponds to the transitional bronchioles as defined by Weibel [1] and gives 29 240 terminal airways, approximately the number of acini in the lung [24]. The number of bifurcations from trachea to terminal bronchiole ranges from 10 to 25, this is achieved by imposing that parent airways branch into daughters with different 'Horsfield orders' (generations counting up from the terminal nodes). The radii and lengths are Gaussian random variables with mean taken from the data in [23] and a standard deviation of either 10 or 20% of the mean to approximate the variation in healthy lungs, which we label H10 and H20 respectively.

Image-based (CT) models
Finally, we also use four image-based (IB) models derived from CT images. These use volume-filling branching algorithms [25,26] to populate the distal generations of airways in each lobe of the lung, which cannot be resolved from CT. The radius of these distal airways is assumed from an exponential fit to lung cast data, as in [26], such that where D(g) is the airway diameter, g is the Horsfield order of the airway, N is the Horsfield order of the most distal ancestor airway taken from CT and its diameter D N . The parameter R d H is taken to be 1.15 as in [26]. The four IB models used here (labelled IB1-4) are from CT scans of children with cystic fibrosis and normal lung function (measured by spirometry) from a previously published study [27] (IB2 is shown in figures 4 and 5 below). The number of airways successfully acquired from CT depends on the image resolution, and was approximately between three and five generations of bifurcations from the trachea, and was acquired using the open-source software Pulmonary Toolkit [28]. The algorithm to generate these models used the C++ library stl_reader [29] and is detailed in electronic supplementary material, text S1.3. As a simple first approximation, the branch resistance is assumed to follow Poiseuille's Law and so r = 8πνl/a 4 where ν is the viscosity of air. To compute the spectra of the operators we use the Eigen [17] and Spectra [30] libraries in C++. The source codes of the algorithms used are available in online repositories [31,32].

Results
In the following section, we first study the spectra of the Maury operator for a number of idealized and image-based airway models and compare with the spectra of the internal Laplacian operator. Then, in the final results section, we demonstrate how (2.24) can be used to produce a lowdimensional approximation of simulations and characterize realistic lung networks.

Comparison of Maury and internal Laplacian
operator mode structures Figure 2 shows the relative accuracy of the inverse approximation for each operator (from (2.19) and (2.10)) plotted against the number of eigenmodes used in the approximation. The Maury decomposition is more variable between networks, and in particular performs better (i.e. a smaller number of modes are required to capture more of the behaviour) for the image-based airway networks. In particular, the network labelled IB4 only requires a single mode to approximate 80% of the Maury inverse. This mode has a small eigenvalue and hence represents a 'short-circuit' that, given appropriate boundary conditions, may account for the majority of the flow on the network. By contrast, the Frobenius norm of the inverse Laplace operator is largely comprised of high-resistance modes, and so the corresponding eigenvectors are centralized on high-resistance motifs in the network. The mode reconstruction is more accurate in the image-based networks than the Horsfield networks for either operator, which is due to the increased variance of resistance in these networks; however, the effect is stronger for the Maury operator. We now consider fixed boundary conditions for the pressure-drop DP term ¼ e jT j , such that flow is directed from the mouth to the alveolar region. Therefore, any flow heterogeneity that results is due to the distribution of resistance in the networks. To comprehensively cover a range of possible networks, we apply these conditions to the Weibel-type model in §2.4.1 with N = 9 divisions (512 terminal nodes), representative of approximately 1/64th of the conducting airways [1].
As the network asymmetry parameter A is changed from A = 0 (symmetric) to A = 0.98 (all-but-one path having very large resistance), the dominant vectors of the Maury matrix (ranked by |q k |) transition from highest resistance modes (smallest index) to the lowest resistance modes (largest index), see figure 3a. This characterizes the transition from a flow solution that is homogeneously spread over all the   terminal nodes and one that is concentrated on a small collection of terminal nodes that terminate the low-resistance paths in the tree. The A = 0 case is considered in more detail in electronic supplementary material, text S1.4.
For all values of A, there is a small subset of modes in the Maury spectrum that dominate the flow solution (figure 3a). By contrast, the Laplace spectrum does not demonstrate any strong ordering of mode dominance (ranked by |p k |) for the majority of A values and clear separation of modes occurs only at small values of A (figure 3b). For larger A values, the Maury spectrum outperforms the Laplace spectrum in terms of the number of modes required to reconstruct the solutions (figure 3c). However, the difference is small in human lungs, where it is estimated that A ≈ 0.35 [21]. The Maury matrix performs best at extreme values of A, but for all A values a 75% accuracy can be achieved using at most 10 eigenmodes (approx. 2% of the spectrum).

Practical application to simulations of lung ventilation mechanics
In this section, we consider the simple mechanical model of regular tidal breathing introduced in §2.3. In this case, the largest eigenvalue modes of the Maury matrix dominate ventilation as they encode the high resistance paths in the tree (as demonstrated systematically in electronic supplementary material, text S2.1).

Dimensionality reduction in simulations and characterization of realistic airway networks
First, we consider the Horsfield networks H10 and H20 ( §2.4.2), choosing parameters for elastance κ and timescale τ that are representative of an average adult (listed in figure 4). In both cases, all eigenvalues of the Maury matrix satisfied μ k /(κτ) < 0.1. Therefore, the ventilation is approximately homogeneous in these networks and the airway resistance has little effect on ventilation.
To generate networks with high-resistance paths, we systematically applied constrictions to the Horsfield network by selecting a fraction f of the airways (at random) in a given Horsfield order g and reducing their radius by a factor between 50% and 95% (drawn from a uniform random distribution between these limits). The networks generated are idealized examples of the effects of obstructive lung disease, with different networks for disease centred on proximal (g = 15), central (g = 10), or distal (g = 5) airways. This is similar to constrictions applied in models of CF and asthma [33][34][35]. In the constricted networks, we found that as the number of random constrictions is increased (either by increasing f or reducing g), the number of independent large-eigenvalue modes of the Maury matrix rises (figure 4a). This agrees with our interpretation of these modes being centred on high-resistance paths in the tree.
The VH parameter σ V from (2.25) does not depend strongly on g in this airway model, and so its variation in figure 4b is mainly due to the range of f values used (larger f results in more constricted airways and so greater VH). Regardless of f and g values, we see that using only largeeigenvalue modes in (2.25) gives a very good approximation of the simulated VH for all values of f and g. These reconstructions use at most 1% of the full Maury spectrum, and so represent a significant reduction in dimensionality.   [36][37][38].
As κ is reduced, airway resistance has a more significant contribution to the dynamics in (2.20), and so VH is increased due to the heterogeneity of network resistance. The variation between the networks is characteristic of the different network structures, with IB1 and IB2 showing significantly more VH than IB3 and IB4. The Maury spectral approximation (using 1% of the modes) in figure 4c performs marginally worse for smaller values of κ, however, overall the agreement is very strong. The Maury eigenvalues and eigenvectors are unchanged by the parameter sweep in κ and so estimation of the ventilation distribution from equation (2.24) can be repeated for a range of κ values at negligible computational cost. As shown in electronic supplementary material, table S1, the time to compute 1% of the eigenspectrum (as used here) is comparable to the time for a single simulation of the ventilation. Therefore, we were able to predict the relationship between VH and lung elastance accurately and more efficiently using this method.
An explicit example of this reduced dimensionality approximation is given in figure 4d and 4e for the network IB2 (corresponding to the highlighted simulation in 4c), which shows distributions of ventilation ΔV i (see (24)) over terminal units computed by full simulation (figure 4d) and the reduced model (figure 4e). The heterogeneity in ventilation is a result of the airway geometry and topology generated by the space-filling branching algorithm, as well as the branch radii which scales with its parent airway from the CT image, and is captured very efficiently by the reduced model.

Characterization of realistic airway networks
Large-eigenvalue (high-resistance) eigenvectors of the Maury matrix are centralized on those terminal nodes most affected. Therefore, the eigen-decomposition of the Maury matrix is a powerful tool for visualizing and ranking the least ventilated regions, regardless of the chosen value of lung elastance κ. Figure 5 visualizes a subset of the largest eigenvalue eigenvectors of the IB2 network by the size of spheres superimposed onto its terminal nodes. Comparing this to figure 4d, we see that these are localized on the least ventilated terminal nodes. Furthermore, the largest eigenvalue modes highlight the regions of poorest ventilation, as the eigenvalue is a measure of the resistance of the mode. The modes in figure 5 are selected by their value of e T jT jv k , which is the coefficient in the mode sum in (2.24), so that they are all important in this reconstruction (for small enough κτ).

Discussion
In this paper, we showed that in many realistic examples the Maury operator is a more efficient tool than the Laplacian operator for dimensionality reduction of the flow calculation on an airway network. Under fixed-pressure boundary conditions, it is not possible a priori to calculate only the eigenmodes which dominate the flow solution without computing the whole spectrum. However, when extended to simulations of lung ventilation, we found it is only the largest eigenvalue modes that significantly contribute to the extent of VH. These modes are generally the most efficient to compute using methods based on power iteration, as used here. Therefore, this provides a method for efficient prediction and characterization of poorly ventilated regions of the lung based on the geometry and topology of the airway tree. The method has significant advantages over previous approaches of dimensionality reduction in airway models, which involve replacing parts of the tree with symmetric models [39,40] based on 'trumpet' models [41] of the airway geometry. These methods remove the asymmetry of the smaller airways, and sacrifice some of the complexity of the system to reduce dimensionality, which is not the case for the spectral graph theory methods used here.
The main limitation of this approach is the restriction to fixed linear resistance relations. This does not fully reflect the dynamics of ventilation in reality where airway resistance changes dynamically and nonlinearly due to inertial effects [42]. Additionally, airways themselves are compliant [43], and dynamic bronchoconstriction occurs in lung conditions such as asthma changing the network resistance over time [44]. In cases where airway compliance becomes a dominant factor the Maury operator will cease to provide a useful description of this system. Nonetheless, the Maury matrix still has the same interpretation for a single instance of the resistance network, and so an interesting future study would be to study the dynamics of the Maury eigenspectrum in a system where airway resistance varies over time.
The results shown here are also limited to artificial lung airway geometries, using a combination of idealized cast-based models and algorithm-generated image-based networks. Therefore, the simulations and results in this paper do not predict how lungs will behave in specific states of health or disease. The models used do not contain enough information about the small airways or tissue compliance to make patient-specific predictions of VH, which was measured in the original study by 3 He MRI [27]. Spectral methods, such as dynamic mode decomposition, have proven useful in similar contexts to generate low-dimensional characterizations of fluid dynamics simulations for machine learning image analysis applications [45]. In a similar way, the method here produces a low-dimensional representation of the airway tree resistance, which we aim to use in inverse modelling methods. For example, in future work, we plan to fit CT-based models to 3 He ventilation MRI data to develop patient-specific models of lung structure and function (similar to the approach in [34]), and this research will make it more feasible by reducing the computational time required to estimate the ventilation distribution for a candidate airway network over a range of parameters.
A further application of these findings is in the characterization of complex airway tree models, such as those acquired from micro-CT imaging [46]. Micro-CT can only be performed on excised lungs, so it would not be possible to directly relate these to function. However, it will be interesting to characterize airway networks associated with various lung diseases. This work provides an efficient method for this as 'problematic' lung regions can be visualized using the large-eigenvalue eigenmodes of the Maury matrix only.
The methods in this paper are generic to resistance trees, and hence we expect them to have applications outside of lung mechanics. First, other resistance trees in biology, such as the pulmonary vasculature, can be analysed using the same principles. Second, in discrete mathematics, the steady diffusion problem is identical to the linear resistance problem, therefore diffusion on tree networks can also be studied within this framework. For example, transport of gas in the respiratory zone of the lungs is dominated by diffusion on the tree network of acinar ducts, and so in future work we aim to use this formulation to produce lowdimensional representations of models of these networks.
In summary, this paper has demonstrated how spectral graph theory provides a powerful tool for dimensionality reduction in the analysis of lung ventilation. This is achieved by reformulating the physical model for gas transport using a linear operator (the Maury matrix) that captures patterns of ventilation heterogeneity more efficiently than the traditional Laplacian operator, by exploiting the tree-like structure of the network. This approach shows potential for the description of a wider class of transport processes on tree-like networks.
Data accessibility. The source code and executable used for generation of the image-based airway networks are published here https://dx.doi. org/10.5281/zenodo.3709073 with accompanying documentation to be found in the github repository https://github.com/CarlWhitfield/Space_filling_tree. The source code and executable used for the spectral-graph theory computations are published here https:// dx.doi.org/10.5281/zenodo.3708114 with accompanying documentation to be found in the github repository https://github.com/ CarlWhitfield/Spectral_decomposition. The artificial and imagebased airway models used can be accessed in plain-text format and vtk format, as well as the airway centreline data generated from CT data is published here https://dx.doi.org/10.5281/zenodo.3709105.