INTERFACE FOCUS

rsfs.royalsocietypublishing.org

Research

CrossMark

click for updates

Cite this article: Hedges KL, Clark AR, Tawhai MH. 2015 Comparison of generic and subject-specific models for simulation of pulmonary perfusion and forced expiration. Interface Focus 5: 20140090. http://dx.doi.org/10.1098/rsfs.2014.0090

One contribution of 11 to a theme issue 'Multiscale modelling in biomechanics: theoretical, computational and translational challenges'.

Subject Areas:

bioengineering, computational biology, biomedical engineering

Keywords:

pulmonary perfusion, acute pulmonary embolism, forced expiration, wave speed limitation

Author for correspondence:

Merryn H. Tawhai

e-mail: m.tawhai@auckland.ac.nz

Comparison of generic and subject-specific models for simulation of pulmonary perfusion and forced

expiration

Kerry L. Hedges, Alys R. Clark and Merryn H. Tawhai

Auckland Bioengineering Institute, University of Auckland, Private Bag 92019, Auckland, New Zealand

The goal of translating multiscale model analysis of pulmonary function into population studies is challenging because of the need to derive a geometric model for each subject. This could be addressed by using a generic model with appropriate customization to subject-specific data. Here, we present a quantitative comparison of simulating two fundamental behaviours of the lung—its haemodynamic response to vascular occlusion, and the forced expiration in 1 s (FEVj) following bronchoconstriction—in subject-specific and generic models. When the subjects are considered as a group, there is no significant difference between predictions of mean pulmonary artery pressure (mPAP), pulmonary vascular resistance or forced expiration; however, significant differences are apparent in the prediction of arterial oxygen, for both baseline and post-occlusion. Despite the apparent consistency of the generic and subject-specific models, a third of subjects had generic model under-prediction of the increase in mPAP following occlusion, and half had the decrease in arterial oxygen over-predicted; two subjects had considerable differences in the percentage reduction of FEV1 following bronchoconstriction. The generic model approach can be useful for physiologically directed studies but is not appropriate for simulating pathophysiological function that is strongly dependent on interaction with lung structure.

1. Introduction

Multiscale models of lung structure and function are becoming more sophisticated, for example now including mechanotransduction across the cellular to whole organ scales [1], coupling of microcirculatory and large vessel fluid dynamics [2,3] and all simulated within subject-specific imaging-based model geometries. Three-dimensional computational fluid dynamics (CFD) is now used routinely to study complex flow phenomena in the central airways [4,5] and pulmonary blood vessels [6], with some emerging model linkages to tissue or cellular-scale function [7]. CFD models are extending to larger domains to capture the impact of subject-specific geometry, and some studies are starting to consider flow within deforming airway or vessel boundaries [8]. The emphasis on subject-specificity and increased structural detail is in contrast to the many generations of models in the respiratory field that have used the simplest airway or vessel structures possible (typically symmetric—or single path— geometry [9], or regular asymmetry [10]) to simulate physiological function, for example the heterogeneous distribution of shear stress during bronchocon-striction [11], or the relationship between forced expiratory flow (FEF) and airway narrowing [12]. Spatially distributed geometric models that directly map to individual subject anatomy are clearly needed to interpret three-dimensional imaging, e.g. to connect the spatial location of imaged 'ventilation defects' to the distribution of airway resistances [13]; but it is not clear that subject-specific geometry is necessary for other applications. One could argue that the inclination to build yet more detail into the new generation of multiscale models is an academic response to a scientific challenge rather than that it is

THE ROYAL SOCIETY

PUBLISHING

© 2015 The Author(s) Published by the Royal Society. All rights reserved.

needed to address physiological questions. On the other hand, experimentalists question the power of modelling studies that only use a handful of structure-based models, and/or that are 'too simplistic' to represent the massively complex structure and function of the lung. The current study presents a quantitative analysis of multiscale structure-based models for pulmonary perfusion and forced expiration, to determine the extent to which subject-specific geometry is required for valid modelling study outcomes in these areas.

2. Multiscale simulation of pulmonary perfusion

Pulmonary blood flow is not distributed uniformly to the lung tissue. Imaging and microsphere studies have suggested that both gravity and asymmetric vascular branching geometry contribute to the characteristic distribution of blood in the lung and its heterogeneity. Human imaging (for example, using external scintillation counters [14], SPECT [15] or MRI [16]) reveals a gradient in the distribution of perfusion, with a tendency for increasing flow in the direction of gravity. Animal studies using destructive techniques [17] suggest that in-plane heterogeneity is at least as significant as the gravitational gradient, leading to the suggestion that the branching geometry influences blood flow distribution to at least the same extent as gravity. It is very difficult—perhaps impossible—to reconcile these different observations without using a multiscale computational model. Microsphere studies cannot be conducted on humans; hence there are species differences in arterial tree asymmetry that influence the characteristic distribution of blood [18], in addition to any effects of posture via the subject being supine during imaging or administration of microspheres.

Capturing the important mechanisms that determine the distribution of pulmonary blood flow in the human lung or in other species requires a model that includes realistic blood vessel asymmetry and spatial distribution of branching, deformation of the lung tissue in response to gravity and posture, and distensible capillaries that respond appropriately to the pressures that surround them. The model must further capture appropriate interactions between the elastic blood vessels and the surrounding parenchymal tissue, and ideally be coupled to a model of ventilation that can respond to the same regional distribution of tissue volumetric changes and stress distribution, for estimation of respiratory gas exchange.

Several computational models have been proposed that include one or more of these features (e.g. [19-21]); however, they are generally limited in scope, as they are designed to investigate function at a single scale, or—for computational expediency—they use simplifying assumptions that restrict them from more general applications. Clark et al. [3] presented a multiscale biophysically representative model that includes each of the components necessary for examining the relative contribution of the important mechanisms in lung perfusion distribution. With respect to its structure, the model couples subject-specific pulmonary arterial and venous geometries [22] to 'ladder-like' models for the distensible and recruitable pulmonary microcirculation [23] via space-filling vascular trees [3]. The models interact with the tissue in which they are embedded, with vessel transmural pressures estimated from soft tissue mechanics models of lung tissue deformation [24]. The model is consistent with

