We analyzed airflow in two species of ankylosaur, the nodosaurid Panoplosaurus mirus (ROM 1215) and the ankylosaurid Euoplocephalus tutus (AMNH 5405). CT data and initial 3D models were obtained from previous work conducted by Witmer & Ridgely [ 20 ]. To aid with soft-tissue reconstruction, we further looked at other specimens of Euoplocephalus (AMNH 5403, ROM 1930) and Panoplosaurus (CMN 2759) as well as the related species: Edmontonia rugosidens (AMNH 5381), Ankylosaurus magniventris (AMNH 5214), Pinacosaurus grangeri (ZPAL MgD-II/1), and Kunbarrasaurus ieversi (QM F18101).

Model construction

Segmentation. Osteological evidence for a complete nasal septum (septal sulcus or partially mineralized septum) meant that left and right nasal passages acted independently of each other. This allowed for modeling of one side of the nose only, which saved on computational costs. Initial segmentation of the airways produced a rough approximation of the nasal passage in life, complete with a rostral and caudal vestibular loop [20] and an enlarged olfactory recess (Fig 2). We refer to the enlarged, looping area of the nasal passages in both ankylosaurs as the nasal vestibule. This demarcation of nasal passage anatomy in sauropsids is typically determined by the placement of the duct for the nasal gland [105,106]. Unfortunately, the duct and its ostium are both soft-tissue structures that do not leave impressions on the bone. As an alternative, we used the region of the nasal passage where nasal cavity diameter suddenly increased [110]. This area is known to correlate with the terminus of the nasal vestibule in extant sauropsids.[111] Further, support for this interpretation comes from the tube-shaped structure of the nasal vestibule in extant reptiles. The morphology of the looping portions of the ankylosaur airways best fits this description. The nasal vestibule and the CNP were unattached to each other in the original airway segmentations [20], requiring further segmentation and attachment to produce a contiguous surface model. The nasopharyngeal duct was not segmented in the original models and required segmentation. Although Witmer and Ridgely [20], as well as Miyashita et al. [21], discussed the presence of well-preserved olfactory turbinates within the olfactory recess of these dinosaurs, neither study published images of the segmented structures. For our study, these structures and their effect on the airway (i.e., the impressions they left on the digital airway cast) were segmented using the program Avizo 7.1 (FEI Visualization Sciences Group, Burlington, MA). As with the initial airway segmentation, the final product was a cast of the inside of the nasal cavity, revealing the potential space in which air could reside within the nasal passage. Examination of the CT data within the olfactory recess revealed both the presence of mineralized olfactory turbinates as well as the outer boundary to the nasal capsule. These data provided insight into the limit of the airway in life, which was substantially more restricted than initial segmentations suggested (Fig 21). Extra space lateral to our interpretations of the nasal passage wall was interpreted as housing the antorbital sinus. Its placement near the olfactory chamber, adjacent to the CNP, is consistent with antorbital sinus placement in extant archosaurs [14]. In preparation for meshing, the segmented models were cleaned of segmentation artifacts using the program Geomagic 10 (3D Systems Geomagic, Rock Hill, SC). PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 21. Segmentation of the airway in Panoplosaurus mirus (ROM 1215). (A) Skull in left lateral view. Line represents the location of (B) axial CT section showing preserved olfactory turbinates. (C) Diagram of CT image showing caliber of airway vs. entire nasal cavity. (D) Segmented olfactory turbinate in same plane as CT image. https://doi.org/10.1371/journal.pone.0207381.g021

Fleshy nostril placement and soft palate. Although ROM 1215 and AMNH 5405 both contained well preserved nasal passages, the terminal regions of the nose—the fleshy nostril and choana—were not preserved. To aid with these soft-tissue reconstructions we turned to other specimens of the same species along with well-preserved specimens of related ankylosaurs to better determine the location of the nostril and choana.

