Supplementary Text: Spectral graph theory eﬃciently characterises 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.


S1.1 Simulation of periodic ventilation
As detailed in section 2.3 of the main text, to simulate ventilation we define the boundary conditions ∆P term (t) = −κV(t) − P pl (t)e |T | (S2) All of the simulations presented in the text use a sinusoidal pleural pressure P pl (t) = P pl0 + P q(t) = q (0) + q (s) sin(2πt/τ ) + q (c) cos(2πt/τ ), (S5) To ensure no net flow over a periodic cycle, q (0) = 0 must be satisfied. It follows from equation (2.2) in the main text (q = diag(r) −1 N T P) that ∆P (0) int = 0 and ∆P (0) term = 0. Therefore, from (S2), we find We define P pl ≡ −κ (V FRC + V T /2) / |T | such that V FRC is the minimum lung volume during tidal breathing (in the limit of zero airway resistance). We also define P (s) pl = −κV T / (2 |T |) so that V T is the tidal volume (in the limit of zero airway resistance). Substituting equations (S3)-(S6) into equations (2.6) and (2.7) (in the main text) and applying the boundary condition (S1) gives We solve this system of equations using SparseLU decomposition solver in the Eigen C++ libraries by representing these equations as a single matrix system of the form Ax = b as follows where the final two rows are the boundary conditions in (S2).

S1.2 Computation of spectra
The spectra of the idealised networks are computed using the Lanczos algorithm which requires only repeated multiplication of a vector by the operator in question. In general the matrix R has |T | 2 non-zero entries and so it is computationally inefficient to perform such multiplications directly. Instead, the multiplication operation y = Rv can be computed efficiently for an arbitrary vector v using the sparse decomposition R = Tdiag(r)T T . Then, performing the sparse multiplication operations w = T T v, and y = Tdiag(r)w completes the operation. The number of operations required for each multiplication roughly scales as |E| N where N is the number of airway generations in the network.
The Horsfield and CT-based networks have ∼ 30, 000 terminal nodes so we compute the full spectrum using a shift-solve method (see Spectra documentation for implementation [10]). To do this efficiently, for the matrix A ∈ {R, L int } we use the following algorithm.
4. Otherwise, increment n = min(n + n 0 + 1, N e ), where N e is the total number of modes to be found, and return to step 1.
8. Otherwise, increment n = n + n 0 and return to step 6.
9. If m < N e , return to step 5.
When considering the Maury matrix R, step 5 in the above algorithm uses the operation y = (xI − R) −1 v for an arbitrary constant x and vector v, which is expensive to solve for if using R directly. However, we can emulate this operation efficiently by solving equations (2.6) and (2.7) in the main text for the case where the terminal edge resistances are modified such that pressure-drop is given by P ext − P 1 e = v. Then, the operation is equivalent to y = −Q ext for this solution.

S1.3 Generation of CT image-based airway models
The image-based models used in this paper are based on the algorithms developed in [5] and [6]. Airway centreline models were generated automatically from the CT data from the cohort in [7] using Pulmonary Toolkit (PTK) [8]. Registration of the lobes was achieved semi-automatically, by correcting the initial guess generated by PTK to coincide with lobe fissures, which were identifiable in all images. The lobe surfaces were outputted and stored in .stl files. The branching algorithm was written in C++ and used the stl reader header library [9]. The algorithm is initialised with the following steps 1. A point density ρ is calculated by dividing the total volume of the lobes by 30,000 (approximate number of acini).

2.
A 3D grid of points is generated with spacing ρ 1/3 , and the position of each point is randomly perturbed by a uniformly distributed random number between −0.01ρ 1/3 and 0.01ρ 1/3 . This perturbation avoids potential overlap with the .stl gridpoints, which can cause problems in determining which points are inside or outside a given surface, as well as adding a small random element to the generated tree.
3. Each point is determined 'inside' or 'outside' of each lobe by casting a line and counting the number of surfaces crossed. Points deemed outside the whole lung are discarded.
4. The terminal nodes of the CT centreline model are identified and assigned to each lobe using the same method.
Starting from the terminal nodes of the CT centreline tree, the space-filling branching algorithm proceeds as follows: 1. The grid points within a lobe are assigned to the nearest open terminal node of the tree.
2. Each point cloud is split by a plane through the mean position of the points in the cloud. The normal of the plane is the cross product of the direction of the terminating branch and the direction from the terminal node to the point cloud mean position. In some cases, these directions may be degenerate, in which case the parent branch of the terminal branch is used instead of the terminal branch direction.
3. New branches are generated along the direction from the terminal node towards the centre of mass of each cloud and terminated 40% of the way. 4. The branching angle of each branch with the parent branch is checked, and if it exceeds z the branch is rotated to have angle π/3.