previous imaging studies as well as higher resolution microsphere studies, predicting a gravitationally dependent distribution overlaid with significant heterogeneity. The importance of coupling geometry, flow and tissue mechanics becomes apparent when considering the relative contributions of gravity and structure to perfusion distribution. Gravity acts via the hydrostatic pressure gradient 'weighting' the blood within the vessels, and by deforming the delicate lung tissue by its own weight (causing a gravitationally directed gradient in the density of pulmonary capillaries). In the integrated model of Clark et al. [3], these two gravitational effects make a similar contribution to the gravitational gradient of flow; i.e. the pressure head on the fluid in the pulmonary vessels is as significant as the so-called 'Slinky effect' [16]. Asymmetry of branching dimensions causes heterogeneity in resistances, so heterogeneous frictional losses and variability in iso-gravitational blood flow. The variability is species-dependent (greater variability in animals with higher branching asymmetry [18]). At the microcirculatory level, blood pressure variability is amplified by the balance of arteriole, venule and air pressures across the capillary bed. This is consistent with the 'zonal' model for pulmonary perfusion [25], which describes how the balance of pressures at the capillary level limits the flow of blood. By including multiple spatial scales (conduit vessels and microcirculation) with their distinctive functional behaviours, and by coupling to the key interactions that the pulmonary circulation has with other lung component functions (lung tissue, gas exchange), this model has been able to show consistency between independent experimental studies for perfusion distribution and response to hypoxia [26].

The integrated model can be parametrized to individual patients or animals using subject-specific imaging and haemo-dynamic data. The challenge considered here is whether a single geometric model can sufficiently and appropriately capture the contribution of vascular structure to key metrics of lung function (mean PAP and arterial blood gases), at baseline and following embolic occlusion in acute pulmonary embolism (APE), such that only one generic structure is required (with customization) to simulate multiple subjects.

2.1. Methods to simulate acute pulmonary embolism in subject-specific and generic models

The current study evaluates simulated function in a 'generic' geometric model that is parametrized to imaging data from a group of subjects that were used in a previous modelling study [27]. Summary methods for simulating perfusion, ventilation and steady-state oxygen transfer (with model equations) are provided in appendix A.

Clark et al. [27] describe a multiscale model analysis of haemodynamic and gas exchange response in a cohort of 12 subjects who were diagnosed with APE. A structure-based finite-element model of the pulmonary arteries, veins, airways and embolus distribution was created for each subject using volumetric computed tomography pulmonary angiograms (CTPAs) that were acquired during routine clinically indicated examination for APE. Imaging data were used to define the model geometry to the level of the first sub-segmental branches. Vessels and airways beyond this level were generated using a volume-filling branching algorithm [28]. The unstrained radius of each vessel was defined based on its Strahler order and data from previous morphometric

Table 1. Haemodynamic and gas exchange quantities predicted by subject-specific and generic models under baseline (non-occluded) conditions. P-values are given for paired Student's f-tests. Quantities shown are mPAP, PVR, alveolar oxygen partial pressure (PAO2), arterial partial pressure of carbon dioxide and oxygen returning to the systemic circulation from the lungs (PaCO2 and PaO2), and the difference between arterial and alveolar oxygen partial pressures (P(a-A)O2).

predicted quantity subject-specific generic p-value

mPAP, Pa (mmHg) 1943 + 104 (14.61 + 0.78) 1987 + 69 (14.94 + 0.52) 0.189

PVR, Wood units 0.36 + 0.06 0.38 + 0.04 0.189

PAO2, mmHg 105.46 + 1.98 105.39 + 1.77 0.833

PaCO2, mmHg 37.55 + 0.45 37.23 + 0.23 0.026

PaO2, mmHg 95.98 + 2.23 99.59 + 1.63 4.51 x 10"6

P(a-A)O2, mmHg 9.47 + 1.34 6.08 + 1.11 8.50 x 10"5

Table 2. Haemodynamic and gas exchange quantities predicted by subject-specific and generic models, post-embolic occlusion of arteries. P-values are given for paired Student's t-tests. Quantities shown are mPAP, PVR, PAO2, arterial partial pressure of carbon dioxide and oxygen returning to the systemic circulation from the lungs (PaCO2 and PaO2), and the difference between arterial and alveolar oxygen partial pressures (P(a-A)O2).

predicted quantity subject-specific generic p-value

mPAP, Pa (mmHg) 2662 + 385 (20.02 + 2.90) 2476 + 384 (18.61 + 2.89) 0.169

PVR, Wood units 0.78 + 0.23 0.67 + 0.22 0.169

PAO2, mmHg 105.04 + 2.71 106.33 + 3.15 0.075

PaCO2, mmHg 62.39 + 15.88 56.53 + 16.63 0.427

PaO2, mmHg 62.91 + 20.27 77.64 + 17.16 0.011

P(a-A)O2, mmHg

42.13 + 20.77

28.69 + 16.27

studies [3]. The elastic vessels were free to distend in response to transmural pressure (Ptm), including changes in blood pressure. The three-dimensional distribution of emboli was determined via a semi-automated procedure [29]. The embolus volumes were calculated in a similar manner to the calculation of vessel dimensions by Lee et al. [30]. The emboli (and their size) were then mapped to the arteries in which they were located.