Fleshy nostril. For Panoplosaurus, we used the well-preserved skull of CMN 2759 to determine the location of the nostril. CMN 2759 preserved the rostral wall of the bony narial aperture, which was comprised of the rostral-most cranial osteoderm, most similar to the median nasal caputegulum of Euoplocephalus [112]. These anatomical structures strongly suggested that the nostril of Panoplosaurus deviated laterally (Fig 22A and 22C). Laterally-facing nostrils are common among diapsids and such a placement in Panoplosaurus was not unexpected. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 22. Nostril placement in ankylosaur models. All skulls in right lateral and rostral views. For Panoplosaurus mirus (A, C) we used (A) ROM 1215 as our base model with nostril placement informed by (C) CMN 8530. For Euoplocephalus tutus (B, D) we used (B) AMNH 5405 as our base model with the skull of (D) ROM 1930 informing us on the limits to the extent of the nostril. Asterisks in A and B denote location of fleshy nostril in our models. https://doi.org/10.1371/journal.pone.0207381.g022 For Euoplocephalus, AMNH 5405 did not present a well-preserved rostral-most portion of the cranium. To assist with nostril positioning, we compared bony narial aperture shape in AMNH 5405 with ROM 1930 (Fig 22B and 22D). In the latter, preservation of the rostral-most region of the cranium revealed an enlarged bony narial aperture that faced forward suggesting terminal nostril placement. However, the width of the bony narial aperture encompassed both the rostral-most portion of the cranium as well as a lateral portion of the cranium. Thus, it remains possible that the nostril in Euoplocephalus could have deviated laterally, or had a rostrolateral combination thereof. This position would be consistent with fleshy nostril placement in Ankylosaurus (AMNH 5214) where dermal ossification is extensive and indicates unambiguously that the nostril was located rostroventrolaterally [77]. In contrast, Anodontosaurus lambei (CMN 8530), a close relative of Euoplocephalus, shows a well constrained bony narial aperture that would limit the nostril to a terminal position on the snout [112]. All known skulls of Euoplocephalus that preserve the rostral-tip of the snout show a much less constrained bony narial aperture. This could indicate that the caputegula covering the cranium were less extensive in this species and that terminal fleshy nostrils were present but were constrained only by soft-tissues. AMNH 5405 has a fossa on the premaxillae (apertura nasalis [17]) that has previously been interpreted as a portion of the nasal vestibule ([17,20], Fig 21B and 21D). This fossa is wide enough that it could have easily housed a laterally deviating nasal vestibule that terminated rostrally in a laterally-facing nostril. For the purposes of our analysis we fit Euoplocephalus with such a nostril, with the caveat that the osteological evidence for it was equivocal (Fig 22B and 22D). Such a position is consistent with the general finding in amniotes that fleshy nostrils tend to be rostroventrally situated within the nasal vestibule [113].

Soft-tissue correction. Airways segmented from the skulls represented the outer limits of the nasal passage, referred to here as the bony-bounded (BB) airway (Figs 15C and 16C). In life, soft tissues within the nasal passage such as the nasal capsule cartilages, mucosa, nerves, and vasculature would have been present and would have occupied space within these well-constrained nasal passages (Fig 26). Previous surveys of airway calibers in the nasal passages of extant amniotes found that airway caliber (the distance spanned between mucosal walls) does not exceed 10 mm in diameter regardless of the size (0.1–600kg) or phylogenetic position (from squamate reptiles to artiodactyl mammals) of the animal [33,80]. This constraint appears to be dictated by the biophysical limitations of diffusion, which effectively works across only very small distances [71]. The thermoregulatory and olfactory functions of the nasal cavity are both diffusion-dependent functions [29,31,126]. In contrast to data from extant animals, the average BB airway calibers in Panoplosaurus and Euoplocephalus were 15.8 mm and 22.9 mm, respectively. These were significantly larger calibers than that observed in the mucosa-lined airways of extant amniotes. To bring airway calibers within the range observed in extant amniotes, we imported airway models of Panoplosaurus and Euoplocephalus into the 3D modeling and animation program Maya (Autodesk, San Rafael, CA) where the airways were compressed using the program’s 3D deformation tools. Airway calibers were reduced until the average diameter was ~10 mm, which is the upper limit observed in extant amniotes, making this a conservative estimate for ankylosaurs. Airway compression followed the natural contours of the nasal passage such that the refined airways resembled a more compressed version of the original segmentation (Figs 15 and 16). These soft-tissue-corrected airways are here referred to as the ST airways (Figs 15D and 16D). An extension off the choana was added to all models tested. This extension served to replicate the connection of the nasal airway to the larynx and trachea. This extension was added to address technical aspects of the software (see below) and to ensure that fully developed airflow would be present at the choana during expiration, thus removing any potential heat flow artifacts created by having the program initialize airflow at the choana during expiration (Figs 15 and 16). PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 26. Generic airway diagram for diapsids. Note the much more constricted airway in the soft-tissue nasal passage (A) vs. the emptier, bony-bounded nasal passage typically preserved in fossils (B). https://doi.org/10.1371/journal.pone.0207381.g026 To further test the hypothesis that the convoluted airways in our two ankylosaur taxa were conferring a heat-transfer benefit, we digitally manipulated duplicates of our ST airways to remove either the length or the convolutions from the nasal vestibule (Figs 17 and 18). One version had a nasal vestibule that extended the straight-line distance between the nostril and the CNP. This model, referred to as the basic airway (Figs 17B and 18B), was used to represent what a simplified or plesiomorphic nasal vestibule would look like. The second version of the ST airway retained the total length of the nasal vestibule but had the curvature of the airway removed (Figs 17A and 18A). This model was referred to as the straightened airway. It represented the effects of airway distance alone on heat transfer through the nasal passage.

