## Abstract

Both the development of accurate models of lung function and their quantitative validation can be significantly enhanced by the use of functional imaging techniques. The advent of hyperpolarized noble gas magnetic resonance imaging (MRI) technology has increased the amount of local, functional information we can obtain from the lung. In particular, application of ^{3}He to measure apparent diffusion coefficients has enabled some measure of lung microstructure and airspace size within the lung. Models mimicking image acquisition in hyperpolarized gas MRI can improve understanding of the relationship between image findings and lung structure, and can be used to improve the definition of imaging protocols. In this paper, we review the state of the art in hyperpolarized gas MRI modelling. We also present our own results, obtained using a Monte Carlo approach and a realistic alveolar sac geometry, which has previously been applied in functional lung studies. In this way, we demonstrate the potential for models combining lung function and image acquisition, which could provide valuable tools in both basic studies and clinical practice.

## 1. Introduction

Chronic obstructive pulmonary disease (COPD) is an umbrella term encompassing a range of respiratory diseases that affect millions of people around the world. The term includes chronic bronchitis, emphysema, chronic obstructive airways disease, chronic airflow limitation and some cases of chronic asthma. COPD was ranked as the fourth most common cause of death worldwide in 2000 (Pauwels *et al*. 2001), and is predicted to become third by 2020 (Murray & Lopez 1996). Despite the impact COPD has on society, disease assessment and accurate monitoring of its progression has been relatively poor. This is mostly due to the lack of tools for regional and longitudinal screening (Woods *et al*. 2006).

The traditional way to assess the severity of COPD is by monitoring integrated indices of global lung function (Pineda *et al*. 1984). Spirometry, the standard test that uses airflow obstruction as a diagnostic criterion, measures the ratio of forced expiratory volume in 1 second or forced vital capacity. Pulmonary function tests (PFTs) are insensitive to regional changes (Alm *et al*. 2007). Imaging is increasingly being used to assess pulmonary disease, with computed tomography (CT) being the most commonly used method to assess the lung in detail. Chest radiographs are insensitive in detecting changes due to COPD. CT, and more specifically high-resolution CT (HRCT), and multislice CT are more sensitive and specific. Areas of emphysema are identified by decreased attenuation. Several studies have shown an excellent correlation between CT detection of emphysema and histology (Gevenois *et al*. 1995; Bankier *et al*. 1999). However, it is uncertain how valuable the technique is in patients with mild emphysema. There are two main disadvantages of CT in studying patients with lung disease. The first is the presence of ionizing radiation, preventing its use in longitudinal studies. The second is its inability to provide dynamic functional imaging.

Functional lung imaging modalities have been intensively developed over recent years. Functional information can be obtained using nuclear scintigraphy (Howarth *et al*. 1999) and single photon emission CT (Newman *et al*. 2003). However, these methods have very poor spatial resolution and convey very limited information about lung microstructure, in addition to exposing the patient to ionizing radiation. The method that has gained the most attention recently is hyperpolarized gas magnetic resonance imaging (MRI; Fain *et al*. 2007), which has been reported to be more sensitive than CT in COPD evaluation (Ley *et al*. 2004). In comparison with other functional imaging modalities, hyperpolarized gas offers better spatial resolution and signal-to-noise ratio, and is harmless for patients. The future for hyperpolarized gas MRI is in longitudinal studies aimed at monitoring regional COPD progression as well as for use in clinical trials aimed at the development of pharmaceutical agents in which repeated measurements are required (Woods *et al*. 2006; Beckmann *et al*. 2007).

Our primary focus in COPD evaluation is on emphysema, which is characterized by the loss of elasticity of the lung tissue as a result of the destruction of tissue supporting the alveoli. The known result is the enlargement of airspaces as the wall separating the adjacent alveoli degrades. Hyperpolarized gas MRI exploits this effect by measuring gas diffusivity in lung spaces, thus providing a non-invasive method to quantify changes in regional lung microstructure at the alveolar level. This effect was first suggested by Chen *et al*. (1999) and Saam *et al*. (2000) and later derived by Yablonskiy *et al*. (2002). Using different MRI pulse sequences, a range of measurements can be performed including spin density, dynamic ventilation, oxygen-sensitive imaging and diffusion-weighted (DW) imaging (Pineda *et al*. 1984; Eberle *et al*. 2000; Schreiber *et al*. 2000; Kauczor *et al*. 2002; Ley *et al*. 2004). In this paper, we will focus on the latter technique, as it is the most relevant one to assess pulmonary airspaces at the levels of bronchioles and alveoli.

The aim of this study was to apply the state-of-the-art DW-MRI simulation methodology described by Grebenkov (2007) and Sukstanskii & Yablonskiy (2008) on an anatomically realistic alveolar geometry (Burrowes *et al*. 2004). By using a geometry that mimics real anatomy, we expect to be able to link image findings to real anatomical changes in COPD, and to make it possible to produce models combining hyperpolarized gas MRI acquisition and lung function at the microscopic level (which is the original purpose for which the geometrical model by Burrowes *et al*. (2004) was used). Our ultimate goal of these kinds of studies is a better understanding of the relationship between lung function and image measurements.

## 2. State of the art

### (a) Hyperpolarized gas MRI technology