The current study uses the model geometry (excluding emboli) from subject S1 in Clark et al. [27] as a generic model. This subject was representative of the population (a 46-year-old male with body mass index (BMI) of 30 kg m"2). However, the selection of a generic lung shape for simulations is somewhat arbitrary as the lung size is then customized to the individual. The model geometry is customized to the other 11 subjects by scaling the lung volume (including the lengths of airways and vessels), and using patient measurements of airway and vessel diameters down to the sub-segmental level. The patient-specific distributions of emboli, as defined in the original study, were placed within the generic model geometry.

Tissue deformation, perfusion and ventilation were simulated with the lung in the upright posture. Simulation requires setting boundary conditions for the left atrial pressure, cardiac output and minute ventilation. Baseline cardiac output and tidal volume were estimated for each subject using height, weight and age, and assuming a respiratory rate of 12 breaths per minute and heart rate of 65 beats per minute. The metabolic rate for each subject at rest in the upright posture was calculated [31], and translated to an oxygen uptake rate and minute ventilation [32]. Cardiac output was calculated from the oxygen uptake rate [33]. Systemic oxygen demand was

assumed the same pre- and post-occlusion, and ventilation was assumed to remain at baseline conditions post-occlusion. The perfusion simulation predicts the distribution of vessel and capillary distension and capillary recruitment, and the post-occlusion redistribution of flow which has a gravitational preference [2]. Mean pulmonary artery pressure (mPAP) is predicted as a consequence of the distribution of vessel size. Arterial and alveolar partial pressures of respiratory gases are output by the gas transfer model that is described in appendix A(e). The gas transfer depends on the ratio of minute ventilation to cardiac output, and the regional matching of ventilation and perfusion at the acini.

2.2. Results for the generic model

Table 1 lists haemodynamic and gas exchange quantities predicted by the models. The mPAP, pulmonary vascular resistance (PVR) and PAO2 are not significantly different between the generic and subject-specific models; however, PaCO2, PaO2 and P(a-A)O2 are significantly different between the two types of model at baseline. Table 2 lists the same quantities for post-embolus simulations. As in the baseline conditions, mPAP, PVR and PAO2 are not significantly different between the two types of model, whereas PaO2 and P(a-A)O2 are. In contrast to baseline, PaCO2 is not significantly different despite a large difference in the mean values. This reflects the large variability in PaCO2 in the post-embolus models in comparison to negligible variability at baseline.

Figure 1 shows that despite the apparently good prediction of mPAP by the generic model, there are four of 12 subjects for whom the generic model under-predicts the increase in mPAP in APE. In figure 2, the decrease in PaO2

eu w << -

. <3 E

2 4 6 8 10 12

increase in mPAP, generic model (mmHg)

Figure 1. Increase in mPAP from baseline to post-embolus, for the generic model against the corresponding subject-specific simulation. The dashed line is the line of identity.

' ® £ sc ae 2 a

10 20 30 40 50 decrease in PaO2, generic model (%)

Figure 2. Percentage decrease in arterial oxygen partial pressure (PaO2) from baseline to post-embolus, for the generic model against the corresponding subject-specific simulation. The dashed line is the line of identity; the solid line shows a linear fit to the data with R2 = 0.502.

following occlusion is over-predicted by the generic model for 5-6 of the 12 subjects.

previous section, comparing forced expiration in a set of subject-specific models with a single generic model geometry parametrized to the subject-specific data.

3. Multiscale simulation of forced expiration

The forced expiration in 1 s (FEV1) and its ratio to forced vital capacity (FVC) are the most widely used indices of lung function. While its lack of sensitivity and specificity is acknowledged, forced spirometry remains the most convenient and repeatable test of lung function. Lambert et al. [34] presented a computational model for simulation of forced expiration through a symmetric (single path) geometric model of the conducting airway tree (summarized in appendix A(f)). The model uses wave speed theory to limit the maximum flow of air through collapsible airways. An intrinsic limitation of the single-path model is that all airways in a generation have identical geometry and function; hence the impact of heterogeneity in parallel pathways cannot be studied. This can partially be addressed by using regular airway asymmetry, as in Polak & Lutchen [35]; however, this does not provide a link to spatial locations of pathology, and the ability to customize the model to different subjects is limited. Recent imaging [36-38] and modelling [39,40] studies have highlighted the importance of regional differences in lung structure and function in asthma, suggesting that 'clustering' of airway closure is important in developing ventilation defects. It is unknown whether the location and magnitude of this clustering is important in interpreting forced expiration. Here we take a similar approach to the perfusion model of the

3.1. Simulation of forced expiration in explicitly branched trees

The following methods describe implementation of the Lambert et al. [34] model within a branched (not single-path) geometry. In a similar approach to the analysis in the previous section, models were derived for and parametrized to six normal young adult subjects. The models include inter-subject variability in the shape and size of the lungs and airways, the distribution of branches within the lungs, the rate of reduction of diameter with decreasing branch order and the driving pressure for expiration. One of the subject models was then selected to be a 'generic' model, and this was re-parametrized and simulation results compared against the original subject-specific models. Again, the selection of a 'generic' subject from a population defined as 'normal' is arbitrary and the subject selected in this case is a 29-year-old male with normal radiology and pulmonary function test (PFT) results.

Volumetric multi-detector row computed tomography (MDCT) imaging and PFTs from three males and three females were acquired retrospectively from a study of normal healthy volunteers conducted at the University of Iowa Comprehensive Lung Imaging Center. Subjects had no history of lung disease or smoking, FEV1 > 80% predicted, FEV1/FVC > 70%, normal radiological lung appearance, and were all

aged less than 40 years. Imaging consisted of a spiral scan of the full lung in the supine posture with volume at 55% of vital capacity (VC). Images were acquired using a Siemens Sensation 64 MDCT scanner, with scan parameters of 120 kV, 100 mA, a pitch of 1.2, slice width 0.6 mm and slice interval 0.6 mm. A finite-element model of the lung and airway geometry was derived for each subject from the volumetric imaging as described in Tawhai et al. [28] and summarized in appendix A(a). Airway diameters for the lung supine at FRC were initialized using data from image segmentation (for the segmented airways); for branches that were generated using a volume-filling approach, radii were assigned using a fixed Horsfield diameter ratio (RDH). RDH was initially set for each subject to give anatomical deadspace as estimated from their height [41]. To accommodate a change in Ptm associated with change in lung volume or posture, the diameter at a given Ptm was related to the initialized diameter by

a(Ptm)

D(Ptm) = D(Pjna)

a(Pinit)'

where D(Ptm) is airway diameter at Ptm, Pinit is Ptm during supine imaging (estimated using the continuum model in appendix A(b)), D(Pinit) is the initial airway diameter defined from imaging or by the rate of reduction of diameter with order, a(Ptm) is the area ratio at Ptm, and (P^t) is the airway area ratio at P^. The continuum model is used to update the elastic recoil pressure (Pe) and hence the local Ptm, at each airway and acinus location for each lung volume.

To drive the exhalation, the distribution of Pe and the pressure that is exerted by the respiratory muscles (Pm) are combined to generate a driving pressure, Pd. Pm is assumed to decay exponentially with lung volume (V), following:

Pm(V) = Pm(max) exp (fci(V - Vpef)'2), (3.2)

where Pm(max) is the maximum muscle pressure, VPEF is the volume at peak expiratory flow and k1 and k2 are numerical constants. k1 was set proportional to each subject's FEF25-75 and k2 was set to 1.9 for all subjects. Pm(max) and RDH were fitted independently for each subject to minimize the sum of differences between model and data PEF and forced expiratory flow estimated between 25% and 75% of VC (FEF25-75).

The Lambert et al. [34] model is used to calculate flow during expiration to 75% of VC, or for a minimum of 1 s if FEV1 was greater than 75% VC. Starting with the model expanded to total lung capacity, the distribution of Pe is calculated at each airway and acinus. Equation (A 11) is solved for the pressure drop along the airways, starting at the model periphery with Pd as a boundary condition and assuming total (lung) expiratory flow rate of 0.01 ls_1. Solution then progresses through the airways towards the trachea. Sites of flow limitation are identified (defined as U greater than 99% of c), and the flows in a limited airway and all of its subtended airways are held fixed over the time interval for computation (taken here as 0.05 s). The total expiratory flow is incremented in steps of 0.011 s_1 and pressure drop in all non-limited airways recomputed, until the maximal flow for the current geometry is reached. Lung volume is then decreased by the computed expired volume, Pe is recomputed, and the process repeated. This continues until the minimum expiration volume is achieved.

To evaluate the potential of a generic model geometry for simulating forced expiration, one of the male subject models was selected for further analysis. The parametrization of

12 ^ 10 8

^ : FVC

0 2 4 6 8

expired volume (L)

Figure 3. MEFV curve simulated using a subject-specific parametrized model for a representative subject. The horizontal line denotes 'PEF'; the vertical line denotes 'FVC'.

3.5 4.5 5.5

FEVi in subject-specific model

Figure 4. Forced expiratory volume in 1 s (FEVn) from simulation of forced expiration in subject-specific geometric models of five subjects, against FEV1 in models for the same subjects but using a generic geometry. Simulations are for baseline (non-constricted) geometry. The solid line is the line of identity.

Pm(max) and RqH to minimize the error in PEF and FEF25_75 was repeated for the individual subject data and the generic model geometry. Using the new parametrization values, FEV1 was simulated following homogeneous constriction of all airways by 50% (reduction in cross-sectional area by 50%), and constriction of only the airways in the right upper lobe by 50%.

3.2. Results for subject-specific and generic models

Figure 3 illustrates a maximum expiratory flow volume (MEFV) curve simulated for one subject-specific model (the subject selected to act as a 'generic' model). The measured PEF and FVC are indicated on the figure. Note that FVC is not simulated as this potentially includes airway closure which is not included in the model. Simulations for other subjects have comparable MEFV curves. The RMSE for the predicted and measured FEV1 is 1.76%.

A paired Student's t-test gives p = 0.510 for FEV1 from the subject-specific compared with generic model. The RMSE between the models is 3.28%. The correspondence between the generic and subject-specific results for individual subjects is shown in figure 4. For all but one subject, the parametrization lies on the line of identity. In this single male subject, the generic model over-predicts FEV1 by 7%.

Figure 5 indicates the general closeness of the generic model simulations to the subject-specific models, for 50% constriction of the right upper lobe (p = 0.248) and for 50% constriction of all airways (p = 0.238). In both cases, the same female subject has FEV1 under-predicted by the generic model in comparison to the subject-specific result (22% and 18% for all airways constricted and upper lobe constriction, respectively).

Figure 6 considers post-constriction FEV1 as a percentage of baseline FEV1. Constriction of the upper lobe results in a

ад 3

♦ all airways x upper lobe

1 2 3 4 5 6

FEVj in subject-specific model

Figure 5. FEVn from simulation of forced expiration following bronchoconstriction in subject-specific geometric models of five subjects, against FEVn in models for the same subjects but using a generic geometry. Results are shown for constriction of the right upper lobe airways by 50 (50% reduction in cross-sectional area), and constriction of all lung airways by 50%. The solid line is the line of identity.

smaller decline in FEV1 than constriction of all of the airways. Three of the five subjects have consistent percentage reduction in FEV1 in both cases; two of the five have considerable differences. In particular, the female subject for whom the generic model under-predicted FEV1 has a difference of 13% and 15% (between generic and subject-specific) for the all airway and upper lobe constrictions, respectively.

4. Discussion

This study shows that a generic model structure at baseline or with patient-specific imposed arterial occlusions is sufficient to predict an average 'global' haemodynamic response that is consistent with a set of subject-specific models, when the generic model is parametrized using patient data (including embolus location). That is, the mPAP and PVR are not significantly different between the two types of model. By contrast, pulmonary arterial oxygen is not well predicted under either baseline or embolic conditions.

PVR reflects the spatial distribution of unstrained vessel size and vessel distension ( plus capillary recruitment) in response to Ptm. The key factors that affect the baseline Ptm distribution are the gradient of tissue elastic recoil pressure (acting externally on the blood vessels) and the hydrostatic pressure gradient acting on the blood. The model assumes identical tissue material properties for all subjects, and the subject cohort has a reasonably similar 'height' of lung (cranio-caudal dimension, not shown). The elastic recoil pressure distribution is therefore similar across subject models, and a change to different model geometry will have little effect unless the cranio-caudal dimension is substantially different. Similarly, because the elastic properties of the vessels are assumed identical for all subjects, the effect of the hydrostatic pressure gradient acting on the weight of the blood is also expected to be similar between subjects. The unstrained vessel size is parametrized to the subjects' imaged and segmented vessels, with rate of reduction of diameter with order (RDH) in the non-imaged model vessels consistent with morphometric studies. The PVR in subjects with larger central vessels is therefore lower than for those with smaller vessels, and mPAP is higher (results given by Clark et al. [27]). These characteristics are easily represented in a generic model by scaling of the lung volume and assignment of vascular diameters.

Of concern is that four of the 12 subjects considered here have their increase in post-occlusion mPAP considerably under-predicted (up to fivefold) by the generic model. The clinical data that accompanied the original study did not include measurement of mPAP, so there is uncertainty that the subject-specific models are 'correct' in their prediction and the generic models in error; however, the focus of this study is to determine whether the generic model parametrization can reproduce results from subject models that have different branching geometry. The under-prediction of mPAP increase suggests that capillaries were more readily recruited in the generic model during arterial occlusion. Capillaries recruit in a gravitational pattern (greater recruitment in non-dependent tissue) [2]. A difference in the tissue distribution along the gravitational (cranio-caudal) axis could cause a different pattern of recruitment, and this would depend on the location of emboli.

The PaO2 reflects the ratio of minute ventilation to cardiac output, and the regional distribution of ventilation and perfusion. PaO2 is not well predicted by the generic model for either scenario, being over-predicted at baseline and post-occlusion. Minute ventilation and cardiac output are unchanged between the subject-specific and generic simulations, therefore, the elevated PaO2 suggests a difference in regional ventilation -perfusion matching. The spatial matching of ventilation and perfusion via the airways and arterial trees is essentially the same for any of these models (the airways and arteries are 'matched' in their paths); the difference in distributions likely arises from a difference in the distribution of tissue over the lung height. For example, a broader based lung would have a higher proportion of gravitationally dependent tissue that has a higher (on average) compliance. This would influence regional tissue expansion during inspiration, and the ability of capillaries to recruit when large vessels are occluded.

The parametrized generic model for forced expiration is capable of reasonable predictions of FEV1 when the average results over all subjects are considered. The difference between simulation of baseline or post-constriction function is not significantly different between the generic and subject-specific models (as indicated by t-tests); however, in both cases there is one subject that does not match well. This is further highlighted by the discrepancy in the ratio of post-constriction FEV1 to FEV1 at baseline for at least one subject (F1): when simulated within the generic geometry, this subject was far more sensitive to bronchoconstriction. Similar to the perfusion model, forced ventilation depends on the interaction of the lung tissue with airway resistance. If a subject model has relatively high resistance at baseline (e.g. female), then it will potentially be more sensitive to a difference in lung tissue distribution (via lung shape).

Geometric model creation is often a time-consuming process, whether considering several airways in three-dimensional detail, or the entire pulmonary circulation as a network of one-dimensional tubes. Simplifying the procedure of complex model creation would streamline the process of simulating function in large numbers of subjects. Alternatively, demonstrating that a single generic model geometry is sufficiently representative of the normal population would allow modellers to simulate multiple subjects with minimal geometric adaptation of the generic model. Our results here show that parametrization to a single model geometry can be an adequate approach for some limited simulation studies, but they also suggest that in general a generic model is not sufficiently predictive of subject-specific function to be confident in using it for simulating the functions considered here.

4 О О

N—( 70

-A— generic, all airways -e- subject-specific, all airways generic, upper lobe subject-specific, upper lobe

F3 M2 M3

subject

Figure 6. Simulated FEV1 post-constriction as a percentage of baseline FEV1. Results are compared for subject-specific models and a generic model geometry. The uppermost curves are for constriction of the right upper lobe by 50%; the lower curves are for constriction of all airways by 50%. Subjects are labelled by gender (F, female; M, male).

Ethics statement. Imaging, and subsequent use of data, was approved by the University of Iowa Institutional Review Board and Radiation Safety Committees.

Funding statement. This work was supported by National Institutes of Health grant HL14494.

Appendix A

The following sections summarize all of the components of the mathematical models that are used in this study. Full details can be found in the cited literature.

upright lung volume (i.e. including air and tissue). The finite-element model allows the lung to slide within a rigid thoracic wall cavity under gravity loading, while enforcing the lung surface to remain in contact with the cavity surface. Deformation of the lung in any posture (e.g. upright, supine, prone) is simulated by changing the direction of the gravity vector, and enforcing a shape change of the cavity to match imaging data. The deformation solution from the continuum model is used to estimate tissue volume distributions and the pressures acting external to the airways, blood vessels or capillaries. These values are used as input to simulation models for ventilation [43] and perfusion [3] distributions.

(a) Structure-based models

Finite-element models for the lung lobes, airways and blood vessels are derived to fit segmented data for each subject using the methods described in Tawhai et al. [28] and Burrowes et al. [22]. A finite-element model with tri-cubic Hermite interpolation functions is geometry-fitted to the surface data of each lobe. Airway or blood vessel centrelines are converted to a one-dimensional finite-element mesh with a tree-like connectivity, and the image-based radius is associated with each airway/vessel. Additional airways or vessels are generated using a volume-filling algorithm [28] to fill the lung volumes with a 'tree' from the trachea, pulmonary artery or pulmonary veins, to the approximately 32 000 terminal bronchioles. Each acinus in the model is supplied by one terminal bronchiole, and one accompanying artery and vein.

(b) Lung soft tissue mechanics

Lung tissue deformation under the influence of gravity or shape change of the thoracic wall is simulated using a finite-element model [24]. The lung tissue is assumed to be a nonlinearly elastic compressible and homogeneous continuum, with stress -strain relationship as originally described by Fung et al. [42], but without separation of tissue elasticity and surface tension. The model uses the strain energy density function (W), where

W = |exp(fl/? + b/2),

where J1 and J2 are the first and second invariants of the Green -Lagrangian finite strain tensor and j, a and b are constant coefficients. A theoretical zero stress state is assumed at 50% of the

(c) The blood flow model

Detail on the components of the blood flow model can be found in [3,23,44]; summary information is provided here. The relationship between pressure and flow in each pre-capillary pulmonary artery and vein is described by

128 |xL

Q + PbgL cc«^

where DP is the pressure drop along each vessel, m is blood viscosity, L is vessel length, D is vessel diameter, Q is the vessel volumetric flow, p is blood density, g is gravitational acceleration and 0 is the angle of the vessel with the direction of gravity. Equation (A 2) is formulated for each pre-capillary artery or vein, and solved simultaneously with an equation for conservation of mass for each vessel bifurcation.

The compliance of the extra-capillary vessels is modelled using a pressure-area relationship from Krenz & Dawson [45]. In this model, the vessels are assumed to be distensible up to a maximum Ptm of 1.9 kPa (14 mmHg); beyond this limit, the vessels are inextensible [46]. Over the extensible range

D = (D0)(aPtm + 1), (A3)

where D is strained vessel diameter, D0 is unstrained vessel diameter (at Ptm = 0), and a = 1.49 x 10"4 Pa s 1 is vessel compliance [45].

The pulmonary capillary blood is modelled as a thin 'sheet' that is bounded by a compliant endothelium, following the classic sheet-flow model of Fung & Sobin [47] such that

H3dPtm

where A is the alveolar surface area, S is the proportion of A that is capillaries, f is a friction factor, lc is the average path length from arteriole to venule, H is the thickness of the capillary bed (defined by Fung & Sobin [48] over the height of the lung).

(d) The ventilation model

The distribution of ventilation is simulated by coupling one-dimensional models for airway flow in a distributed branching geometry with a deformable elastic tissue unit at the periphery of each terminal airway [43]. The size and location of the acini is initialized by the continuum model described previously. The elastic behaviour of the terminal units (acini) is derived from the lung tissue continuum model, under the assumption of isotropic unit expansion. The volume-dependent elastic recoil pressure (Pe) and compliance (C ) of each acinus are calculated as

= U (3a + b)(X2

3(3a + b)2(12 - 1)2 (3a + b)2(12 - 1)2

where A is the isotropic stretch from the undeformed reference volume V0, and g = 3/4(3a + b)(A2 — 1). Other coefficients are as per the continuum model.

The acini are coupled to the airways using an equation of motion that balances airway resistance with peripheral tissue compliance:

Paw = RawQ Pe

where Paw is the pressure in the terminal conducting airway, Va is the volume of the acinus, Raw is the resistance of the terminal airway and Q is the rate of flow in that airway. The airway resistance includes the contribution of airway bifurcations to frictional energy loss [49]. To simulate ventilation, the air pressure at the 'entry' to the model was assumed zero (relative to atmospheric) and an oscillating muscle pressure was prescribed at the model periphery.

(e) The gas transfer model

The steady-state gas transfer model of Kapitan & Hemple-man [50] is used to estimate the oxygen (O2) partial pressure in each acinus, for time-averaged rates of acinar ventilation and perfusion. This model assumes that the partial pressure of O2 in the pulmonary capillaries (PcO2) is the same as the partial pressure in the adjacent alveolar air (PAO2) by end inspiration, i.e. that the capillary transit time is sufficiently long to allow equilibration, and there is no diffusion limitation [51]. Then,

Vi ■ P1O2 — VA ■ PAO2 = Qk(CcO2 — QO2), (A 8)

where VI is the rate at which air enters the acinus (lmin"1), VA is the rate at which air leaves the acinus (lmin"1), Q is acinus blood flow (lmin"1), k (constant) accounts for the difference between partial pressure in inspired air and air at body temperature and pressure, PjO2 is the inspired air O2 partial pressure (kPa, or mmHg). CjO2 is the O2 content of

inspired air (j = I), alveolar air (j = A), end-capillary blood (j = C), or mixed venous blood (j = v).

Capillary blood O2 partial pressure and content (in ml O2 per 100 ml) are related by

CCO2 = 15 x 1.34 x p(PcO2) + 0.03PCO2, (A 9)

where the first term on the right hand side represents the O2 bound to haemoglobin, the second term represents O2 dissolved in blood plasma, and p(PcO2) is oxygen saturation. p(PcO2) is a function of PcO2 [52] following:

LKTsPCO2(1+KTsPCO2)3+KRsPCO2(1+KRsPCO2)3

p(PcO2)=

L(1+KxoPcÜ2)4+(1+KRa-PcO2)4

(A 10)

where Kx = 10 x 103 l mol"1,

KR = 3.6 x 106 l mol"1, L =

171.2 x 106 and s (O2 solubility) = 1.43 x 10"4 mol l"1 kPa"

(f) Airway flow limitation

The following section summarizes the Lambert et al. [34] model for forced expiration, which is based on wave speed theory [53]; this theory limits the maximum rate of flow through collapsible tubes to the speed at which disturbances in the flow can propagate upstream (i.e. the 'wave speed'). Following Lambert et al. [34], the pressure gradient along an airway (dP/dx) is defined as

dP _ —f (x) dX - (1 — UVc^

(A 11)

where f(x) represents the viscous pressure drop along the airway, U is local velocity and c is the wave speed. Here f(x) can be described by an empirical equation [54]:

f (x) = PV (« + bRe),

(A 12)

where m is fluid viscosity, V is flow through the airway, Re is the Reynolds number for the airway, and a and b are constants (1.5 and 0.0035, respectively, based on measurements in human airway casts [54]). c2 is defined as

c2 - ridAW (A13)

where A is airway cross-sectional area and p is gas density. dA/dP is the airway pressure-area relationship given by

for Pb

1 - (1 - «0)1 -

for Ptm > 0,

(A 14)

here a is the ratio of airway cross-sectional area to maximal cross-sectional area, a0 is a at zero transmural pressure (Ptm), n1 and n2 are dimensionless coefficients, and P1 and P2 are equations for the vertical asymptotes of the curves given by

«0«!

-n2(1 - a) « 0

(A 15)

(A 16)

c ss :

References

1. Lauzon A, Bates J, Donovan G, Tawhai M, Sneyd J, Sanderson M. 2012 A multi-scale approach to airway hyperresponsiveness: from molecule to organ. Front. Physiol. 3, 191. (doi:10.3389/fphys. 2012.00191)

2. Burrowes KS, Clark AR, Tawhai MH. 2011 Blood flow redistribution and ventilation-perfusion mismatch during embolic pulmonary occlusion. Pulmon. Circul. 1, 365-376. (doi:10.4103/2045-8932.87302)

3. Clark AR, Tawhai MH, Burrowes KS. 2011 The interdependent contributions of gravitational and structural features to perfusion distribution in a multi-scale model of the pulmonary circulation. J. Appl. Physiol. 110, 943-945. (doi:10.1152/ japplphysiol.00775.2010)

4. Yin Y, Choi J, Hoffman EA, Tawhai MH, Lin CL. 2013 A multiscale MDCT image-based breathing lung model with time-varying regional ventilation.

J. Comput. Phys. 244, 168-192. (doi:10.1016/j.jcp. 2012.12.007)

5. Katz I et al. 2014 Using helium-oxygen to improve regional deposition of inhaled particles: mechanical principles. J. Aerosol Med. Pulmonary Drug Deliv. 27, 71-80. (doi:10.1089/jamp.2013.1072)

6. Sohrabi S, Zheng J, Finol EA, Liu Y. 2014 Numerical simulation of particle transport and deposition in the pulmonary vasculature. J. Biomech. Eng. 136, 121010. (doi:10.1115/1.4028800)

7. Madas BG, Balashazy I, Farkas A, Szoke I. 2011 Cellular burdens and biological effects on tissue level caused by inhaled radon progenies. Radiat. Prot. Dosimetry 143, 253-257. (doi:10.1093/rpd/ ncq522)

8. Tang B, Pickard S, Chan F, Tsao P, Taylor C, Feinstein J. 2012 Wall shear stress is decreased in the pulmonary arteries of patients with pulmonary arterial hypertension: an image-based, computational fluid dynamics study. Pulmonary Circul. 2, 470 - 476. (doi:10.4103/2045-8932.105035)

9. Weibel ER. 1963 Morphometry of the human lung. Berlin, Germany: Springer.

10. Horsfield K, Dart G, Olson DE, Filley GF, Cumming G. 1971 Models of the human bronchial tree. J. Appl. Physiol. 31, 207-217.

11. Gillis HL, Lutchen KR. 1999 How heterogeneous bronchoconstriction affects ventilation distribution in human lungs: a morphometric model. Ann. Biomed. Eng. 27, 14-22. (doi:10.1114/1.161)

12. Lambert R, Beck K. 2004 Airway area distribution from the forced expiration maneuver. J. Appl. Physiol. 97, 570-578. (doi:10.1152/japplphysiol. 00912.2003)

13. Tgavalekos N, Tawhai MH, Harris RS, Venegas J, Lutchen KR. 2005 Identifying airways responsible for heterogeneous ventilation and mechanical dysfunction in asthma: an image-functional modeling approach. J. Appl. Physiol. 99, 23882397. (doi:10.1152/japplphysiol.00391.2005)

14. West JB, Dollery CT. 1960 Distribution of blood flow and ventilation-perfusion ratio in the lung,

measured with radioactive carbon dioxide. J. Appl. Physiol. 15, 405-410.

15. Petersson J et al. 2007 Posture primarily affects lung tissue distribution with minor effect on blood flow and ventilation. Respir. Physiol. Neurobiol. 156, 293-303. (doi:10.1016/j.resp.2006.11.001)

16. Hopkins SR, Henderson AC, Levin DL, Yamada K, Arai T, Buxton RB, Prisk GK. 2007 Vertical gradients in regional lung density and perfusion in the supine human lung: the Slinky effect. J. Appl. Physiol. 103, 240 - 248. (doi:10.1152/japplphysiol. 01289.2006)

17. Glenny RW, Lamm WJE, Albert RK, Robertson HT. 1991 Gravity is a minor determinant of pulmonary blood flow distribution. J. Appl. Physiol. 71, 620-629.

18. Burrowes KS, Hoffman EA, Tawhai MH. 2009 Species-specific pulmonary arterial asymmetry determines species differences in regional pulmonary perfusion. Ann. Biomed. Eng. 37, 2497 - 2509. (doi:10.1007/s10439-009-9802-2)

19. Parker JC, Cave CB, Ardell JL, Hamm CR, Williams SG. 1997 Vascular tree structure affects lung blood flow heterogeneity simulated in three dimensions. J. Appl. Physiol. 83, 1370-1382.

20. Bshouty Z, Younes M. 1990 Distensibility and pressure-flow relationship of the pulmonary circulation. I. Single-vessel model. J. Appl. Physiol. 68, 1501-1513.

21. Qureshi MU, Vaughan GD, Sainsbury C, Johnson M, Peskin CS, Olufsen MS, Hill NA. 2014 Numerical simulation of blood flow and pressure drop in the pulmonary arterial and venous circulation. Biomech. Model. Mechanobiol. 13,1137-1154. (doi:10.1007/ s10237-014-0563-y)

22. Burrowes KS, Hunter PJ, Tawhai MH. 2005 Anatomically-based finite element models of the human pulmonary arterial and venous trees including supernumerary vessels. J. Appl. Physiol. 99, 731-738. (doi:10.1152/japplphysiol.01033. 2004)

23. Clark AR, Burrowes KS, Tawhai MH. 2010 Contribution of serial and parallel micro-perfusion to spatial variability in pulmonary inter- and intra-acinar blood flow. J. Appl. Physiol.

108, 1116-1126. (doi:10.1152/japplphysiol.01177. 2009)

24. Tawhai M, Nash N, Lin C, Hoffman E. 2009 Supine and prone differences in regional lung density and pleural pressure gradients in the human lung with constant shape. J. Appl. Physiol. 107, 912-920. (doi:10.1152/japplphysiol.00324.2009)

25. West JB, Dollery CT, Naimark A. 1964 Distribution of blood flow in isolated lung; relation to vascular and alveolar pressures. J. Appl. Physiol. 19, 713-724.

26. Burrowes KS, Clark AR, Wilsher ML, Milne DG, Tawhai MH. 2014 Hypoxic pulmonary vasoconstriction as a contributor to response in acute pulmonary embolism. Ann. Biomed. Eng. 42, 1631-1643. (doi:10.1007/s10439-014-1011-y)

27. Clark AR, Milne D, Wilsher M, Burrowes K, Bajaj M, Tawhai M. 2013 Lack of functional information explains the poor performance of 'clot load scores' at predicting outcome in pulmonary embolism. Respirat. Physiol. Neurobiol.

28. Tawhai MH, Hunter PJ, Tschirren J, Reinhardt JM, McLennan G, Hoffman EA. 2004 CT-based geometry analysis and finite element models of the human and ovine bronchial tree. J. Appl. Physiol. 97, 2310-2321. (doi:10.1152/japplphysiol.00520. 2004)

29. Tawhai MH, Clark AR, Wilsher ML, Milne DG, Subramaniam K, Burrowes KS. 2012 Spatial redistribution of perfusion and gas exchange in patient specific models of pulmonary embolism. In IEEE International Symposium on Biomedical Imaging, Barcelona, Spain, p. 1365-1368. Piscataway, NJ: IEEE.

30. Lee Y-C, Clark A, Fuld M, Haynes S, Divekar A, Hoffman E, Tawhai MH. 2013 MDCT-based quantification of porcine pulmonary arterial morphometry and self-similarity of arterial branching geometry. J. Appl. Physiol. 114, 1191-1201. (doi:10.1152/japplphysiol.00868.2012)

31. Levine JA, Schleusner SJ, Jensen MD. 2000 Energy expenditure of nonexercise activity. Am. J. Clin. Nutr. 72, 1451-1454.

32. Weir JBdV. 1949 New methods for calculating metabolic rate with special reference to protein metabolism. J. Physiol. 109, 1-9. (doi:10.1113/ jphysiol.1949.sp004363)

33. Stringer WW, Hansen JE, Wasserman K. 1996 Cardiac output estimated noninvasively from oxygen uptake during exercise. J. Appl. Physiol. 82, 908-912.

34. Lambert RK, Wilson TA, Hyatt RE, Rodarte JR. 1982 A computational model for expiratory flow. J. Appl. Physiol. 52, 44-56.

35. Polak A, Lutchen K. 2003 Computational model for forced expiration from asymmetrical normal lungs. Ann. Biomed. Eng. 31, 891-907. (doi:10.1114/1. 1588651)

36. Tzeng Y, Lutchen K, Albert M. 2008 The difference in ventilation heterogeneity between asthmatic and healthy subjects quantified by using hyperpolarized 3 He MRI. J. Appl. Physiol. 106, 813-822. (doi:10. 1152/japplphysiol.01133.2007)

37. de Lange EE, Altes TA, Patrie JT, Gaare JD, Knake JJ, Mugler III JP, Platts-Mills TA. 2006 Evaluation of asthma with hyperpolarized helium-3 MRI: correlation with clinical severity and spirometry. Chest 130, 1055-1062. (doi:10.1378/chest.130.4. 1055)

38. Costella S, Kirby M, Maksym G, McCormack D, Paterson N, Parraga G. 2012 Regional pulmonary response to a methacholine challenge using hyperpolarized 3He magnetic resonance imaging. Respirology 17, 1237-1246. (doi:10.1111/j.1440-1843.2012.02250.x)

39. Anafi R, Beck K, Wilson T. 2003 Impedence, gas mixing, and bimodal ventilation in constricted lungs. J. Appl. Physiol. 94, 1003-1011.

40. Venegas JG et al. 2005 Self-organized patchiness in 45. asthma as a prelude to catastrophic shifts. Nature

434, 777-782. (doi:10.1038/nature03490)

41. Hart M, Orzalesi M, Cook CD. 1963 Relation between anatomic respiratory dead space and body size and 46. lung volume. J. Appl. Physiol. 18, 519-522.

42. Fung YC. 1990 Biomechanics: motion, flow, stress, and growth. New York, NY: Springer.

43. Swan AJ, Clark AR, Tawhai MH. 2012 A 47. computational model of the topographic

distribution of ventilation in healthy human lungs. 48. J. Theoret. Biol. 300, 222 - 231. (doi:10.1016/j.jtbi. 2012.01.042)

44. Burrowes KS, Tawhai MH. 2010 Coupling of lung 49. tissue tethering force to fluid dynamics in the pulmonary circulation. Int. J. Numer. Methods

Biomed. Eng. 26, 862-875.

Krenz GS, Dawson CA. 2003 Flow and pressure 50. distributions in vascular networks consisting of distensible vessels. Am. J. Physiol. Heart Circul. Physiol. 284, H2192 — H2203.

Zhuang FY, Fung YC, Yen RT. 1983 Analysis of 51. blood flow in cat's lung with detailed anatomical and elasticity data. J. Appl. Physiol. 55, 1341-1348. 52.

Fung YC, Sobin SS. 1969 Theory of sheet flow in lung alveoli. J. Appl. Physiol. 26, 472 - 488. Fung YC, Sobin SS. 1972 Elasticity of the pulmonary alveolar sheet. Circul. Res. 30, 451-469. (doi:10. 53. 1161/01.RES.30.4.451)

Pedley TJ, Schroter RC, Sudlow MF. 1970 Energy losses and pressure drop in models of human 54.

airways. Respir. Physiol. 9, 371-386. (doi:10.1016/ 0034-5687(70)90093-9)

Kapitan K, Hempleman S. 1986 Computer simulation of mammalian gas-exchange. Comput. Biol. Med. 16, 91-101. (doi:10.1016/0010-4825(86)90034-X)

West JB. 2000 Respiratory physiology: the essentials, 6th edn. Philadelphia, PA: Lippincott Williams and Wilkins.

Monod J, Wyman J, Changeaux J. 1965 On the nature of allosteric transitions: a plausible model. J. Mol. Biol. 12, 88-118. (doi:10.1016/S0022-2836(65)80285-6)

Dawson SV, Elliot EA. 1977 Wave-speed limitation on expiratory flow—a unifying concept. J. Appl. Physiol. 43, 498 - 515.

Reynolds D, Lee J. 1981 Steady pressure-flow relationship of a model of the canine bronchial tree. J. Appl. Physiol. 51, 1072-1079.