## Abstract

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 (FEV_{1}) 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 FEV_{1} 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 bronchoconstriction [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 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 haemodynamic 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 studies [3]. The elastic vessels were free to distend in response to transmural pressure (*P*_{tm}), 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 P_{A}O_{2} are not significantly different between the generic and subject-specific models; however, P_{a}CO_{2}, P_{a}O_{2} and P(a-A)O_{2} 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 P_{A}O_{2} are not significantly different between the two types of model, whereas P_{a}O_{2} and P(a-A)O_{2} are. In contrast to baseline, P_{a}CO_{2} is not significantly different despite a large difference in the mean values. This reflects the large variability in P_{a}CO_{2} 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 P_{a}O_{2} following occlusion is over-predicted by the generic model for 5–6 of the 12 subjects.

## 3. Multiscale simulation of forced expiration

The forced expiration in 1 s (FEV_{1}) 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 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.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, FEV_{1} > 80% predicted, FEV_{1}/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 (R_{D}H). R_{D}H was initially set for each subject to give anatomical deadspace as estimated from their height [41]. To accommodate a change in *P*_{tm} associated with change in lung volume or posture, the diameter at a given *P*_{tm} was related to the initialized diameter by
3.1
where *D*(*P*_{tm}) is airway diameter at *P*_{tm}, *P*_{init} is *P*_{tm} during supine imaging (estimated using the continuum model in appendix A(b)), *D*(*P*_{init}) is the initial airway diameter defined from imaging or by the rate of reduction of diameter with order, *α*(*P*_{tm}) is the area ratio at *P*_{tm}, and (*P*_{init}) is the airway area ratio at *P*_{init}*.* The continuum model is used to update the elastic recoil pressure (*P*_{e}) and hence the local *P*_{tm}, at each airway and acinus location for each lung volume.

To drive the exhalation, the distribution of *P*_{e} and the pressure that is exerted by the respiratory muscles (*P*_{m}) are combined to generate a driving pressure, *P*_{d}. *P*_{m} is assumed to decay exponentially with lung volume (*V*), following:
3.2
where *P*_{m(max)} is the maximum muscle pressure, *V*_{PEF} is the volume at peak expiratory flow and *k*_{1} and *k*_{2} are numerical constants. *k*_{1} was set proportional to each subject's FEF_{25–75} and *k*_{2} was set to 1.9 for all subjects. *P*_{m(max)} and *R*_{D}*H* 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 (FEF_{25–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 FEV_{1} was greater than 75% VC. Starting with the model expanded to total lung capacity, the distribution of *P*_{e} 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 *P*_{d} as a boundary condition and assuming total (lung) expiratory flow rate of 0.01 l s^{−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.01 l 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, *P*_{e} 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 *P*_{m(max)} and *R*_{D}*H* to minimize the error in PEF and FEF_{25–75} was repeated for the individual subject data and the generic model geometry. Using the new parametrization values, FEV_{1} 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 FEV_{1} is 1.76%.

A paired Student's *t*-test gives *p* = 0.510 for FEV_{1} 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 FEV_{1} 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 FEV_{1} 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 FEV_{1} as a percentage of baseline FEV_{1}. Constriction of the upper lobe results in a smaller decline in FEV_{1} than constriction of all of the airways. Three of the five subjects have consistent percentage reduction in FEV_{1} in both cases; two of the five have considerable differences. In particular, the female subject for whom the generic model under-predicted FEV_{1} 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 *P*_{tm}. The key factors that affect the baseline *P*_{tm} 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 (R_{D}H) 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 P_{a}O_{2} reflects the ratio of minute ventilation to cardiac output, and the regional distribution of ventilation and perfusion. P_{a}O_{2} 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 P_{a}O_{2} 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 FEV_{1} 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 FEV_{1} to FEV_{1} 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.

## 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.

**(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
A 1
where *J*_{1} and *J*_{2} are the first and second invariants of the Green–Lagrangian finite strain tensor and *ξ*, *a* and *b* are constant coefficients. A theoretical zero stress state is assumed at 50% of the 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.

**(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 A 2

where Δ*P* is the pressure drop along each vessel, *μ* is blood viscosity, *L* is vessel length, *D* is vessel diameter, *Q* is the vessel volumetric flow, *ρ*_{b} is blood density, *g* is gravitational acceleration and *Θ* 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 *P*_{tm} of 1.9 kPa (14 mmHg); beyond this limit, the vessels are inextensible [46]. Over the extensible range
A 3
where *D* is strained vessel diameter, *D*_{0} is unstrained vessel diameter (at *P*_{tm} = 0), and *α* = 1.49 × 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
A 4
where *A* is the alveolar surface area, *S* is the proportion of *A* that is capillaries, *f* is a friction factor, *l*_{c} 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 (*P*_{e}) and compliance (*C*) of each acinus are calculated as
A 5
A 6
where *λ* is the isotropic stretch from the undeformed reference volume *V*_{0}, and . 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:
A 7
where *P*_{aw} is the pressure in the terminal conducting airway, *V*_{a} is the volume of the acinus, *R*_{aw} 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 & Hempleman [50] is used to estimate the oxygen (O_{2}) partial pressure in each acinus, for time-averaged rates of acinar ventilation and perfusion. This model assumes that the partial pressure of O_{2} in the pulmonary capillaries (*P*_{c}O_{2}) is the same as the partial pressure in the adjacent alveolar air (P_{A}O_{2}) by end inspiration, i.e. that the capillary transit time is sufficiently long to allow equilibration, and there is no diffusion limitation [51]. Then,
A 8
where *V*_{I} is the rate at which air enters the acinus (l min^{−1}), *V*_{A} is the rate at which air leaves the acinus (l min^{−1}), *Q* is acinus blood flow (l min^{−1}), *k* (constant) accounts for the difference between partial pressure in inspired air and air at body temperature and pressure, *P _{j}*O

_{2}is the inspired air O

_{2}partial pressure (kPa, or mmHg).

*C*

_{j}O_{2}is the O

_{2}content of inspired air (

*j*=

*I*), alveolar air (

*j*=

*A*), end-capillary blood (

*j*=

*C*), or mixed venous blood (

*j*=

*v*).

Capillary blood O_{2} partial pressure and content (in ml O_{2} per 100 ml) are related by
A 9
where the first term on the right hand side represents the O_{2} bound to haemoglobin, the second term represents O_{2} dissolved in blood plasma, and *ρ*(*P*_{C}O_{2}) is oxygen saturation. *ρ*(*P*_{C}O_{2}) is a function of *P*_{C}O_{2} [52] following:
A 10
where *K*_{T} = 10 × 10^{3} l mol^{−1}, *K*_{R} = 3.6 × 10^{6} l mol^{−1}, *L* = 171.2 × 10^{6} and *σ* (*O*_{2} solubility) = 1.43 × 10^{−4} mol l^{−1} kPa^{−1}.

**(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 (d*P*/d*x*) is defined as
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]:
A 12
where *µ* is fluid viscosity, 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]). *c*^{2} is defined as
A 13
where *A* is airway cross-sectional area and *ρ* is gas density. d*A*/d*P* is the airway pressure–area relationship given by
A 14
here *α* is the ratio of airway cross-sectional area to maximal cross-sectional area, *α*_{0} is *α* at zero transmural pressure (*P*_{tm}), *n*_{1} and *n*_{2} are dimensionless coefficients, and *P*_{1} and *P*_{2} are equations for the vertical asymptotes of the curves given by
A 15
and
A 16

## Footnotes

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

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