Two agents are currently used in hyperpolarized gas MRI applications: inert noble gas isotopes He-3 and Xe-129. Both gases have nuclear magnetic moments that can be polarized using spin-exchange optical pumping (Bouchiat *et al*. 1960; Honig *et al*. 2000). Although most published studies have been conducted using He-3, driven by its better spatial resolution, there is a growing interest in Xe-129. This is due to both economical reasons (cost of He-3) and the solubility of Xe-129 in tissue and blood, which makes it possible to measure gas diffusion across the alveolar–capillary membrane and into blood plasma and erythrocytes (Månsson *et al*. 2003; Ruppert *et al*. 2004; Abdeen *et al*. 2006).

Noble gas isotopes are hyperpolarized using optical pumping techniques (Happer *et al*. 1984). In contrast to proton imaging, where the small thermal equilibrium nuclear polarization re-establishes itself after the application of a radiofrequency (RF) pulse, the non-equilibrium polarization of the noble gas is irreversibly depleted by the action of the RF pulses and through *T*_{1} relaxation. Hence, hyperpolarized gas MRI is commonly performed using low flip-angle MRI pulse sequences that gradually use up the finite magnetization. Also, dedicated RF amplifier and transmit/receive coils, tuned to the Larmor frequency of the gas, are required.

DW-MRI is traditionally realized using a pulsed-gradient spin-echo sequence (Stejskal & Tanner 1965). However, in hyperpolarized gas MRI, spoiled gradient-recalled echo sequences with a bipolar DW gradient are employed due to short transverse relaxation times and potential difficulties in achieving a uniform 180° RF pulse. The bipolar field gradient can be applied on any of the axes, allowing diffusion to be measured in any arbitrary direction. However, most of the experimental studies measure diffusion only in one direction. The underlying assumption is that, while lung microstructure is anisotropic at the microscopic level, each MRI voxel is so large, compared with these structures, that the amount of diffusion can be represented by a scalar (Schreiber *et al*. 2005) referred to as the apparent diffusion coefficient (ADC).

### (b) Clinical applications

By 2005, over 1000 patients and volunteers had been studied with He-3 imaging (van Beek & Wild 2005), confirming it to be a safe technique for investigating lung disease. Its sensitivity in the detection of early disease has been assessed, compared with induced histological changes in small animal models, and it has been trialled as an investigative tool, albeit in small numbers, in a variety of pathological conditions. The diseases investigated include asthma, cystic fibrosis, emphysema, lung cancer, lung transplantation and pulmonary embolic disease.

In healthy subjects, the distribution of helium was relatively homogeneous within the airspaces (de Lange *et al*. 1999). A comparison of healthy non-smokers and smokers demonstrated an increase in the mean ADC and heterogeneity in the smokers (van Beek & Wild 2005). Several studies (Kauczor *et al*. 1996; Saam *et al*. 2000; Ley *et al*. 2004) have been performed on patient and animal models with varying degrees of severity of emphysema using ventilation distribution and ADC maps. There appears to be relatively good correlation between PFTs, CT, histology and ADC measurements, with hyperpolarized helium being more sensitive for disease detection than CT and PFTs (Salerno *et al*. 2002; Peces-Barba *et al*. 2003; Woods *et al*. 2006; Tanoli *et al*. 2007). Studies on healthy and emphysematous patients have shown good reproducibility of the measurements. The application of hyperpolarized gas on lung cancer, fusing He-3 MRI to CT for functional weighted intensity-modulated radiotherapy, has also been reported in a feasibility study of six patients (Ireland *et al*. 2007).

Lung transplant recipients with graft-versus-host disease, manifest as bronchiolitis obliterans, have been studied by Zaporozhan *et al*. (2004), with He-3 imaging detecting disease earlier than spirometry and HRCT. In a small group of eight children with cystic fibrosis, comparison of He-3 imaging with functional measures of airway obstruction showed good correlation and also that children are able to tolerate helium imaging (Koumellis *et al*. 2005). He-3 MRI has also been used to assess regional patterns of airway closure in asthma (Altes *et al*. 2001; de Lange *et al*. 2006). Response to therapy has been investigated in a small number of patients with asthma (Samee *et al*. 2003).

It has also been shown that different imaging protocols may provide additional information, and a protocol consisting of He-3 ventilation imaging and three-dimensional perfusion MRI enabled the detection of ventilation defects, gas trapping and three-dimensional perfusion (Panth *et al*. 2002).

Owing to its potential to derive regional ventilation–perfusion ratios using oxygen-dependent signal changes, hyperpolarized gas MRI may be of use in evaluating pulmonary vascular disease. It has been suggested that the method may help in planning for thromboendarterectomy and treatment for patients with chronic thromboembolic pulmonary hypertension (van Beek & Wild 2005).

### (c) DW-MRI