Boundary conditions. Prior to volumetric meshing, the airway models were assigned a series of boundary conditions comprising a set of criteria that described this region of the model to the CFD program. Boundary condition assignment was done to elicit physiologically realistic airflow within the nasal passage (i.e., pressure driven air movement between nostril and choana). These conditions consisted of a pressure inlet located at the fleshy nostril and a pressure outlet located at the end of the artificial trachea (Fig 27A). During expiration, the assignment of these boundary conditions (inlet and outlet) was swapped. A series of impermeable wall boundaries covered the rest of the nasal passage model. Wall boundaries were demarcated based on anatomical location (Fig 27A). This was done to better control for regional variation in heat transfer across the nasal passage. Each wall boundary was considered rigid and incorporated a “no-slip” condition that states that air at the fluid-solid interface would be static, an assumption based on known properties of fluid movement through enclosed structures [71]. Note that amniote nasal passages do not have a truly rigid boundary layer between the mucosa and the air field. Boundary layers act as obstacles to diffusion-based processes, thus it is beneficial for amniotes to have means of reducing the size of these boundary layers. In extant amniotes, there is a mucociliary “conveyor belt” comprised of ciliated epithelium that beats in unison towards the nostril or choana [127,128]. This conveyor belt serves to move mucous across the mucosa of the nasal passage. This movement has the potential to reduce the boundary layer between the mucosa and the air field, which would aid in diffusion of heat and odorant molecules across the air-mucosa boundary. However, the speed of cilial movement is extremely slow (≤ 1 cm/min, [45]) compared to airflow, and its effects on airflow and heat transfer can be considered negligible for the purposes of our study. Thus, our use of a no-slip boundary condition should not hamper or otherwise reduce the quality of our results. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 27. Mesh example for Panoplosaurus mirus (ROM 1215). (A) Nasal passage was assigned a series of boundary conditions (color-coded). Black line indicates location of (B) axial cross section illustrating the distribution of volumetric cells within the nasal passage. https://doi.org/10.1371/journal.pone.0207381.g027

Meshing. Volumetric meshing was performed using the meshing program ICEM CFD (ANSYS Inc., Canonsburg, PA). Models consisted of a tetrahedral-hexahedral (tet-hex) hybrid mesh. First, an unstructured tetrahedral mesh was constructed using the robust OCTREE method [129]. After mesh reconstruction, the core of the mesh was deleted and flooded with hexahedra wherever possible. This hybrid construction offered the versatility of tetrahedra for unstructured mesh reconstruction [130,131], coupled with the computational efficiency of hexahedra [132,133]. To better resolve wall boundary effects on heat transfer, we incorporated a prism layer along the wall of the nasal passage. This layer consisted of four cells that grew in size from the wall by 1.5 times their parent cell, producing a combined thickness of approximately 0.7 mm (Fig 27B).

Computational fluid dynamic analysis. Meshes were imported into the CFD program Fluent (ANSYS Inc., Canonsburg, PA) for analysis. To determine the appropriate fluid dynamic model to apply to our airway models, we first took cross sections throughout the nasal passages to determine the dimensionless Reynolds and Womersley numbers for the airway. The Reynolds number is a staple of fluid dynamic analyses [71,134], representing the mathematical relationship between the viscous and inertial forces within a fluid. Low Reynolds numbers (< 2000) indicate that viscous forces dominate the system and that orderly (i.e, laminar) flow will be the dominant flow type expected. Reynolds numbers between 2000–4000 indicate a transition zone in which laminar flow may be punctuated with periods of turbulence that manifest in the form of secondary flows such as Von Kármán trails [71]. Reynolds numbers above 4000 indicate that inertial forces dominate the system and that flow will be chaotic or turbulent [71,134]. We used the following equation to calculate the Reynolds number for the airway [135]: (3) where Q = volumetric flow rate (m3/sec), P = the wetted perimeter in meters [136], and v = the kinematic viscosity of air at 15°C (1.412e-5 m2/sec). We used the dimensionless Womersley number [137] to determine the steadiness of the flow field, or how often airflow was able to produce a steady, parabolic flow profile. This profile is related to the size of the airway, viscosity of the fluid, and frequency of the oscillation (i.e., breathing rate). Womersley numbers < 1 indicate a quasi-steady flow field that can be modeled as time independent, or steady-state. As the Womersley number climbs above unity, steadiness decreases. At Womersley numbers > 10 the oscillation of the flow is too high for flow to completely develop [138] and is considered unsteady, thus requiring a transient or time-step-based modeling approach. We calculated the Womersley number using the following equation [45]: (4) where f = the frequency of oscillation (Hz), and Dh = hydraulic diameter of the airway. Hydraulic diameter was calculated using the following equation [136]: (5) where A = the area of the cross section measured (m2).

