Assessing the relationship between movement and airflow in the upper airway using computational fluid dynamics with motion determined from magnetic resonance imaging

Clinical Biomechanics xxx (xxxx) xxx–xxx

Contents lists available at ScienceDirect

Clinical Biomechanics journal homepage: www.elsevier.com/locate/clinbiomech

Research paper

Assessing the relationship between movement and airflow in the upper airway using computational fluid dynamics with motion determined from magnetic resonance imaging Alister J. Batesa,b,c,⁎, Andreas Schuhd, Gabriel Amine-Eddinee, Keith McConnella, Wolfgang Loewb,f, Robert J. Fleckf, Jason C. Woodsa,b,f, Charles L. Dumoulinb,f, Raouf S. Amina a

Division of Pulmonary Medicine, Cincinnati Children's Hospital Medical Center, Cincinnati, OH, USA Imaging Research Center, Cincinnati Children's Hospital Medical Center, Cincinnati, OH, USA Department of Bioengineering, Imperial College London, UK d Department of Computing, Imperial College London, UK e Siemens Industry Software Computational Dynamics Ltd, London, UK f Department of Radiology, Cincinnati Children's Hospital Medical Center, Cincinnati, OH, USA b c

A R T I C L E I N F O

A B S T R A C T

Keywords: MRI CFD Movement Upper airways Sleep apnoea Image registration

Background: Computational fluid dynamics simulations of respiratory airflow in the upper airway reveal clinically relevant information, including sites of local resistance, inhaled particle deposition, and the effect of pathological constrictions. Unlike previous simulations, which have been performed on rigid anatomical models from static medical imaging, this work utilises ciné imaging during respiration to create dynamic models and more closely represent airway physiology. Methods: Airway movement maps were obtained from non-rigid image registration of fast-cine MRI and applied to high-spatial-resolution airway surface models. Breathing flowrates were recorded simultaneously with imaging. These data formed the boundary conditions for large eddy simulation computations of the airflow from exterior mask to bronchi. Simulations with rigid geometries were performed to demonstrate the resulting airflow differences between airflow simulations in rigid and dynamic airways. Findings: In the analysed rapid breathing manoeuvre, incorporating airway movement significantly changed the findings of the CFD simulations. Peak resistance increased by 19.8% and occurred earlier in the breath. Overall pressure loss decreased by 19.2%, and the proportion of flow in the mouth increased by 13.0%. Airway wall motion was out-of-phase with the air pressure force, demonstrating the presence of neuromuscular motion. In total, the anatomy did 25.2% more work on the air than vice versa. Interpretations: Realistic movement of the airway is incorporated into CFD simulations of airflow in the upper airway for the first time. This motion is vital to producing clinically relevant computational models of respiratory airflow and will allow novel analysis of dynamic conditions, such as sleep apnoea.

1. Introduction Computational fluid dynamic (CFD) studies of airflow within the human airways have revealed much about the effect airflow has on patients' wellbeing. Previous CFD studies (Bates et al., 2014; Calmet et al., 2016; Roth et al., 2017) have explored the effect of turbulence and unsteadiness on breathing, while others (Bates et al., 2016a, 2016b; Mimouni-Benabu et al., 2012; Mihaescu et al., 2009; Palomar & Chandra, 2011; Brouns et al., 2006; Hamilton et al., 2015) have demonstrated how CFD can reveal the mechanisms by which various airway pathologies degrade subjects' health. Mylavarapu et al. ⁎

(2013) showed the utility of CFD as a surgical planning tool, comparing the predicted outcome of several virtual surgeries. However, none of the above studies incorporate the effect of airway motion into their analysis. Some previous studies that do model airway motion have been based on static images taken at particular times through the breath. One study (Yin et al., 2013) took images during three breath holds to demonstrate different pressure losses in moving airways (compared to static), and another (Subramaniam et al., 2016) used images at different applied airway pressures to estimate structural properties. Such studies illustrate the feasibility and necessity of modelling moving airways, but

Corresponding author at: Division of Pulmonary Medicine, Cincinnati Children's Hospital Medical Center, Cincinnati, OH, USA. E-mail address: [email protected] (A.J. Bates).

http://dx.doi.org/10.1016/j.clinbiomech.2017.10.011 Received 18 July 2017; Received in revised form 5 October 2017; Accepted 10 October 2017 0268-0033/ © 2017 Elsevier Ltd. All rights reserved.

Please cite this article as: Bates, A.J., Clinical Biomechanics (2017), http://dx.doi.org/10.1016/j.clinbiomech.2017.10.011

Clinical Biomechanics xxx (xxxx) xxx–xxx

A.J. Bates et al.

Fig. 1. The geometry and motion of selected planes throughout the breath cycle. The grey image shows the geometry at the start of the breathing manoeuvre (t = 0 s). The inset images on the right show cross sectional planes at three locations at six instants through the breath, as colour coded (shown with the anterior direction upmost in the plane). The planes are at fixed vertical locations and do not move with the anatomy as it moves up or down. The upper plane, at the tip of the soft palate contains holes at 0, 1 and 2 s, which is where the tip of the soft palate impinges through the plane. The middle plane shows a secondary region on the left-hand side, anterior to the epiglottis. The graph inset on the left-hand side shows the cross-sectional area of each plane through the duration of the breathing manoeuvre.