DW-MRI is used to measure the physical diffusion of hyperpolarized gas within the lung airspaces. In any given medium, the diffusion of gas molecules is characterized by Brownian motion resulting in root-mean-square displacement (2*D*Δ*t*)^{1/2} along a single axis, where Δ*t* is the diffusion time and *D* is the diffusion coefficient. The free diffusion coefficient (when no barriers are present) for diluted He-3 and Xe in air at atmospheric pressure and 37°C is equal to 0.88 and 0.14 cm^{2} s^{−1}, respectively (Miller *et al*. 2007). In short-regime DW-MRI measurements with a diffusion time of 2 ms, the molecule free displacement is 0.59 mm, making it sensitive to lung microstructure at the alveolar level (alveolar inner ‘diameter’ being 0.125 mm on average; Haefeli-Bleuer & Weibel 1988).

Historically, sensitivity to molecular diffusion in MRI has been quantified using the so-called *b*-values, which describe the degree of diffusion weighting. The *b*-value depends on the diffusion gradient strength and the diffusion time, both of which are experimentally controllable. Equation (2.1) describes the situation for a gradient-echo MRI sequence with bipolar sensitizing gradient pulses with no gradient ramp time(2.1)where *γ* is the gas gyromagnetic ratio (*γ*=2.037894×10^{8} s^{−1} T^{−1} for He); *g* is the maximum strength of the bipolar gradient; *δ* is the gradient duration; and *Δ* is the gradient separation.

In DW imaging experiments performed with low *b*-values (1–4 s cm^{−2}), the phase accumulated by diffusing molecules is small and the phase distribution is well approximated by a Gaussian function (Stepišnik 1999). At these low *b*-values, motion can be described using the ADC(2.2)where *S*_{0} is an image acquired with no diffusion weighting.

Considering the local anisotropy of airway trees, it is reasonable to expect anisotropic diffusion. However, experiments either measuring diffusion in three different directions or calculating a full diffusion tensor from six directions have shown no evidence of anisotropy (Schreiber *et al*. 2005). This may be attributed to the coarse resolution of hyperpolarized He images, with a typical voxel size of 5×5×20 mm, where airway orientations are being averaged. For the currently available image resolution in short time regime experiments, the assumption of independence of ADC values from the direction of measurement seems to be valid.

Owing to the restriction of gas diffusion by the lung microstructure, the signal decay is no longer mono-exponential when higher *b*-values are applied, which indicates that the accumulated spin phase is not Gaussian. Higher *b*-values can be obtained by increasing experimental diffusion time or strength of the sensitizing gradient. In long-range experiments, where diffusive displacements can be measured over time scales as long as *T*_{1}, molecule displacements are in the range of centimetres. Larger displacements are affected by airway connectivity on the acinar level (Bartel *et al*. 2008). Resulting ADC values are smaller than short-range ADCs (approx. 10 times smaller in healthy lungs; Woods *et al*. (2004)).

Imaging at higher and multiple *b*-values reveals a complex non-exponential signal decay in pulmonary airspaces. One approximation that has been used in the analysis of hyperpolarized gas MRI data is sometimes referred to as diffusion kurtosis imaging. The non-Gaussian signal decay can be described using (Jensen *et al*. 2005)(2.3)where ADK is apparent diffusion kurtosis.

ADK is a dimensionless metric of the departure from Gaussian diffusion behaviour (Jensen *et al*. 2005). Some authors have suggested that kurtosis is particularly sensitive to diffusion over longer distances—along longitudinal axes of bronchi and bronchioles (Trampel *et al*. 2006). In Jacob *et al*. (2007), it has been shown that ADK is related to diffusion anisotropy and should be sensitive to changes in tissue organization. However, the same authors also demonstrate that equation (2.3) will describe the signal decay only over a limited range of *b*-values. A different approach is to use models of lung geometry to derive a description for the diffusion attenuation of the DW hyperpolarized gas MRI signal.

### (d) Models of lung geometry

Secondary pulmonary lobules are the smallest units of lung structure separated by connective tissue septa (Weibel & Taylor 1988), with a size of approximately 1–2.5 cm in diameter. Secondary pulmonary lobules are normally built from between 3 and 24 acini (Itoh *et al*. 1993). A single acinus is functionally recognized as the largest airway unit in which all airways participate in gas exchange (Rodriguez *et al*. 1987), with a size ranging from 6 to 10 mm (Osborne *et al*. 1983). The smallest structures taking part in gas exchange are alveoli, with approximately 480 million of them, on average, in the human lung (Ochs *et al*. 2004).

The presence of alveoli is first encountered past the terminal bronchioles distal to the trachea and located approximately at the 14th–15th generation of the bifurcating airway tree (Haefeli-Bleuer & Weibel 1988). Within an acinus, clusters of alveoli increase in number and gradually cover the acinar airway entirely. Intra-acinar airways are divided into three types: respiratory bronchioles (partly covered with alveoli), alveolar ducts (entirely covered with alveoli) and alveolar sacs (blind-ending terminations of alveolar ducts). There are 12 generations recognized within an acinus. The longitudinal path, connecting transitional bronchiole to alveolar sac, has an average length of 8.25 mm (5.3 mm for paths ending in generation 6 and 12 mm for paths ending in generation 12).

Airway generations, respiratory bronchioles, alveolar ducts and sacs decrease in diameter moving towards the terminal air spaces, while alveolar diameters increase (Rodriguez *et al*. 1987; Hammersley & Olson 1992). By comparison, the alveolar duct outer diameter remains constant and equal to approximately 700 μm (Haefeli-Bleuer & Weibel 1988). From the point of view of hyperpolarized gas MRI simulations, the most interesting parts of the lung are single airways covered with alveoli or whole gas exchange units (acini). This is because acini are functionally separated units where the gas moves by diffusion and where most of the inhaled gas resides (95%).