5.
If the branch terminates outside the lobe surface, the branch is iteratively shortened until this is not longer the case.
6. If the branch is < 1mm in length, or its point cloud has 1 point in it, the branch is terminal and its terminal node is classed as 'closed'. The nearest point in its point cloud is removed from reassignment in the next round.
7. If open terminal nodes remain, return to step 1.
Once the branching is complete, diameters are assigned to the generated airways using equation (2.27) from the main text, so 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 .

S1.4 Spectral properties of the Maury matrix on symmetric dyadic trees and sub-trees
Consider the case where a subset of the terminal nodes T j ⊂ T of the tree are indistinguishable by symmetry, i.e. all paths from their lowest common ancestor (LCA) edge e j are identical in resistance structure. By symmetry, T j can be replaced by a single terminal node connected to v i (where e j = {v i , v i }) by a single edge with resistance equal to the effective resistance of the replaced sub-tree. In fact, approximating distal sub-trees as symmetric is a method of dimensionality reduction in these systems [1,2] which has its origin in so-called "trumpet" models of the airways [3]. However, we find that this is a specific example of the more general method of dimensionality reduction proposed here based on the Maury spectrum.
Airway trees are almost always dyadic (each branching point being a bifurcation), so we explicitly consider that case here. In the dyadic case the number of terminal nodes in the subtree where N + 1 is the number of generations of airways in the sub-tree, i.e. there are N bifurcations between the edge e j and the terminal nodes T j . If we assume the terminal nodes are ordered such that all v i ∈ T \ T j are listed first, and v i ∈ T j are last then R becomes The matrix A ∈ R (|T |−2 N )×(|T |−2 N ) only contains terms corresponding to resistance paths that terminate outside of the sub-tree nodes V j . The off diagonal blocks contain the resistance of paths shared between terminal nodes not in the sub-tree (v i ∈ T \ T j ) and those that are (v i ∈ T j ). These are all the same for all terminal nodes in the symmetric subtree (v i ∈ T j ) as they share the same common ancestor and, by definition, are more distantly related to the other terminal nodes (v i ∈ T \ T j ). Therefore, these blocks can  lines indicate paths to terminal nodes whose corresponding entry in the eigenvector is +1 or −1 respectively (see equation (S16)).
The block Y ∈ R 2 N ×2 N contains resistance paths that terminate in the symmetric tree S j and so has the structure Y = r 0 e 2 N e T 2 N + r 1 where I ∈ R 2 N ×2 N is the identity matrix. The resistances r n are the resistance of any given edge e j in generation n ∈ {0, . . . , N } of the subtree branches E j , as all the airway in a generation are identical by symmetry.
It follows from the definition above that the all-one vector e 2 N is an eigenvector of Y. All other eigenvectors of Y can be given in (non-normalised) form as The indices g and m are the airway generation and index respectively corresponding to the edge e j in the subtree S j as shown in figure S1(b). The vectors in (S16) satisfy e 2 N · y g,m = 0 and so, padded with zeros, are also eigenvectors of the matrix R in (S14) as shown below R 0 s n,m = A0 + (e T 2 N y g,m )r a (r T a 0)e 2 N + Yy g,m = µ g 0 y g,m , (S17) where µ g is the eigenvalue of Y corresponding to the eigenvector y g,m . These eigenvectors of R can also be expressed asv where t j1 and t j2 are the j 1 th and j 2 th columns of the matrix T (that indicate which terminal nodes are descended from the edges e j1 and e j2 , the daughter edges of e j ). Given symmetric boundary conditions on the symmetric sub-tree (∆P i = const for all v i ∈ T j ) these eigenvectors do not contribute to the flow solution at all (see equation (2.17) in the main text), reducing the effective dimensionality of this block from 2 N terminal nodes to 1, and so is equivalent to replacing the symmetric sub-tree with a single edge by Kirchhoff's law (i.e. a trumpet representation).
The above extends trivially to the case where an entire N + 1 generation tree is symmetric, and so the whole tree is fully characterised by a single mode. As was shown in [4], the degeneracy of the 2 n modes associated with each generation n = 0, . . . , N − 1 means there are N + 1 eigenvalues in total which can be labelled µ symm,pf = 2 N R eff ≡ µ 1 , (S19) where R (g) trunc = g g =0 r g /2 g the effective resistance of the tree truncated at generation g and R eff ≡ R (N ) trunc is the effective resistance of the whole tree. The eigenvalue µ symm,pf = µ 1 is the Perron-Frobenius eigenvalue and in the case of symmetric boundary conditions is the only one that contributes to the flux. All others again satisfy e · v k = 0. This also extends to symmetric n-adic trees as shown in the the following section, although we focus on the dyadic case in this paper.
To summarise, the spectrum of the Maury matrix, unlike the Laplacian spectrum, accounts for symmetry reductions exactly. However, these symmetric trees are an idealisation, and so dimensionality reduction by assuming symmetry is restricted to "nearly-symmetric" cases [1]. The spectral graph theory approach proposed here has no such restriction, and can be extended to arbitrarily asymmetric cases, as we show in the results section.