omitting the upper airway. Yin et al. (2013) used 4DCT to derive boundary conditions for CFD simulations. However due to the ionising radiation, only static CT images of the upper airway were obtained, and therefore the upper airway was assumed rigid in their simulations. The purpose of this study is to develop a new method to allow simulation of airflow in the upper airway with airway wall movement derived from MRI. We hypothesise that CFD models which incorporate more accurate motion of the upper and descending airways will provide new clinically relevant information (e.g. dynamic maps of airway sizes and resistance to flow throughout the breath) to aid in the treatment of airway diseases in which motion plays a significant role, such as obstructive sleep apnoea (OSA) and tracheomalacia. This technique will allow the relationship between aerodynamic forces on the airway walls and the motion of the airway to be categorised as either passive (when the pressure force is in the same direction as the airway motion) or active (when the pressure force is in the opposite direction to the airway wall motion, so cannot be the cause of motion). As airway movement can be significant in both heavy breathing and several pathological conditions, simulations which incorporate accurate movement increase clinical relevance by allowing analysis of airflow across a larger range of conditions than was previously possible. An understanding of the specific geometric changes of the airway in response to the various forces it experiences during respiration may, in the future, lead to structural models of the airway. For example, the specific structural properties of the tracheal rings and trachealis could inform tracheomalacia diagnoses by highlighting the reasons behind tracheal collapse.

do not exploit the full potential of novel 4D imaging techniques. Airway wall motion interpolated between static images neglects hysteresis of airway motion (Boldea et al., 2008), as well as high-speed movements, such as the vocal cords in asthmatics (Low et al., 2011) or disease related to collapsing airways. Many studies (Zhao et al., 2013a; Wang et al., 2012; Zhu et al., 2012; Pirnar et al., 2015; Ashaat & Al-Jumaily, 2016; Zhao et al., 2013b) have used fluid structure interaction (FSI) modelling to incorporate upper airway motion into CFD simulations. This approach combines flow modelling of respiratory airflow with structural modelling of the surrounding tissues to predict changes in airway shape and airflow. However, this approach assumes that the motion of the airway walls is entirely induced by internal pressure forces and neglects the effects of muscle tone and active motion. Studies that investigated the relationship between upper airway motion and the breathing cycle have found that the two are not in phase (Schwab et al., 1993a, 1993b; Cheng et al., 2011) and that upper airway dilator muscles balance low intraluminal pressure during inspiration, while the genioglossus tongue muscle moves to dilate the airway before airflow starts (Cheng et al., 2011). Airway collapse tends to happen at end expiration, rather than at peak inspiration, when the airway pressure would be most below atmospheric pressure, and therefore most likely to cause collapse. As FSI does not take active motion into account, this study adopts a new approach that incorporates airway motion observed through real-time imaging, regardless of whether the cause of motion is in response to aerodynamic forces or neuromuscular motion. Miyawaki et al. (2016a) incorporated motion into their study to demonstrate the importance of dynamic airway imaging and showed that the information provided by imaging could change the computed pressure loss within the airways by 24%, with rigid airways under-estimating pressure loss. Their moving model was based on image registration and incorporated the airway below the glottis. Jahani et al. (2015) used four-dimensional computed tomography (4DCT) images recorded over several breath cycles and retrospectively ordered to create a series of images over the breath cycle. However, the authors note 4DCT is not always obtainable due to its use of ionising radiation, a problem that can be obviated by the use of MRI, a non-ionising modality. In these models, Miyawaki et al. (2016b) demonstrated a difference of 18% in particle distribution to lobes between simulations based on static and dynamic images. Mead-Hunter et al. (2013) similarly found large discrepancies in idealised and rat lung models, again