#### (i) Lung models

In order to build computational models, suitable geometric simplifications have been developed in the literature ranging from single branches to full acinar models. In Haefeli-Bleuer & Weibel (1988), airways were depicted as simple branches covered with alveoli (figure 1). This simple model was adopted by Yablonskiy *et al*. (2002), where, for the first time, the relationship between ADC values and airway geometry was derived. Airway dimensions were adopted from Haefeli-Bleuer & Weibel (1988) as follows: the internal airway radius *r* falls from 250 to 135 μm and the outer radius *R* (including the sleeve of the alveoli) remains constant at 350 μm (figure 1*a*). The same authors used a similar model for computer simulation of gas diffusion by Sukstanskii & Yablonskiy (2008), simplifying it to a periodic cylinder structure, where each alveolus covers one-quarter of an annular ring. The cylinder has an infinite length and consists of 350 μm long segments (figure 1*b*). In finite-difference simulations presented by Fichele *et al*. (2004*a*), a similar three-dimensional geometry was used, but this time the alveoli are spherical in shape (figure 1*c*).

Although even simple cylindrical models have been shown to produce ADCs consistent with clinical findings (Sukstanskii & Yablonskiy 2008), they represent a single airway and are far from modelling the realistic lung tissue organization.

Fichele *et al*. (2004*b*) suggested that more realistic models than simple cylinders should be used to describe airway and alveoli geometries. A grape-like structure is proposed, with branching angles resembling those of lung tissue (figure 1*d*). Emphysema was simulated by creating additional pathways between branches. However, difficulties in generating this kind of structure in three dimensions restricted their investigation to two-dimensional models.

Two models of the whole acinar geometry have been described. The Kitaoka model (Kitaoka *et al*. 2000) uses a three-dimensional labyrinth filling a cubic volume to represent the acinus (figure 2*a*). This structure has a total surface area similar to the real acinus, and resembles it also in its dichotomic non-symmetric branched structure. Some degree of emphysema may be simulated by removal of random cell walls that will lead to creation of air-trapping loops (Grebenkov 2007).

The second group of acinar models is focused on representing intra-acinar branching patterns (figure 2*b*). These are ‘tree models’, applying an asymmetrical multiple-branch-point model derived from the anatomical measurements by Haefeli-Bleuer & Weibel (1988). These models have been used for simulations by Verbanck & Paiva (1990, 2007).

Although both described models resemble real lung structure in terms of volume filling and branching, they represent only airways, and thus do not give any insight into the alveolar geometry, which is the main interest of our study. The models are rather limited to the study of emphysema, their main use being in the study of airway obstructions.

More realistic lung tissue models are based on morphological images. In Miller *et al*. (2007), a model of lung tissue was generated using histology images of normal and emphysematous rabbit lungs. In order to define disease progression, surface-to-volume ratios were calculated from each model. This could be easily expended to three-dimensional geometries, but authors restricted their models to only two-dimensional lung tissue representation. Probably the most realistic model to date was presented by Tsuda *et al*. (2008), where finite-element three-dimensional reconstruction of a rat pulmonary acinus, based on high-resolution synchrotron X-ray tomography (resolution 1.4 μm^{3} per voxel), was used (figure 2*c*).

#### (ii) Voronoi lung model

An alternative alveolar sac model was developed by Burrowes *et al*. (2004) by applying a Voronoi meshing technique to represent the alveolar geometry. This model was originally developed to represent the three-dimensional volume-filling structure of alveoli, together with the corresponding capillary network wrapped around the alveolar surface. We will refer to this as the Voronoi model.

From a collection of nodes defined in a volume, a Voronoi diagram partitions the space into small regions in such a way that each ‘cell’ contains the points that are closer to the corresponding node than to any other nodes in the volume. In practice, the Voronoi mesh is implemented using a Delaunay triangulation, as shown in figure 3. By applying this meshing technique, geometries such as the one shown in figure 4 can be generated. More details about the procedure are given in §3.

The geometric model developed by Burrowes *et al*. (2004) and applied in the current study consists of a single alveolar sac (19 alveoli; figure 4*a*) but can, in theory, be extended to represent larger tissue volumes, e.g. an acinar unit. The capillary model developed over this alveolar sac was previously applied to investigate regional variation of capillary perfusion and cellular transit times within a lung (Burrowes *et al*. 2004).

### (e) Simulations of hyperpolarized gas diffusion

The most important phenomena in the hyperpolarized gas MRI experiments are free induction decay (FID) and echo signal formation in restricted geometries. The FID signal may be obtained by solving the Bloch–Torrey equation (Torrey 1956), for which an exact solution was presented by Kenkre *et al*. (1997). However, in the case of restricted diffusion and in the presence of constant field gradient, calculating an exact solution becomes very complicated. The solution to the problem has been proposed by Caprihan *et al*. (1995) and is known as the ‘multiple propagator approach’. It is based on the division of the gradient pulse into successive short time intervals. This solution was reformulated into a matrix form (Callaghan 1997) that allows the calculation of the signal in one-dimensional cylindrical and spherical geometries as shown by Sukstanskii & Yablonskiy (2002). However, to investigate gas diffusion in more complicated geometries, computer simulations are necessary.