Physiological variables. Following our previous methodology [80], we used the same phylogenetically corrected allometric equation for resting respiration in birds [49] to estimate flow rates during inspiration and expiration (Table 8). Mass estimates for both ankylosaurs were taken from the supplementary data in Brown et al. [139]. Estimated masses for the two ankylosaurs were over 20 times larger than the largest animal in the Frappell et al. dataset ([49], an 88 kg ostrich). Using the equations from Frappell et al. [49] required extending the regression line well beyond the initial scope of the data. To alleviate the effects of this approach we used the lightest mass estimates for Panoplosaurus and Euoplocephalus as provided by Brown et al. ([139], Table 8). PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Table 8. Respiratory values for the ankylosaurs Panoplosaurus mirus and Euoplocephalus tutus. https://doi.org/10.1371/journal.pone.0207381.t008 We performed a cross-sectional analysis of the airways following the methodology of Bourke et al. [80]. Respiratory flow variables input into the Reynolds and Womersley equations above indicated that both taxa would have had steady flow, but air would have been largely transitional if not completely turbulent within the nasal passages during simulated respiration. The latter result was at odds with previously published results on resting respiration in amniotes. In extant amniotes measured during restful breathing, laminar flow dominates the air field [48,140–142]. Laminar flow is less energetically expensive (due to lower resistance) than turbulent flow and can be expected in animals that are not undergoing strenuous exercise. Our two ankylosaurs did not show this pattern, suggesting that the regression equations we used may not be viable at such large body masses or that ankylosaur body plans were not appropriate for the avian-based equation of Frappell et al. [49]. To reduce the effects of poor flow rate choice, we recalculated airflow rates by rearranging the Reynolds equation to solve for flow rate: (6) We refer to this new method as the Reversed-Reynolds approach. This approach provides an upper limit to volumetric flow through the nasal passage under relaxed conditions. We set the Reynolds number at 2000 (corresponding to the transition into turbulence) as our constant and recalculated flow rate across the nasal passages. The lowest flow rate obtained from cross-sectional analysis using this Reversed-Reynolds equation was chosen as the new flow rate (Table 8). This approach, the first of its kind as far as we know, ensures that laminar flow dominates the air field under simulated respiration. Flow rate has been shown to be the primary variable affecting the efficiency of the nasal passage at managing heat flow [29,70]. With this in mind, we chose to run the ST and BB models under both flow rate estimates. The initial flow rate, estimated from the regression by Frappell et al. [49], was deemed the “high flow rate” condition. Our revised flow rate estimate, calculated using our Reverse-Reynolds method, was deemed the “low flow rate” condition. For the high flow rate condition, we used the Wilcox two-equation κ–ω turbulence model [143] with low Reynolds corrections and shear stress transport. For the low flow rate condition, we used the standard laminar viscosity model for continuity, momentum, and energy. (7) (8) (9) Where = the velocity vector, ρ = the density of air at a given temperature, T = T(x,y,z) temperature at a given time, Cp = the specific heat, and k = the thermal conductivity for air at a given temperature.