2. Methods 2.1. MRI acquisition and image registration Static and ciné MRI was acquired on a Philips Ingenia™ 1.5T scanner (Philips Medical Systems, Best, Netherlands) using a 16-channel head and neck vascular coil. Static, high-spatial-resolution images (0.35 × 0.35 × 0.8 mm) were captured in the coronal plane using proton-density volume isotropic turbo spin echo acquisition (PDVISTA). High contrast images (3D ultrashort echo-time (UTE)) were captured and non-rigidly registered to the PD-VISTA image to delineate between the airway lumen and surrounding tissue. These images were segmented using a user-guided active contour technique (Yushkevich 2

Clinical Biomechanics xxx (xxxx) xxx–xxx

A.J. Bates et al.

→ d ′i =

et al., 2006) in ITK-Snap 3.4.0 (Penn Image Computing and Science Laboratory, USA, www.itksnap.org), yielding a high-resolution airway surface extending from the external nose and mouth to the main stem bronchi (Fig. 1). Ciné MRI was acquired using an enhanced T1 high-resolution isotropic volume excitation approach (3D e-THRIVE), which yielded 3D image-volumes at a spatial resolution of 2.37 × 2.37 × 3.00 mm and a temporal resolution of 620 ms. Non-rigid image registration was performed to extract airway movement from the ciné images using a joined, deformable motiontracking algorithm implemented in MIRTK 1.1 (Department of Computing, Imperial College London, UK, https://mirtk.github.io/). This technique produced a motion-field through the duration of the breathing manoeuvre that covered the region of the extra-thoracic airway and head. This motion-field was then applied to the high-resolution surface segmented from the static images, yielding a high-resolution surface with accurate and realistic movement throughout the breathing manoeuvre. The geometry is shown in Fig. 1. Three crosssectional planes, located at key sites within the airway (the tip of the soft palate, the tip of the epiglottis and the glottis) are shown overlaid at one-second intervals to demonstrate the motion at each location. The planes grow from their initial size (red plane and boundary) to a maximum at 2 s (blue plane), before reducing in volume. The static MRI images were obtained while the subject breathed restfully. The image acquisition was gated to the end of inspiration by respiratory bellows, hence this image represents an average anatomy at this phase of the breath taken over the 6 min of the MR scan. This image was then registered to the first in the series of ciné images to provide a high-resolution surface in the initial position of the breathing manoeuvre. As both images represent the period at the end of one exhalation, little motion was required. This moving surface provided the geometry for the CFD simulations, while the rigid geometry from the static MR images was used to determine the effect of incorporating movement in airway analysis.

N

∑j=1

→ f j (rij ) λj + → α

(1) → where N is the number of control points, λj the expansion coefficient, and fj(rij) a radial basis function of the form

f j (rij ) =

rij2 + c j2

(2) → → Here, rij = | x i − x j | represents the magnitude of distance between two vertices and cj the basis constant which is set to zero in STAR-CCM → +. For large displacements, the constant vector α in Eq. (1) is used to → satisfy the additional constraint on λj , such that N

∑j=1

→ λj = 0

(3)

Eqs. (1) and (3) are solved together to provide Cartesian compo→ → nents for all λj and α . From these vectors, using multi-quadratic interpolation theory (Hardy, 1990) an interpolation field,

d′ (→ x)=

N

∑j=1

→ f j (r ) λj + → α

(4) ′

for the whole domain is generated, where d is the calculated displacement applied at each node defined by the position vector → x . This interpolation field is used to move all floating mesh vertices in the domain. After each time-step, if the volume of any cell in the mesh became negative, or the quality metrics dropped below certain thresholds (cell quality < 1 × 10− 5 (a function of the relative geometric distribution of the centroids of the neighbouring cells and the orientation of the cell faces (STAR-CCM+ User Guide, 2016)), face validity < 0.73 (an areaweighted measure of the direction of a cell's face-normal vectors relative to the position of the cell centroid (STAR-CCM+ User Guide, 2016))), the simulation was reverted back to its solution state and topology at the previous time-step. This ensures that the morphing procedure does not degrade the mesh to the point whereby the flow solution would suffer. A new mesh is then generated from the topology, with the solution state mapped over from old to the new mesh. The simulation is then allowed to continue with the above process. Since data-mapping introduces numerical artefacts associated with interpolating solution data from the old to new mesh, additional inner iterations were performed by the flow solvers in subsequent time-steps to asymptotically converge residuals back down to normal levels of the simulation. The continuity and momentum Navier-Stokes equations for a moving mesh are of the finite volume form (STAR-CCM+ User Guide, 2016)

2.2. Subject and breathing manoeuvre A healthy, 31-year-old, male volunteer was imaged with informed consent and IRB approval. During imaging, the volunteer wore a size-5 anaesthesia face mask (Vital Signs, GE Healthcare, Totowa, New Jersey, USA), connected to an MRI-compatible pneumotachograph (ADInstruments MLT300L Respiratory Flow Head and Spirometer) (Eichinger et al., 2007) to record the breathing flow rate. This data was synchronised with pulses from the MRI scanner at each cine image acquisition. The volunteer was asked to perform a spirometry manoeuvre (a deep breath in followed by a fast breath out) at non-maximal effort and was free to breathe through his nose, mouth or any combination of the two. Previous tests showed that this manoeuvre produced a similar range of motion in the healthy, awake volunteer to that seen in sleeping patients with sleep apnoea.

∂ ∂t

⎯v )·d→ a =0 ∫V∼ ρdV∼ + ∮ ρ (→v − ⎯→ g

(5)

and

∂ ∂t

⎯v ) ⊗ → v )·d→ a = − ∮ pI ·d→ a + ∮ T·d→ a ∫V∼ ρ→v dV∼ + ∮ (ρ (→v − ⎯→ g (6)

2.3. Computational flow simulation

respectively, where t is time, V is the volume of each cell in the mesh, ρ ⎯v is the mesh velocity as calis the air density, → v is the air velocity, ⎯→ g → culated from the mesh displacements in Eq. (1), a is a vector representing the surface of each mesh cell, I is the identity matrix, and T is the viscous stress tensor.

2.3.1. Moving wall CFD methodology The Navier-Stokes equations were used to calculate the air pressure and velocity throughout the domain, bounded by moving walls, using STAR-CCM+ 11.04.012 (Siemens PLM Software, Plano, TX, USA). Control points were distributed across the topology of the surface spaced 2 mm apart and a series of co-ordinates representing the new position of the control points at each time-step throughout the breath were provided to the mesh morpher. Automatic thinning was applied to reduce the number of control points and speed up the execution time of the morpher. → The displacement d ′i , for each control point i is expanded as (STARCCM+ User Guide, 2016),

2.3.2. Mesh generation The same meshing approach was used for both the initial mesh and all subsequent re-meshes throughout the morphing procedure. The volume mesh was comprised of polyhedral cells in the core and five layers of prismatic cells to resolve the boundary layer at the walls. A polyhedral cell is an unstructured mesh topology, whereby each cell can have an arbitrary number of faces. This generally increases cell connectivity to neighbours as well as the degree of freedom available to the 3

Clinical Biomechanics xxx (xxxx) xxx–xxx

A.J. Bates et al.

Fig. 2. The surface of the airway model at six instants through the breathing manoeuvre. The model extends from the mask worn by the subject to the main stem bronchi. The anatomy has changed position at each time-point: the soft palate moves anteriorly between 1 and 4 s. The tongue moves back initially, then forward again by the end of the breathing manoeuvre, changing the calibre of the oral airway throughout the breath. The interior of each image is coloured by volume rendering of the velocity field at each instant. The velocity throughout the geometry is shown, with higher values represented with high opacities while lower values are more transparent. The top three images are during inhalation, while the lower three represent exhalation. An animation of the moving geometry with volume rendered velocity is included in the Supplementary materials.

mesher to resolve topological features. As a result, higher quality cells can be obtained with low skewness and smooth volumetric growth, characteristics that facilitate ease of convergence. Consequently, the polyhedral mesh technology generally shows better convergence than a tetrahedral mesh for a given number of cells (Peric, 2004). The number of cells in the mesh averaged 3.6 million cells, with a mean cell volume of 0.064 mm3. A mesh and time-step convergence study was carried out with mesh base sizes ranging from 0.5 to 2.0 mm and time-steps between 2 × 10− 3 s and 1 × 10− 4 s. The simulation with a mesh base size of 1 mm (1.1 million cells) and a time-step of and 1 × 10− 3 s gave peak inhalation pressure and resistance values to within 0.26% of the finest simulation (3.6 million cells), demonstrating convergence at this point. Meshes below this resolution had significantly different peak pressure values (up to 22.4% different in the coarsest case).

As the name implies, LES aims to resolve large-scale motions, and model instead only the filtered small-scale motions of the turbulent flow field. In STAR-CCM+, this filtering is a function of local grid size, and each increase in mesh refinement leads to greater resolution of the small-scale motions (Zhiyin, 2015). A consequence of this means that a truly mesh independent LES solution becomes the equivalent of a DNS, hence no mesh independency study is really required. Instead, typical approaches aim to resolve at least 80% of the turbulent kinetic energy (Pope, 2000), with a level of mesh size refinement chosen that sufficiently resolves the Taylor length scale (Agarwal et al., 2014; Addad et al., 2008). The Taylor length scale, λ, marks the point at which fluid viscosity begins to significantly affect the dynamics of turbulent eddies and is approximated as

2.3.3. Flow calculation methodology and boundary conditions The Navier-Stokes equations (Eqs. (5) & (6)) were solved using the large eddy simulation (LES) techniques that have previously been verified against near-direct numerical simulation (DNS) results (Bates et al., 2017).

where ν is the kinematic viscosity, k the turbulent kinetic energy, and ε the turbulence dissipation rate. These quantities were readily computed from a pre-curser URANS simulation, with measurements taken at specific locations within the airway to determine local mesh size requirements for the LES computation.

λ≈

4

10ν

k ε

(7)

Clinical Biomechanics xxx (xxxx) xxx–xxx

A.J. Bates et al.

curves) and static (dashed curves) geometries further reveals the effect motion has on airway resistance. Initially (0–0.65 s), there is little difference between the resistances in the two cases as the two geometries are identical at time (t) = 0 s, and there is little movement in the first half of inhalation (see the relatively similar planes for 0 and 1 s in Fig. 1 and the small initial change in airway volume shown in Fig. 4). Once the airway starts to move, the resistance of the mouth and nose increases sharply from the static case, with a maximum increase of 76% at 1.15 s. This is due to the tongue, which initially moves backwards, narrowing the rear oral airway and forcing more air through the higherresistance path of the nose. In the static simulation, the proportion of air travelling through the oral (as opposed to nasal) airway is relatively constant through inhalation at 53%, in the dynamic case, this value initially falls to 50%, before rising to 66% later in the inhalation (t = 1.38 s). The rise in the proportion of airflow travelling through the mouth (Fig. 5) coincides with the increase in airway volume shown in Fig. 4 (blue curve). The opening of the jaw to widen the oral airway causes this increase in airway volume and increase of oral airflow (Fig. 5). Once this has occurred, the resistance of the nose and mouth in the moving geometry is significantly less (84% at peak exhalation, t = 3.22 s) than in the static case and remains so throughout exhalation. The resistance of the pharynx is much higher in exhalation than inhalation, and more markedly so in the static case than the dynamic due to the way its changing shape influences the impingement of the jet formed by the narrowing at the glottis, as can be seen in the video in the Supplementary material. The resistance of the whole geometry decreased by 10.7% in inspiration and increased by 19.8% during expiration.

Second order schemes were used for both spatial and temporal discretisation with a time-step of 1 × 10− 4 s. No-slip conditions were applied at all airway walls, as well as surfaces of the subject's face and mask. The mask inlet was held at atmospheric pressure, while the massflow rate recorded by the pneumotachograph was applied at the ends of the bronchi, split 55–45% to the right and left respectively, based on the subject's lung volumes. At both the mask inlet and distal bronchi, any incoming flow was uniformly distributed across the inlet surface. By the time the flow has reached the thoracic inlet and the area of interest, any effects of an incorrect inlet spatial velocity profile will have dissipated (Tadjfar & Smith, 2004). The instantaneous resistance is calculated as the pressure loss between two locations divided by the flow rate through them. For localised regions, the two locations are the beginning and end of that section, for the overall resistance, it is the mask inlet to the bronchi. 2.3.4. Wall force calculations The aerodynamic forces on the face on one mesh element of the airway wall Fw, are made up of the internal pressure at the wall p, which always acts on a wall in the direction normal to the surface → n , and the τ , and can be expressed as wall shear stress →

→+→ Fw = A (pn τ)

(8)

where A is the area of the face. The dot-product of this force and the wall velocity, Uw give a power measurement quantifying the rate of work by the aerodynamic forces on the airway wall, if positive, or the rate of work done by the airway wall on the air if negative;

P = Uw·Fw

(9)

Summing these values over every face on the airway surface mesh gives the total power exchange between the airway wall and the air.

3.1.2. Force/motion correlation The correlation between the motion of the airway wall and the direction of the force the airflow applies to the wall can reveal the causes of airway movement. Fig. 6 shows the airway at two instants from the breathing manoeuvre - peak inhalation and exhalation. The red arrows show the motion of the airway walls at this time (e.g. the posterior aspect of the tongue is moving down and forwards at peak inhalation). The blue arrows show the aerodynamic force applied to the airway at the same two instants. These are the vector sum of the pressure and wall shear stress forces and largely point into the lumen during the low pressure of inspiration and outwards during exhalation. Taking the dot-product of the force and motion vectors (Eq. (9)) reveals the correlation between the movements of the airway wall (determined by registration of the MR images) and the aerodynamic forces it is subjected to (determined by CFD simulation). The third pair of images in Fig. 6 shows this correlation, with positive values indicating the movement and force vectors point in the same direction and negative showing they point in opposite direction. At both peak inhalation and peak exhalation, little of the airway surface motion is in same direction as the force vectors (red) and this is true for intermediate times too. Most of the airway movement is in the opposite direction to the applied aerodynamic force (blue) or there is no discernable relationship between the two (white). This last case represents when the force and motion are at 90° to each other, the airway wall is not moving, or there is no force applied.

3. Results & discussion 3.1. Flow analysis The position of the airway surface is shown at one-second intervals throughout the breathing manoeuvre in Fig. 2, coloured by a volume rendering of the air's velocity field. The motion of the tongue, soft palate, epiglottis and glottis are all shown by the different positions of the geometry at each instant. During inhalation (upper three images), highvelocity regions occur in the mouth, nasopharynx, in a jet emanating from the constriction posterior to the epiglottis and the glottal jet. The incoming flow rates at 1 and 2 s are similar (1.78 and 1.87 l/s respectively), yet there is much more high-velocity flow at 1 s. The motion of the jaw and tongue has widened the airway at 2 s, as can be seen from the large planes from this instance in Fig. 1 (blue). This widening of the airway allows the air to travel at a lower velocity, despite the slightly higher flow rate, thereby reducing resistance (discussed below). During exhalation (lower images, Fig. 2), the region of fastest airflow is at the tip of the epiglottis and around the tip of the soft palate. 3.1.1. Resistance to airflow Fig. 3 shows the resistance to airflow through several different regions of the extrathoracic airway throughout the breath, demonstrating the contribution of each region to the overall resistance to airflow. The variation of resistance through each region with time is due to the airflow rate (black curve, Fig. 3) through the breath and the movement of the airway. During inspiration (0–2.26 s), the major source of resistance to airflow in the extrathoracic airway is the mouth and nose (solid blue curve), with resistance peaking at 4.23 × 10− 3 cm H2O·s/ml after 1.15 s of inhalation. This peak resistance occurs earlier than the peak flow rate (1.72 s into the manoeuvre), indicating that the motion of the airway influences the resistance more than the flow rate at this time. Comparison between the resistance in the moving (Fig. 3, solid

3.1.3. Power transfer In order to determine how much of the airway movement is in phase with the aerodynamic forces over the entire extrathoracic airway surface and throughout the duration of the breathing manoeuvre, the surface integral of the force-motion dot-product is taken. This calculation gives the total power transfer between the airway wall and the intraluminal airflow. The positive and negative contributions of this measure are plotted in Fig. 4. There is little motion in the first 0.5 s of inhalation, perhaps as the volunteer had slightly opened his mouth in preparation for the 5

Clinical Biomechanics xxx (xxxx) xxx–xxx

A.J. Bates et al.

Fig. 3. The resistance to airflow through the breath. The coloured lines represent the resistance (left axis) through each of the regions between the planes shown in the inset (top left). The solid lines show the resistance in the moving wall simulation, while the dashed lines show the resistance in the same regions in the static geometry. The black curve shows the flow rate throughout the breath (right axis). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4. The power transfer between the airway wall and the intraluminal airflow. The solid red curve shows the work done by the airflow on the airway wall (defined as positive). The dashed red curve shows the work done by the extrathoracic airway anatomy on the airflow to move the airway walls (defined as negative). The blue curve shows the airway volume and the black curve shows the total pressure difference between the ends of the main stem bronchi and the atmospheric pressure at the inlet to the mask. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5. The total flow rate of air inhaled and exhaled by the subject during the breathing manoeuvre (blue curve). The red, yellow and purple curves show how this flow is divided between the left and right nostrils, and the mouth respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

(Fig. 6). Therefore, during this period, there is a large increase in out-ofphase power transfer, meaning the anatomy is doing work on the fluid to expand the airway against the pressure force. There is also a small increase in in-phase work during this period, as the tip of the soft palate moves down, partially closing the oral airway.

breathing manoeuvre. Therefore, there is little power transfer between the air and the wall. There is a period of rapid airway expansion between 1 and 2 s, when the airway expands by 10 cm3. However, during inspiration, the intraluminal pressure is below atmospheric and consequently, the force vectors contract rather than expand the airway 6

Clinical Biomechanics xxx (xxxx) xxx–xxx

A.J. Bates et al.

Fig. 6. The movement and aerodynamic forces acting on the airway at two instants from the breathing manoeuvre. The top row shows peak inhalation, the bottom row peak exhalation. The red arrows in the left pair of images show the movement vectors of the airway surface at the respective instants. The blue arrow in the central images show the aerodynamic force vectors (pressure force plus wall shear stress) acting on the airway wall. The right-hand images show the dot-product between these two vectors, indicating the correlation between aerodynamic force and motion in the extrathoracic airways. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

motion and its causes in OSA subjects. The technique implemented here has quantified the effect of incorporating motion into CFD simulations of airflow in the extrathoracic airway. Comparisons between simulations with moving and static geometry demonstrated large differences in key metrics: a 14.6% increase in peak pressure loss (mask inlet to carina) during inspiration, a 19.2% decrease during exhalation and a maximum of 13% change in the distribution of airflow through the oral rather than nasal geometry. These changes in pressure losses are similar to those found by Miyawaki et al. (2016a) lower down the airway. The resistance to inspiratory flow decreased by 10.7% and to expiratory flow, increased by 19.8%. Not only did the extreme values change, but the time at which they occurred changed too, despite the same flow profile, demonstrating that the anatomical movement had a bigger effect on the resistance than the changing flow rate throughout the breath. This effect is shown in the curves representing local resistance in Fig. 3, while the extent of the changes in airway lumen is shown by the cross-sectional planes shown in Fig. 1, the geometries at different instants in Fig. 2 and the airway volume curve in Fig. 4. The airway expanded once the rapid-breathing manoeuvre started, perhaps as a reaction to the manoeuvre and demonstrating that these techniques may allow clinicians to assess the body's control of the airway in response to airflow. Much of the airway wall motion observed in this study was found to be active motion, i.e. not caused by the intraluminal aerodynamic forces, and therefore likely to be caused by neuromuscular activity. Previous studies that have incorporated upper airway movement into the CFD model have used fluid-structure interaction, a technique which only models passive rather than active

At the crossover between inspiration and expiration, the intraluminal forces drop to zero, yet the airway continues to move, again demonstrating movement caused by neuromuscular activity rather than aerodynamic forces. During exhalation, the anatomy and air do similar amounts of work on each other, as different regions of the airway move in and out of phase from the force. The airway volume decreases significantly (by 19 cm3) from 3.6 to 5.0 s, which is the end of the exhalation period, but as the pressure forces in the extrathoracic airway are relatively small during this period, little work is done by the anatomy on the air to collapse the airway. In total, the anatomy does 25.2% more work on the air than the air does to move the anatomy.

4. Conclusions This study describes the first simulations of breathing in the extrathoracic airway to incorporate both passive and neuromuscular movement of the airway. CFD simulations have previously shown clinical relevance in a range of different conditions (Bates et al., 2016b; Hamilton et al., 2015; Mylavarapu et al., 2013), but without incorporating this motion, were unable to model morbidities that include motion, such as OSA or tracheomalacia. This work has demonstrated the feasibility of combining static and ciné MR images to produce airway geometries for CFD simulations that move realistically, with high-spatial and temporal resolution and without the use of ionising radiation. The range of motion exhibited by the volunteer in this study was similar to that seen in pediatric patients with OSA, indicating that this method may be an ideal way to assess 7

Clinical Biomechanics xxx (xxxx) xxx–xxx

A.J. Bates et al.

simulations of gas turbine combustors. In: Combust. Fuels Emiss. 4B ASME (p V04BT04A011). Ashaat, S., Al-Jumaily, A.M., 2016. Reducing upper airway collapse at lower continuous positive airway titration pressure. J. Biomech. 49, 3915–3922. Bates, A.J., Doorly, D.J., Cetto, R., Calmet, H., Gambaruto, A.M., Tolley, N.S., Houzeaux, G., Schroter, R.C., 2014. Dynamics of airflow in a short inhalation. J. R. Soc. Interface 12, 20140880. Bates, A.J., Comerford, A., Cetto, R., Schroter, R.C., Tolley, N.S., Doorly, D.J., 2016a. Power loss mechanisms in pathological tracheas. J. Biomech. 49, 2187–2192. Bates, A.J., Cetto, R., Doorly, D.J., Schroter, R.C., Tolley, N.S., Comerford, A., 2016b. The effects of curvature and constriction on airflow and energy loss in pathological tracheas. Respir. Physiol. Neurobiol. 234, 69–78. Bates, A.J., Comerford, A., Cetto, R., Doorly, D.J., Schroter, R.C., Tolley, N.S., 2017. Computational fluid dynamics benchmark dataset of airflow in tracheas. Data Brief 10, 101–107. Boldea, V., Sharp, G.C., Jiang, S.B., Sarrut, D., 2008. 4D-CT lung motion estimation with deformable registration: quantification of motion nonlinearity and hysteresis. Med. Phys. 35, 1008–1018. Brouns, M., Jayaraju, S.T., Lacor, C., De Mey, J., Noppen, M., Vincken, W., Verbanck, S., 2006. Tracheal stenosis: a flow dynamics study. J. Appl. Physiol. 102, 1178–1184. Calmet, H., Gambaruto, A.M., Bates, A.J., Vázquez, M., Houzeaux, G., Doorly, D.J., 2016. Large-scale CFD simulations of the transitional and turbulent regime for the large human airways during rapid inhalation. Comput. Biol. Med. 69, 166–180. Cheng, S., Butler, J.E., Gandevia, S.C., Bilston, L.E., 2011. Movement of the human upper airway during inspiration with and without inspiratory resistive loading. J. Appl. Physiol. 110, 69–75. Eichinger, M., Puderbach, M., Smith, H.J., Tetzlaff, R., Kopp-Schneider, A., Bock, M., Biederer, J., Kauczor, H.U., 2007. Magnetic resonance-compatible-spirometry: principle, technical evaluation and application. Eur. Respir. J. 30, 972–979. Hamilton, N.J., Kanani, M., Roebuck, D.J., Hewitt, R.J., Cetto, R., Culme-Seymour, E.J., Toll, E., Bates, A.J., Comerford, A.P., McLaren, C.A., Butler, C.R., Crowley, C., McIntyre, D., Sebire, N.J., Janes, S.M., O'Callaghan, C., Mason, C., De Coppi, P., Lowdell, M.W., Elliott, M.J., Birchall, M.A., 2015. Tissue-engineered tracheal replacement in a child: a 4-year follow-up study. Am. J. Transplant. 15, 2750–2757. Hardy, R.L., 1990. Theory and applications of the multiquadric-biharmonic method 20 years of discovery 1968–1988. Comput. Math. Appl. 19, 163–208. Jahani, N., Choi, S., Choi, J., Iyer, K., Hoffman, E.A., Lin, C.-L., 2015. Assessment of regional ventilation and deformation using 4D-CT imaging for healthy human lungs during tidal breathing. J. Appl. Physiol. 119, 1064–1074. Low, K., Lau, K.K., Holmes, P., Crossett, M., Vallance, N., Phyland, D., Hamza, K., Hamilton, G., Bardin, P.G., 2011. Abnormal vocal cord function in difficult-to-treat asthma. Am. J. Respir. Crit. Care Med. 184, 50–56. Mead-Hunter, R., King, A.J.C., Larcombe, A.N., Mullins, B.J., 2013. The influence of moving walls on respiratory aerosol deposition modelling. J. Aerosol Sci. 64, 48–59. Mihaescu, M., Gutmark, E.J., Elluru, R., Willging, P.J., 2009. In: Large eddy simulation of the flow in a pediatric airway with subglottic stenosis. 47th AIAA Aerosp Sci Meet Incl New Horizons Forum Aerosp Expo 3333. Mimouni-Benabu, O., Meister, L., Giordano, J., Fayoux, P., Loundon, N., Triglia, J.M., Nicollas, R., 2012. A preliminary study of computer assisted evaluation of congenital tracheal stenosis: a new tool for surgical decision-making. Int. J. Pediatr. Otorhinolaryngol. 76, 1552–1557. Miyawaki, S., Choi, S., Hoffman, E.A., Lin, C.L., 2016a. A 4DCT imaging-based breathing lung model with relative hysteresis. J. Comput. Phys. 326, 76–90. Miyawaki, S., Hoffman, E.A., Lin, C.L., 2016b. Effect of static vs. dynamic imaging on particle transport in CT-based numerical models of human central airways. J. Aerosol Sci. 100, 129–139. Mylavarapu, G., Mihaescu, M., Fuchs, L., Papatziamos, G., Gutmark, E., 2013. Planning human upper airway surgery using computational fluid dynamics. J. Biomech. 46, 1979–1986. Palomar, A.P., Chandra, S., 2011. FSI analysis of a healthy and a stenotic human trachea under impedance-based boundary. J. Biomech. Eng. 133, 1–12. Peric, M., 2004. Flow simulation using control volumes of arbitrary polyhedral shape. ERCOFTAC Bull. 62, 25–29. Pirnar, J., Dolenc-Grošelj, L., Fajdiga, I., Žun, I., 2015. Computational fluid-structure interaction simulation of airflow in the human upper airway. J. Biomech. 48, 3685–3691. Pope, S.B., 2000. Turbulent flows. J. Turbul. 1, 771. Roth, C.J., Yoshihara, L., Ismail, M., Wall, W.A., 2017. Computational modelling of the respiratory system: discussion of coupled modelling approaches and two recent extensions. Comput. Methods Appl. Mech. Eng. 314, 473–493. Schwab, R.J., Gefter, W.B., Pack, A.I., Hoffman, E.A., 1993a. Dynamic imaging of the upper airway during respiration in normal subjects. J. Appl. Physiol. 74, 1504–1514. Schwab, R.J., Gefter, W.B., Hoffman, E.A., Gupta, K.B., Pack, A.I., 1993b. Dynamic upper airway imaging during awake respiration in normal subjects and patients with sleep disordered breathing. Am. Rev. Respir. Dis. 148, 1385–1400. STAR-CCM+ User Guide. Siemens PLM Software, Plano, TX, USA. Subramaniam, D.R., Mylavarapu, G., McConnell, K., Fleck, R.J., Shott, S.R., Amin, R.S., Gutmark, E.J., 2016. Compliance measurements of the upper airway in pediatric down syndrome sleep apnea patients. Ann. Biomed. Eng. 44, 873–885. Tadjfar, M., Smith, F., 2004. Direct simulations and modelling of basic three-dimensional bifurcating tube flows. J. Fluid Mech. 519, 1–32. Wang, Y., Wang, J., Liu, Y., Yu, S., Sun, X., Li, S., Shen, S., Zhao, W., 2012. Fluid-structure interaction modeling of upper airways before and after nasal surgery for obstructive sleep apnea. Int. J. Numer. Methods Biomed. Eng. 28, 528–546. Yin, Y., Choi, J., Hoffman, E.A., Tawhai, M.H., Lin, C.L., 2013. A multiscale MDCT imagebased breathing lung model with time-varying regional ventilation. J. Comput. Phys.

motion. In the breathing manoeuvre analysed, the airway anatomy did 25.2% more work than the intraluminal airflow in moving the airway wall, showing the importance of incorporating both types of motion into CFD simulations in order to produce physiologically realistic results. In conditions involving excessive airway narrowing such as OSA, the ratio of power transfer between the anatomy and airflow may be an indicator which would allow different phenotypes of OSA to be identified, along with patient-specific surgical treatments. The temporal resolution of fast MRI does not allow for very fast dynamic events such as speech generation to be analysed, although the temporal resolution is sufficient to capture several images through a breath cycle, revealing the major changes in airway calibre. While only one healthy volunteer was considered in this study, the success of the technique represents a step forward in allowing CFD to be used to assess morbidities of the upper and descending airways. Any individual's unique airway anatomy and motion can be modelled similarly, providing a subject-specific approach applicable to a broad range of airway pathologies. The large differences between static and dynamic simulations found in this study all occurred during a breathing manoeuvre using much higher flow rates than would be normal in restful breathing. Future studies should address whether similar results are found in OSA patients with lower breathing flow rates, but similar ranges of motion. The techniques developed in this study are designed to provide information about short periods of breathing such as a few breaths, for example identifying why an OSA patient's airway collapses on some breaths but not others. In analysing a single breath of a particular subject, a 20% difference in resistance represents a significant discrepancy. However, differences of this magnitude are small in comparison to those found between different breathing regimes or across different subjects. For future static simulations, this also raises the question of which static geometry should be used. Many clinical scans are captured at end expiration and this may not be an appropriate static geometry for simulating peak inhalation or exhalation breathing. In order to utilise this method to extract clinically useful information, future studies should investigate the similarity between multiple breaths by the same patient, and how different breathing patterns affect motion and the causes of motion. The range of motion in normal subjects should be investigated and compared to patients with dynamic airway diseases. In conclusion, this study has demonstrated a novel method for computing the airflow of breathing within the extrathoracic airway with prescribed wall movement, calculated from non-rigid registration of ciné MR images. The addition of motion into existing CFD airway models widens their applicability to a much broader range of breathing manoeuvres and pathological conditions than was previously possible. Acknowledgments The authors would like to express their thanks to Brynne M. Williams RT(MR) CNMT and J. Matthew Lanier RT(R)(MR) for their help and expertise in capturing the MR images. The authors acknowledge financial support from Imperial College London through a Junior Research Fellowship grant awarded to the lead author. Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.clinbiomech.2017.10.011. References Addad, Y., Gaitonde, U., Laurence, D., Rolfo, S., 2008. Optimal unstructured meshing for large eddy simulations. In: Qual. Reliab. Large-Eddy Simulations. Springer Netherlands, Dordrecht, pp. 93–103. Agarwal, K.K., Konakalla, S.R., Selvan, S., 2014. On meshing guidelines for large eddy

8

Clinical Biomechanics xxx (xxxx) xxx–xxx

A.J. Bates et al.

treatment response of oral appliances for obstructive sleep apnea using computational fluid dynamics and fluid-structure interaction simulations. Biomed. Biotechnol. Eng. 3A (p V03AT03A051). Zhiyin, Y., 2015. Large-eddy simulation: past, present and the future. Chin. J. Aeronaut. 28, 11–24. Zhu, J.H., Lee, H.P., Lim, K.M., Lee, S.J., Teo, L.S.L., Wang, D.Y., 2012. Passive movement of human soft palate during respiration: a simulation of 3D fluid/structure interaction. J. Biomech. 45, 1992–2000.

244, 168–192. Yushkevich, P.A., Piven, J., Hazlett, H.C., Smith, R.G., Ho, S., Gee, J.C., Gerig, G., 2006. User-guided 3D active contour segmentation of anatomical structures: significantly improved efficiency and reliability. NeuroImage 31, 1116–1128. Zhao, M., Barber, T., Cistulli, P.A., Sutherland, K., Rosengarten, G., 2013a. Simulation of upper airway occlusion without and with mandibular advancement in obstructive sleep apnea using fluid-structure interaction. J. Biomech. 46, 2586–2592. Zhao, M., Barber, T., Cistulli, P., Sutherland, K., Rosengarten, G., 2013b. Predicting the

9