#### (i) Physical background

Spin-bearing gas particles within the lung tissue move following a reflected Brownian motion model (Grebenkov 2007). Because the diffusion coefficient for helium is relatively large, the spatial distribution of gas inside the lung can be assumed to be uniform. After the application of a RF pulse, the spin magnetization is flipped into the transversal plane. A magnetic field gradient in a particular direction ** e**, described by its amplitude

*g*and temporal profile

*f*(

*t*), is applied to encode the nuclei positions. A diffusing gas nucleus along a trajectory

**(**

*r**t*) during total time

*T*accumulates phase equal to(2.4)The transverse magnetization at time

*T*can be written in a complex form: e

^{iφ}. The total signal along direction

**is the result of averaging signals from a single nucleus and can be described as the signal expectation over all possible trajectories (Grebenkov 2007)(2.5)If the confining domain is not isotropic, which is the case in lung tissue, the signal depends on the orientation**

*e***of the applied gradient. In theory, the signal could be averaged over all possible gradient directions in space (Grebenkov 2007)(2.6)whereIn practice, if we assume that each image voxel contains a multitude of structures aligned at different orientations, an approximation can be obtained with a small number of gradient orientations, or even with a single one.**

*e*#### (ii) Analytical solution

In Yablonskiy *et al*. (2002), the link between the alveolar geometry and apparent diffusion in the lung was analysed using a cylinder model. Signal attenuation in the lung was approximated by an explicit relationship obtained for restricted diffusion in an ensemble of isolated infinite cylinders randomly placed in space, based on the notion that longitudinal and transverse motions inside alveolar ducts should be significantly different.

By taking the anisotropy of lung tissue on the acinar level into account, equation (2.2) can be written in the following form:where *α* is the angle between the diffusion gradient and the cylinder axis, and *D*_{T} and *D*_{L} are the transverse and longitudinal diffusion coefficients, respectively.

The total signal from a voxel, under the assumption of uniform distribution of airway orientations, is(2.7)where *Φ*(*x*) is the error function and *D*_{AN}*=D*_{L}−*D*_{T} is the diffusion anisotropy.

From experimental results with multiple *b*-values, one can estimate both diffusion coefficients. In the same study, the final relationship between diffusion and alveolar and airway mean radii was proposed. The results derived from the cylinder model were tested using the finite-difference method described by Fichele *et al*. (2004*b*), showing a rather good agreement. In Tanoli *et al*. (2007), the geometry–signal relationship was applied to experimental results successfully. The same authors in Sukstanskii & Yablonskiy (2008) updated the empirical relationship between *D*_{T}, *D*_{L} and geometry using Monte Carlo simulations.

#### (iii) Monte Carlo simulations

The problem of gas diffusion in lung tissue can be solved numerically using the Monte Carlo method, which allows signal decay in complex geometries to be calculated.

The method was first applied to calculate the signal in one-dimensional cylindrical and spherical geometries to validate already existing analytical solutions (Balinov *et al*. 1993; Linse & Soderman 1995). Although the method is well known and has been applied in a wide variety of applications, including the theory of porous media (Mitra *et al*. 1993; Conradi *et al*. 2004), it was only recently that it has been used to explore MRI signal decay in models of lung tissue geometry. In Grebenkov (2007), Monte Carlo simulations were used to calculate signal behaviour in different realizations of the three-dimensional Kitaoka model.