Environmental conditions. We simulated sea-level air at 15°C and 50% relative humidity (r.h.), which was within the range of expected conditions that these dinosaurs would have experienced in their environment [144,145]. We gave the nasal passages an estimated body temperature of 35°C. This body temperature fell within the range of core body temperatures typically observed in extant large terrestrial mammals, birds, and reptiles (Table 5). Density (1102 kg/m3), thermal conductivity (0.34 W/mK), and specific heat for the mucosal walls (3150 J/kgK) were obtained from the database provided by the Foundation for Research on Information Technologies in Society (IT’IS). Nasal walls were given a thickness of approximately 0.5 mm for heat to conduct through. This distance simulated the distance between the nasal passage and the adjacent blood vessels and was based on observed CT data from extant amniotes [146]. Humidity was simulated by using the species transport option in Fluent, which allows for the incorporation of the mass fractions of water at various temperatures (i.e., relative humidity). The nasal walls were assumed to be at 100% relative humidity during both phases of respiration, following observation on extant animals [147]. Pressure and velocity coupling used the SIMPLEC algorithm along with a node-based discretization gradient. We used a second order accurate spatial discretization scheme for pressure, momentum, turbulence (when applicable), and energy. Models ran until the results obtained from each analysis had reached a specified level of stability and consistency referred to as convergence [148]. This indicated that the numerical process used to solve the problem had asymptotically approached the “true” solution (i.e., airflow and heat transfer through a real nasal passage) given the conditions provided to the program. In CFD, convergence may be determined based on the global imbalances in the values for each node within the mesh between steps, or iterations, of the model. Imbalances, or errors, between each iteration are referred to as residuals [148]. The smaller these residuals are, the smaller the error is and the more converged the solution becomes. Global variables for momentum, pressure, and continuity (the conservation of fluid mass) are generally considered solved when their residuals have fallen below 1.0e-3 [148]. However, for physiological studies such as ours, we applied the stricter criterion of 1.0e-4 [45,80]. For energy (heat flow), convergence is determined when the residuals of error had fallen below 1.0e-6. To further aid in determining convergence we monitored special surfaces placed throughout regions of the model. These surfaces were designed to output data from a single location only (point surfaces). They provided localized measures of convergence and were used in conjunction with standard convergence measures for continuity, momentum, and energy to determine when models had been sufficiently solved.

Mesh independence. To ensure that the results obtained from our analyses were independent of the mesh resolution, we performed a solution-based adaptive mesh refinement (AMR), following previously established protocols [82]. This approach used converged or a mostly converged solution to determine regions of the mesh that had poor resolution. Local refinement was performed in these regions of poor resolution and the analysis was run again. This approach was repeated until refined meshes return values that fell below a pre-determined threshold. For our analysis, we used a threshold of 1% for regional temperature values as determined from point surfaces (Fig 28). PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 28. Final model resolutions used for simulations. Surrounding graphs show sample adaptive mesh comparisons between different model resolutions for temperature and velocity within parts of the nasal vestibule. https://doi.org/10.1371/journal.pone.0207381.g028

Post processing and heat flow measurements. Solved models were exported to the CFD module of Avizo (Avizo Wind) where qualitative and quantitative measurements were taken. For heat flow, we took measurements from cross sections of the nasal passage. These cross sections were taken orthogonal to the curvature of the nasal passage. Multiple measurements were taken from each cross section and the mean was recorded for each cross section.

Caloric costs and savings. Calculating the energetic costs of heating a single bolus of air by 20°C requires knowledge of the mass of air being moved between breaths, coupled with the caloric cost of heating that bolus of air. We calculated our estimated energetic costs using the following equation: (10) where M air = the mass of air in a given tidal volume (g), and Cp = the specific heat capacity of air, which is 0.24 cal/g°C across most physiological temperatures [147]. ΔT is the temperature change the air volume goes through (°C). The mass of air present in a single breath was determined by multiplying the tidal volume of air respired, by its density at 35°C, as shown in the following equation: (11) where V T = the tidal volume (L) and 1.146 = the density of air (g/L) at the estimated body temperature of 35°C. Body-temperature air density will be the limiting factor behind tidal volume within the lungs. To determine tidal volume, we used the equation relating tidal volume to body mass (M) in birds [49]. (12) Much of the heat loss from the nasal passage occurs via evaporation of water off the nasal mucosa [147]. Thus, along with the caloric cost of temperature change within the nasal passage (sensible heat), one must also take into account the caloric cost associated with the phase change of water from a liquid to a gas (latent heat). We used saturated steam table values to determine the latent heat of vaporization for our temperature range. The caloric cost of heating air by 20°C was calculated using the following equation: (13) where ΔM H2O = the absolute difference in the mass of water (kg) at two given humidities. ΔH vap = the latent heat of vaporization for a given temperature (cal/kg). The mass of water was determined by multiplying the mass fraction of water at a given humidity (g/kg), by the mass of air (kg) in a single tidal volume. As heat is a form of energy, the same equations for calculating the sensible heat gain to the air can be used to determine heat loss during expiration. Similarly, the caloric costs associated with the latent heat of vaporization will be the same values as the latent heat of condensation, just with a reversed sign. Thus, the equations used to determine the caloric costs of heating air during inspiration can be used to determine heat savings upon air cooling during expiration. The only change during expiration is the value for ΔH at the expired air temperature (cooler air holds less water), and the assumption that air leaves the nostril at 100% relative humidity regardless of expired air temperature [147]. Energy savings were calculated by taking the calories returned to the body during air cooling and condensation during expiration and dividing it by the initial caloric cost of heating and humidifying the air during inspiration. Water saving were calculated following the method of Schmidt-Nielsen et al. [149].