S1.4.1 Perturbation of the symmetric solution
Now consider the case where a single edge e j in generation g p of the otherwise symmetric dyadic tree is perturbed to have resistance r gp + δr. The Maury matrix is modified by the perturbation such that Rv g,m = µ symm,gvg,m + (t T jvg,m )δr t j .
(S22) Therefore, if t T jv g,m = 0 thenv g,m is also an eigenvector of R and its eigenvector remains unchanged. First, this is satisfied for all of the eigenvectors where g ≥ g p as these modes are by definition orthogonal to t j (see (S18)). Second, this condition is also satisfied if the branch e j is not descended from the internal node associated with the modev g,m (in which case, from (S18) all of the entries that are non-zero in t j will be zero inv g,m and vice versa). Linearising the eigenvalue equation Rv k = µ kvk for perturbations about the affected modes gives where we have used the notationv k =v g,m + δv g,m and µ k = µ symm,g + δµ g,m . The modes need to remain orthogonal to the unperturbed modes so the solutions for δv g,m can be expressed as linear superpositions of the other symmetric modesv g ,m . Furthermore, to retain orthonormalityv T g,m δv g,m = 0. Therefore, we can express the eigenvector perturbations to linear order as where m = m(g ) is the affected branch index for generation g and we have used the indices g = −1 and m = 0 to denote the Perron-Frobenius mode (which is not associated with an internal node). Substituting into (S23) and comparing coefficients we can solve for c g ,m to find δv g,m ≈ (t T jvg,m )δr Equations (S25) and (S26) predict that the more distal eigenmodes (i.e. those associated with nodes fewer generations above the perturbation) will be most strongly perturbed in the linear limit (both their eigenvalue and eigenvectors).

S2.1 Eigenvector centralisation in the Maury matrix locates high-resistance paths in airway trees
As shown in detail in section S1.4, the Maury matrix has a number of special properties when part or all of the airway tree is symmetric. Considering fixed boundary conditions (∆P term = e |T | ) then when the entire tree is symmetric only one mode has a non-zero contribution (q k ) to the flow solution, and so the whole tree acts as a single resistor. If the resistance of a single edge in generation g of this symmetric tree is perturbed to break the symmetry of the whole network, the tree can still be divided into g + 1 symmetric sub-trees, which can each be represented as single resistors as shown in [1] and figure S2(a). In figure S2(b)   how the eigenvalues of the g + 1 affected modes changes when the resistance of a single airway is perturbed. The initial trajectories of these curves can be predicted by linear perturbation theory in (S25) and (S26).
In figure S2(b) see that a large increase in resistance of a single airway is mostly captured by the largest eigenvalue mode and in figure S2(c) that the associated eigenvector becomes centralised on the nodes downstream of the affected airway. This corroborates with our interpretation of the largest eigenvalue modes of R indicating high resistance paths in the tree.
Thus large eigenvalue modes have eigenvectors centralised on high-resistance paths in the tree (when such paths exist) and small eigenvalue modes will be centralised on low-resistance paths (when they exist).
However, only a small number of these modes will be important to the flux on the network depending on the boundary conditions. When the external pressure drop ∆P term is perpendicular to a mode (∆P T termvk = 0) that mode will not contribute to the flow on the network.

S2.2 Matrix dimensions and computation times
Table S1 summarises some of the network sizes and simulation times used in the text for context. In the Horsfield networks (H10 and H20) we demonstrate the cases with no constrictions applied (f = 0), and hence the fewest large eigenvalue modes, and the case with the largest number of large eigenvalue modes in figure   4(a) of the main text.  Table S1: Each simulation or computation time is the mean of 3 runs ± standard deviation on a single 2.6GHz processor. The simulations solve the system of equations in (S12). The last two columns show the time taken to compute the modes related to the largest 1% of eigenvalues of the Maury matrix, and to compute the modes above the cut-off value stated (corresponding to the minimum value of κτ used in the paper).