The same approach was applied by Miller *et al*. (2007) to calculate signal attenuation in realistic two-dimensional lung tissue models generated from histology. The focus of that study was on very short diffusion times and the relationship between diffusivity and surface area-to-volume ratio (*S*/*V*) of the pore space.

Monte Carlo simulations by Sukstanskii & Yablonskiy (2008) were used to derive the relationship between longitudinal and transverse diffusivities and geometrical parameters of airways and alveoli in the ‘cylinder’ model, diffusion time *Δ*=1.8 ms and *b*-values ranging from 1 to 10 s cm^{−2}. Simulations confirm non-Gaussian signal behaviour, showing the dependence of diffusion on *b*-values. The results, which appear to be valid for healthy lungs and mild emphysema, predict strong dependence of measured ADC on the experimentally controllable diffusion time, which suggests that when comparing experimental data obtained with different pulse sequences that phenomenon should be taken into account. The proposed diffusion–geometry relationship is claimed to be valid for characterizing lung tissue from gas MRI experiments with multiple *b*-values, and, indeed, gives a very good fit for the data previously presented by Tanoli *et al*. (2007).

#### (iv) Finite-difference approach

An alternative to Monte Carlo for DW-MRI simulation is the finite-difference method (Zientara & Freed 1980). The technique is described by Fichele *et al*. (2004*b*), where a number of assumptions are made in order to apply the Bloch–Torrey equation to the problem. A number of structures were investigated by simulating diffusion in three-dimensional porous models, two-dimensional ‘tree-like’ and ‘grape-like’ structures. The authors' conclusion was that continuous porous media are not suitable to characterize healthy lung tissue; however, results from two-dimensional ‘healthy’ models were in good agreement with *in vivo* measurements and confirmed the analytical cylinder theory presented by Yablonskiy *et al*. (2002).

A novel approach using finite differences was presented by Verbanck & Paiva (2007). In this case, simulations are focused on long-range diffusivity (order of seconds). Using an acinar branching model, the initial condition of a hyperpolarized gas bolus was specified at different locations of the structure. Concentration distribution was calculated from the gas transport equation of diffusion and convection. Although the theory behind this approach is different from other approaches to simulate hyperpolarized gas MRI, values obtained for long-range ADCs are very similar to experimental results. The simulations show an impact of branching pathways on signal attenuation, and lower ADCs as the result.

## 3. Our simulations

### (a) Model geometry

To generate the Voronoi model, initial points are uniformly distributed within the given volume. Specific Voronoi cells are defined as airway passage cells and their faces were removed to create an open duct. The adjacent cells simulate alveoli. In this way, a honeycomb-like structure (figure 4*b*), similar to the one observed in morphometric studies of lung tissue (Haefeli-Bleuer & Weibel 1988), is generated. In our simulations, we use the alveolar sac geometry, which consists of 19 adjoining alveoli made up of polyhedral elements, as shown in figure 4*a*. Average dimensions of the alveoli in the model are compared with known anatomical values in table 1.

The model represents only a single generation of intra-acinar airways, an alveolar sac, which is open on the one end and covered with a terminal alveolus on the other end. The outer diameter of the sac is 738 μm, the inner is 246 μm and the segment length is 836 μm. The values are close to anatomical measurements by Haefeli-Bleuer & Weibel (1988), where the alveolar sac average outer diameter is 657 μm, the inner diameter is 251 μm and the segment length 1012 μm. Because we expect a large influence of the diffusion along the airway main axis when compared with diffusion in the plane perpendicular to this axis (which we refer to as longitudinal and transversal diffusion in the rest of the paper), we treat the alveolar sac as a single duct, and we build the whole acinar path that is known to have an average length of 8.25 mm (Haefeli-Bleuer & Weibel 1988). For simplification, we assume that the periodical duct structures are aligned in a straight line.

In order to mimic an early stage of emphysema with our model, walls are randomly removed between adjacent alveoli (figure 4*c*), mimicking the disease definition and creating a model that visually resembles histology of emphysematous lung tissue (Woods *et al*. 2006). For simulation of more severe emphysema, we completely remove walls separating the adjacent alveoli.

### (b) DW-MRI simulation

We simulate the DW-MRI signal using Monte Carlo simulations, following the approach previously published by Grebenkov (2007) and Sukstanskii & Yablonskiy (2008). A gradient-echo sequence with a bipolar DW gradient was used as described above. Simulations were performed using a random-walk algorithm. Initially, the acinar geometry was established and stored in the form of a three-dimensional binary image. Diffusing gas molecules were then placed at random positions within the simulation volume. The trajectory of a single molecule is modelled by a sequence of *j* random jumps. At each time step (Δ*t*=*T*/*j*), each particle moves in a random direction at a distance *l*_{j}=(6*D*Δ*t*)^{1/2}, where *D*=0.88 cm^{2} s^{−1} is the free diffusion coefficient of He-3 in air. Reflections were not simulated: if a selected move penetrates a wall, the displacement is rejected and stays at its original location. Though this can introduce errors, in accordance with previous studies (Sukstanskii & Yablonskiy 2008), it has been found that the introduction of more complicated reflection laws does not affect the final results in the case of our simulations.

At each step *j* of its random walk through the magnetic field gradient, the particle gains a phase(3.1)where and are, respectively, the position of the particle and the magnetic field gradient at step *j*.

In our simulation, we use typical short-regime experiment parameters with diffusion time equal to *Δ*=1.9 ms. The number of jumps is set to *N*=3600, corresponding to a time step Δ*t* equal to approximately 1 μs and jump distance *l*=23 μm.

The MRI signal *S* is calculated by averaging signals from *M* individual particles,(3.2)We use a gradient with a simplified bipolar temporal profile,(3.3)where *T*=2*Δ*. The *b*-values for the simplified gradient profile can be calculated from the following formula, under the assumption that gradient duration is equal to gradient separation (*δ*=*Δ*):(3.4)The number of gas molecules in our simulations is 10^{6}. We calculate the MRI signal decay for different gradient strength values ranging from 0 to 20 mT m^{−1}, with an increment of 2 mT m^{−1}.

The choice of random generator can significantly affect the results. We use the L'Ecuyer algorithm with Bays–Durham shuffle and added safeguards, as has been implemented by Press *et al*. (2007).

## 4. Results

### (a) Healthy tissue

Simulations were based on the geometry shown in figure 4*a*, containing a single generation of the alveolar sac (19 alveoli) with the dimensions in table 1. Inner and outer diameters of the model are 246 and 738 μm, respectively. Diffusion time *Δ* was set to 1.9 ms. Longitudinal and transversal diffusivity values were calculated corresponding to orientations along the axis of the structure in figure 4*a* and perpendicular to it, respectively. Average diffusivity values were also computed to simulate ADC values obtained in clinical studies. The ratio between longitudinal and transversal diffusivities gives an idea of anisotropy with respect to individual airways.

The original geometry shown in figure 4*a* has a length of only 836 μm, which is too short to consider that a diffusion time of 1.9 ms corresponds to a displacement of 1000 μm. This was confirmed by preliminary simulations on this model showing average diffusivity values clearly below reported ADC values for healthy tissues. An extended model was built by extending the alveolar sac to lengths of 5.3 and 12 mm to represent acinar paths ending in the 6th and 12th generations, respectively. Figure 5 shows calculated MRI signal decay and diffusivities for different *b*-values ranging from 0 to 6.5 s cm^{−2}, simulated by changing the gradient value *g* from 0 to 20 mT m^{−1} and applying equation (3.4).

### (b) Emphysematous tissue

To mimic different stages of emphysema, we randomly removed selected alveolar walls from the initial geometry. Three emphysema models were generated in this way, corresponding to a removal of 10, 36 and 100 per cent of the internal walls, respectively. The results for average signal decay and diffusion changes in these three models, each 12 mm long, are shown in figure 6.

## 5. Discussion

Our work aims to investigate short-range signal decay in healthy and diseased lungs as may be measured via hyperpolarized gas MRI techniques in order to gain an insight into the underlying microstructures present within the lungs. To do this, we applied previously presented Monte Carlo technique (Grebenkov 2007) within a model that realistically represents lung tissue (Burrowes *et al*. 2004). The model can mimic different stages of emphysema. The results are comparable with clinical imaging findings in both healthy and diseased lungs. Though previous studies have been able to simulate clinical ADC measurements, they have used strong simplifications (two-dimensional geometry, far from realistic alveolar shapes, etc.). By using realistic geometries in our hyperpolarized gas MRI simulations, we introduce the possibility of assessing emphysema progression in greater detail and of investigating other aspects of lung function, such as gas exchange. The presented geometry offers the option to replicate an acinus rather than a single airway.

The quantitative validation of the presented model and comparison with previously presented simulation, though desirable, are extremely difficult. We are unable to conduct MRI experiments on a structure exactly resembling our geometry, and the aim of the model is to represent a general human population airway. Our validation, as in most previously presented studies, has focused on showing that the range of ADC values matches with those observed in clinical investigations. The main advantage of our approach is the ability to simulate different COPD stages in an anatomically realistic way, achieving changes in diffusivity reported in experiments.

Figure 5 shows the results of simulations with two model lengths and different *b*-values. A change in the structure length from 5.3 to 12 mm does not seem to influence signal decay significantly, showing only a slight decrease in values for the longer structure (figure 5*a*). The average diffusion calculated at the *b*-value of 1.6 s cm^{−2} increases by approximately 0.005 cm^{2} s^{−1} for the longer model. Average diffusion values for both models are in the range of 0.135–0.16 cm^{2} s^{−1} (figure 5*b*), which is at the lower end of reported experimental values for healthy subjects in the range of 0.15–0.25 cm^{2} s^{−1} (Sukstanskii & Yablonskiy 2008). Both longitudinal and transversal diffusion coefficients (figure 5*c*,*d*) show a dependence on *b*-values, which is attributed to a different degree of non-Gaussian signal behaviour. As expected, longitudinal diffusivity increases in the 12 mm structure compared with the 5 mm one, whereas transversal values remain almost unchanged. The diffusion values for both models for *b*-values ranging from 0 to 6.5 s cm^{−2} are in the range of 0.26–0.33 cm^{2} s^{−1}, and transversal ones are in the range of 0.073–0.075 cm^{2} s^{−1}.

Reported mean isotropic ADC values measured in healthy human are in the range of 0.15–0.25 cm^{2} s^{−1}, with diffusion times ranging from 5.8 to 1.3 ms. The variability is caused by intersubject differences and different MRI experiment settings such as gradient strength and diffusion time (Sukstanskii & Yablonskiy 2008). The direct comparison of simulation results with the literature is difficult. Most of the authors do not specify exactly the experimental conditions, giving not enough parameters to replicate the experiments with certainty. Our simulation of healthy tissue delivers the results that fall rather at the lower end of the reported interval ranging from 0.16 to 0.14 cm^{2} s^{−1} for *b*-values 1–4 s cm^{−2} and diffusion time 1.9 ms. One of the reasons why our calculated diffusion values are low compared with experiments is because our model represents the terminal acinar duct or alveolar sac. In the first few generations of intra-acinar airways, the alveoli are located much more sparsely and are shallower than those in the alveolar sac. To investigate the effect of this, we removed internal alveolar walls completely in the first 3 mm of the 12 mm long airway. The calculated averaged diffusion with *Δ*=1.9 ms for a *b*-value of 1.6 s cm^{−2} was 0.18 cm^{2} s^{−1}, closely matching experiments on healthy tissue (Saam *et al*. 2000). Future work will include introducing this effect in the model, which will require a gradual variation in alveolar depth and density, as the distance from the acinus entry (respiratory bronchiole) increases.

By gradually removing a fraction of alveolar walls, we were able to simulate different stages of emphysema (figure 6). There is a distinct increase in longitudinal diffusion values as a result of the reduction in the alveolar tissue (figure 6*d*). In the extreme case of 100 per cent internal tissue loss, longitudinal diffusion values reach 0.7 cm^{2} s^{−1}, which is not far from free ^{3}He diffusion in the air (0.88 cm^{2} s^{−1}). Figure 6*c* presents changes in transversal diffusivity, showing no clear dependence on the percentage of reduced alveolar walls. Similar behaviour of transversal values has been shown by Sukstanskii & Yablonskiy (2008). The differences between different models are very small and have a negligible effect on average ADC values, and can probably be attributed to numerical errors.

In MRI experiments, the range of ADC values for emphysema varies depending on disease severity. In Mata *et al*. (2007), changes in ADC values measured after inducing mild and moderate emphysema in rabbits were in the ranges of 5–7.8 and 15–18%, respectively. In our 12 mm model and for a *b*-value equal to 1.6 s cm^{−2}, we obtained increases in average diffusion values of 5.8 and 15 per cent when reducing the alveolar tissue by 10 and 36 per cent, respectively, which compares well with these experiments.

For comparison, the simulation results by Sukstanskii & Yablonskiy (2008) conducted on a simple cylindrical geometry with diffusion time *Δ*=1.8 ms produced the following results. Longitudinal diffusivity at *b*=1.6 s cm^{−2} was approximately 0.33 cm^{2} s^{−1} for a healthy geometry and 0.88 cm^{2} s^{−1} for the emphysematous one. Transversal values for *b*=1.6 cm^{2} s^{−1} were approximately 0.09 cm^{2} s^{−1} in the healthy airway model and 0.082 cm^{2} s^{−1} in the emphysematous one. Our simulations at diffusion time *Δ*=1.8 ms and *b*=1.6 s cm^{−2} produced the following values: longitudinal diffusivity equal to 0.31 cm^{2} s^{−1} for healthy lungs and 0.72 cm^{2} s^{−1} for severe emphysema; transversal 0.079 cm^{2} s^{−1} for healthy and 0.08 cm^{2} s^{−1} for the emphysematous one. In our calculation, in contrast to Sukstanskii & Yablonskiy (2008), we used a simplified gradient profile assuming the ramp time to be equal to zero, which may slightly lower the diffusivity results compared with those previously presented. In Grebenkov (2007), where the simulations were performed on a geometry resembling the acinus, with no alveolar structures, the ADC values were as follows. At *Δ*=1.5 ms, with the same gradient profile, ADCs for healthy and severely damaged acinus were approximately 0.22 and 0.55 cm^{2} s^{−1}, respectively, at *b*=1.6 s cm^{−2}. Our model, for *Δ*=1.5 ms, produced ADC values of 0.17 cm^{2} s^{−1} for healthy lung and 0.3 cm^{2} s^{−1} for severe emphysema. The higher values in Grebenkov (2007) can probably be attributed to lack of alveolar walls and simplified airway shapes. This allows the simulation of more advanced emphysema, enabling neighbouring airways to connect. This should also be possible in our model after extending the single airway model to a more complex structure.

Severe emphysema has been characterized by increases in ADC values up to 250 per cent with respect to healthy subjects (Saam *et al*. 2000). These values are close to free ^{3}He diffusion in the air (0.88 cm^{2} s^{−1}), and can only be simulated by destruction of duct walls, which is not possible in our model in its present state. The complete elimination of alveolar internal tissue in our model resulted in a diffusivity increase of 65 per cent (figure 6), which can be considered as the beginning of severe emphysema. In future work, we will extend our model from a single pathway to a whole acinus, which will give us the potential to simulate severe emphysema, including phenomena such as emphysematous bullae formation in extreme cases. Also, a larger model may enable us to simulate diffusion behaviour for extended time scales.

In our simulations, we used 10^{6} independent particles with random starting positions. A larger number was shown to increase only the signal-to-noise ratio, without significantly affecting average values, while 10^{6} particles provided reasonable simulation times. The extension of the model to larger, more complicated, structures is likely to require high-performance computational tools and/or the definition of multiscale models.

As mentioned in §3*b*, following Sukstanskii & Yablonskiy (2008), we do not explicitly include particle reflection in the walls in our Monte Carlo simulations. For the parameters we have used in our simulations, the approximation does not appear to affect our findings, as was also the case in Sukstanskii & Yablonskiy (2008). However, this is not always a valid approximation and can introduce significant errors with increased time-step duration. Calculating exact reflections carries a high computational cost; however, alternative strategies (such as calculating an alternative step vector that does not intersect a surface) would produce a more realistic approximation while maintaining reasonable simulation times. A detailed study of this aspect will be carried out in our future work.

## 6. Conclusions

We have reviewed the fundamental mechanisms and applications of hyperpolarized gas MRI and, in particular, the state of the art in simulation of DW-MRI, by comparing different lung tissue models. We have also presented our own realistic alveolar sac model and applied existing Monte Carlo methodologies, achieving results that are in the reported range of physical experiments. The presented approach is the first step towards the development of models combining lung function and image acquisition.

For more information, please see http://www.comlab.ox.ac.uk/projects/lungimaging.

## Footnotes

One contribution of 15 to a Theme Issue ‘The virtual physiological human: tools and applications II’.

- © 2009 The Royal Society